Green's function for reversible geminate reaction ... - Semantic Scholar
Green's function for reversible geminate reaction with volume reactivity Svetlana S. Khokhlova and Noam Agmon Citation: J. Chem. Phys. 137, 184103 (2012); doi: 10.1063/1.4764357 View online: http://dx.doi.org/10.1063/1.4764357 View Table of Contents: http://jcp.aip.org/resource/1/JCPSA6/v137/i18 Published by the American Institute of Physics.
Additional information on J. Chem. Phys. Journal Homepage: http://jcp.aip.org/ Journal Information: http://jcp.aip.org/about/about_the_journal Top downloads: http://jcp.aip.org/features/most_downloaded Information for Authors: http://jcp.aip.org/authors
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
THE JOURNAL OF CHEMICAL PHYSICS 137, 184103 (2012)
Green’s function for reversible geminate reaction with volume reactivity Svetlana S. Khokhlova and Noam Agmon The Fritz Haber Research Center, Institute of Chemistry, The Hebrew University of Jerusalem, Jerusalem 91904, Israel
(Received 27 September 2012; accepted 15 October 2012; published online 12 November 2012) The kinetics of a diffusing particle near a reversible trap may be described by an extension of the Feynman-Kac equation to the case of reversible binding, which can occur within a finite reaction sphere. We obtain the Green’s function solution for the Laplace transform of this equation when the particle is initially either bound or unbound. We study the solution in the time-domain by either inverting the Laplace transform numerically or propagating the partial differential equation in the time-domain. We show that integrals of this solution over the reaction sphere agree with previously obtained solutions. © 2012 American Institute of Physics. [http://dx.doi.org/10.1063/1.4764357] I. INTRODUCTION
Chemical reactions in solution are controlled by translational diffusion, which brings the reacting partners into a critical distance from which reaction may occur.1 A fundamental model for reaction involves “contact reactivity,” in which a pair of particles has first to collide in order to react. The reaction is then imposed as a boundary condition at the contact distance σ .2 This model can be generalized to a reversible diffusion-influenced reaction3 kr
A + B C,
(1.1)
kd
in which kr and kd are recombination and dissociation rate coefficients, respectively. In the case of a geminate pair, when we deal with a single C molecule or an isolated single pair of A and B, the solution can be obtained in Laplace space and inverted into the time-domain exactly, approximately, or numerically under a variety of conditions.4–23 This model has been applied extensively to (ultra-fast) excited-state proton transfer reactions in solution,24–40 within cavities,41–43 in ice,44, 45 and within the green fluorescent protein.46, 47 It has also been utilized for analyzing translational diffusion correlation functions in molecular simulations.48, 49 Despite this success, it is worthwhile to consider an alternate model for reversible geminate recombination. As opposed to the “contact reactivity” model, where reaction occurs at a specified A-B distance, r = σ , in the “volume reactivity” model reaction can occur throughout the spherical volume.50 In this model a particle diffuses freely in space, but when inside a “reaction sphere,” r ≤ σ , it can react reversibly with the rate coefficients kr and kd . The volume reactivity model can have different implementations with respect to the distance at which the pair A–B is released upon dissociation. Reference 51 considered the case when the release is at a specified distance, either smaller or larger than σ , and solved the ensuing diffusionreaction equations for the steady-state concentrations. Reference 50 considered the cases that the released particle is uniformly distributed within the reaction sphere or else when it appears at the same distance where it was trapped. This 0021-9606/2012/137(18)/184103/8/$30.00
immobilization/remobilization problem is the one treated here. The problem has been previously studied in Laplace space for a particle initially uniformly bound within the reaction sphere.50 Here we extend the solution to the Green’s function52 (GF) p(r,t|r0 ), namely, the probability that the particle transits from an initial distance r0 to r in threedimensional space during time t, in the presence of a reactive sphere of radius σ . Interestingly enough, the integral of the GF over r for kd = 0 obeys the celebrated Feynman-Kac equation (FKE).53–55 Its solution gives the “residence time” of the particle within the reaction sphere,56 called “occupation time” in random walk theory.57 The present work generalizes the solution, which can be found in the appendix of Ref. 55, to reversible binding and an initial δ-function distribution. This solution, which can only be obtained in Laplace space, is inverted numerically to the time-domain in order to exemplify its behavior for the various initial conditions. II. THE REVERSIBLE FEYNMAN-KAC EQUATION
In the “volume reactivity model,”50 a particle diffusing (with diffusion coefficient D) in three-dimensional space may bind and unbind only within a sphere of radius r = σ (centered at the origin), whose volume is V = 4π σ 3 /3. Thus we limit the discussion to spherical symmetry, where r is the norm of the radius vector. Let p(r, t) denote the probability density to find the particle unbound and at distance r from the origin at time t, while its probability density for being bound at some r < σ is denoted by q(r, t). Unbound and bound populations interchange at rates kr H(σ − r)p(r, t) and kd q(r, t) in the association and dissociation directions, respectively. H(x) is the Heaviside step-function, which is zero for x < 0 and 1 for x ≥ 0. It restricts the association reaction to occur only within the σ -sphere. No such restriction is needed in the dissociation direction because an initially unbound particle cannot bind at r > σ . The dissociating particle always resumes its diffusive motion at the same point to which it was bound. Considered as a model for reversible geminate recombination, here one particle (e.g., a ligand) can diffuse inside the
137, 184103-1
© 2012 American Institute of Physics
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-2
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
volume of the other (e.g., its binding site). For reaction between two small molecules this may seem artificial, yet of little quantitative consequence, because the volume is small. In other cases the present model may be more realistic than the contact reactivity model, for example, reversible binding of a diffusing particle to a circular receptor array.58 Using these notations, the following two coupled differential equations, a partial differential equation (PDE) plus an ordinary one, describe the reversible reaction in the volume reactivity model: ∂p(r, t) = L p(r, t) − kr H (σ − r) p(r, t) + kd q(r, t), ∂t (2.1a) ∂q(r, t) = kr H (σ − r) p(r, t) − kd q(r, t). ∂t
∂ 1 ∂ D(r)r 2 , 2 r ∂r ∂r
(2.2)
D(r) is the particle’s diffusion coefficient (for a pair it is their relative diffusion coefficient), which is hereafter assumed to be independent of r. In the irreversible limit there is no dissociation, so that kd = 0. Then the equations are no longer coupled: it suffices to solve Eq. (2.1a) for p(r, t) whereas q(r, t) could then be obtained by integrating the outgoing flux 4π r2 D∂p(r, t)/∂r with respect to t (for any r ≤ σ ). In this limit, Eq. (2.1a) is known as the Feynman-Kac equation. The boundary conditions (BCs) imposed on it are zero flux at the origin 2 ∂p(r, t) →0 (2.3) r ∂r r→0 and vanishing probability at infinite distances p(r, t) → 0 as
r → ∞.
(2.4)
There is no need to impose BCs on q(r, t) because in the bound state the particle does not diffuse. We solve these equations in Laplace space. The Laplace transformation (LT) of a function f(t), denoted by fˆ(s), is ∞ ˆ f (s) = f (t)e−st dt. 0
Taking the LT of Eqs. (2.1) we obtain ˆ s) = L p(r, ˆ s) − kr H (σ − r) p(r, ˆ s) s p(r, ˆ s) + p(r, 0), + kd q(r, ˆ s) = q(r,
For a particle that is initially unbound and located at distance r0 from the origin, we assume δ(r − r0 ) , 4π r02
p(r, 0|r0 ) =
q(r, 0|r0 ) = 0.
(2.5a)
1 ˆ s) + q(r, 0)]. [kr H (σ − r) p(r, s + kd (2.5b)
Here we still need to insert the initial conditions p(r, 0) and q(r, 0). For the GF, these involve a particle that is either unbound or bound at distance r0 .
(2.6)
By convention, the initial separation is denoted on the right of the vertical line, so that p(r, t|r0 ) depicts the GF in the timeˆ s|r0 ) is the GF in Laplace space. Insertdomain, whereas p(r, ing Eqs. (2.5b) and (2.6) into Eq. (2.5a) gives ˆ s|r0 ) − (s − L)p(r,
δ(r − r0 ) skr H (σ − r) ˆ s|r0 ), p(r, =− 2 s + kd 4π r0 (2.7)
(2.1b)
Here L is the spherically-symmetric diffusion operator in a three-dimensional space L≡
A. Initially unbound particle
which is an ordinary (because s is just a parameter) inhomogeneous differential equation in the variable r. This reversible FKE extends the conventional FKE in two respects. First, it is reversible, though this only amounts to replacing kr by skr /(s + kd ). Second, it applies to a δ-function initial condition, from which the solution for any other initial distribution can be obtained by averaging the GF over this distribution. The spatial integral of p(r, t|r0 ) over r is the “survival probability” of the unbound particle ∞ P (t|r0 ) ≡ 4π p(r, t|r0 )r 2 dr. (2.8) 0
The customary FKE deals with this function,54, 55 which is discussed in Appendix A, whereas Appendix B deals with a uniform initial distribution within the reaction sphere.50 It is convenient to solve Eq. (2.7) for r = r0 when it becomes homogeneous and then match the solutions at r0 . Denote the solutions for r < r0 and r > r0 by pˆ − (r, s|r0 ) and pˆ + (r, s|r0 ), respectively. They must be continuous at r0 , but their derivative is discontinuous there, pˆ − (r, s|r0 ) = pˆ + (r, s|r0 ),
(2.9a)
1 ∂ pˆ − (r, s|r0 ) ∂ pˆ + (r, s|r0 ) = − . (2.9b) ∂r ∂r 4π Dr02 r=r0 − r=r0 + Because Eq. (2.7) is of second order, each of these solutions is a linear combination of two fundamental solutions (see below) with coefficients determinable from the above condition. Due to the discontinuous Heaviside function, there is another matching sphere at r = σ , but there one demands that both the function and its derivative are continuous. Thus we solve the equation separately in three regions of space bounded by the spheres r = r0 and r = σ and match the solutions at these two interfaces. Finally, after solving for ˆ s|r0 ), one can obtain the probability density of the bound p(r, particle from Eq. (2.5b), namely ˆ s|r0 ) = q(r,
kr ˆ s|r0 ) , r ≤ σ. p(r, s + kd
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
(2.10)
184103-3
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
B. Initially bound particle
A. Unbound particle, initially outside: r0 > σ
For a particle that is initially bound and at r0 ≤ σ we impose
For r0 > σ , we consider the solution in the three regions: (1) r < σ , (2) σ < r < r0 , and (3) r > r0 , which are depicted in Fig. 1. In region 3 the BC at infinity eliminates the rising exponential, so that
p(r, 0|r0 ) = 0 ,
q(r, 0|r0 ) =
δ(r − r0 ) . 4π r02
(2.11)
ˆ s|r0 ) = p(r,
Instead of Eq. (2.7), we now obtain ˆ s|r0 ) − (s − L)p(r,
kd δ(r − r0 ) s + kd 4π r02
which has a modified δ-function term. Hence its solution is simply that of Eq. (2.7) times kd /(s + kd ). The probability density of the bound particle is subsequently obtained from Eq. (2.5b), namely, 1 δ(r − r0 ) ˆ s|r0 ) + ˆ s|r0 ) = kr p(r, . (2.13) q(r, s + kd 4π r02
ˆ s|r0 ) = p(r,
ˆ s|r0 ) = p(r,
Given the differential operator in Eq. (2.2), the two fundamental solutions of this equation are exp (± αr)/r, where (3.2) α ≡ s/D.
A sinh(βr), r < σ. r
(3.5c)
1 e−α(r0 −σ ) [α sinh(βσ ) + β cosh(βσ )]−1 , 4π Dr0 (3.6a)
B1 =
1 e−αr0 , 8π Dαr0
B2 =
α sinh(βσ ) − β cosh(βσ ) −α(r0 −2σ ) 1 , e 8π Dαr0 α sinh(βσ ) + β cosh(βσ )
(3.6b)
(3.6c) (3.3)
and the two fundamental solutions are exp (± βr)/r, where β is defined by β = α[1 + kr /(s + kd )]1/2 .
(3.5b)
The coefficients A, B1 , B2 , and C are functions of s and r0 (but not of r), to be determined from the matching conditions at r = r0 and r = σ . Thus we find that
Let us consider Eqs. (2.7) and (2.12) for r = r0 and either r > σ or r < σ . For r > σ they simplify to (3.1)
B1 αr B2 −αr e + e , σ < r < r0 . r r
Finally, in region 1 the vanishing flux at the origin dictates that
A=
III. THE GREEN’S FUNCTION SOLUTION
For r < σ they become kr ˆ s|r0 ) = 0 p(r, L−s 1+ s + kd
(3.5a)
In region 2 both exponentials contribute,
skr H (σ − r) ˆ s|r0 ), (2.12) p(r, =− s + kd
ˆ s|r0 ) = 0. (L − s)p(r,
C exp(−αr), r > r0 . r
(3.4)
When kr = 0, β = α and Eq. (2.7) reduces to free diffusion. Let us now solve these equations separately for r0 > σ or r0 < σ (Figure 1).
C = B2 +
1 eαr0 . 8π Dαr0
(3.6d)
These coefficients obey the following four identities: (i) From continuity of the solution and its derivative jump at r = r0 : Ce−αr0 = B1 eαr0 + B2 e−αr0 , Ce−αr0 =
1 − B1 eαr0 + B2 e−αr0 . 4π αDr0
(3.7a) (3.7b)
(ii) From continuity of the solution and its first derivative at r = σ:
FIG. 1. The spatial arrangement of elements pertaining to the GF solution, which depends on the two distances r and r0 . The initial location of the particle, r = r0 , can be either inside or outside the reactive sphere, r = σ . This divides space, in two different ways, into three regions (indicated by 1, 2, and 3). The particle can react only within the colored sphere, where r ≤ σ .
A sinh(βσ ) = B1 eασ + B2 e−ασ ,
(3.8a)
Aβ cosh(βσ ) = B1 αeασ − B2 αe−ασ .
(3.8b)
B. Unbound particle, initially inside: r0 < σ
For r0 < σ we consider the solution in the three regions: (1) r < r0 , (2) r0 < r < σ , (3) r > σ (see Fig. 1). The solutions are the same as in Eqs. (3.5), except for region 2, where α is
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-4
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
replaced by β, ˆ s|r0 ) = p(r,
C exp(−αr), r
r > σ,
(3.9a)
ˆ s|r0 ) = p(r,
B1 −β(σ −r) B2 β(σ −r) e e + , r r r0 < r < σ,
(3.9b)
ˆ s|r0 ) = p(r,
A sinh(βr), r
r < r0 .
(3.9c)
The matching conditions now give A=
1 α sinh[β(σ − r0 )] + β cosh[β(σ − r0 )] , 4π Dβr0 α sinh(βσ ) + β cosh(βσ ) (3.10a)
B1 =
β − α −ασ e C, 2β
(3.10b)
B2 =
β + α −ασ e C, 2β
(3.10c)
eασ sinh(βr0 ) . 4π Dr0 α sinh(βσ ) + β cosh(βσ )
(3.10d)
C=
These coefficients obey the following four identities: (i) From continuity of the solution and its first derivative at r = σ: Ce−ασ = B1 + B2 ,
(3.11a)
αCe−ασ = β(B2 − B1 ).
(3.11b)
(ii) From continuity of the solution and its derivative jump at r = r0 : A sinh(βr0 ) = B1 e−β(σ −r0 ) + B2 eβ(σ −r0 ) , A cosh(βr0 ) =
(3.12a)
1 + B1 e−β(σ −r0 ) − B2 eβ(σ −r0 ) . 4πβDr0
FIG. 2. The distance dependence of the Laplace transformed GF for an initially unbound particle at three values of r0 : inside (r0 = 0.5), on the surface (r0 = 1), or outside (r0 = 2) the reaction sphere. Lines demonstrate Eqs. (3.5) and (3.6), or (3.9) and (3.10), for s = 1 (full line), 5 (dashed line), and 10 (dashed–dot line). All the figures in this publication use the same parameters set: σ = 1, D = 1, kd = 0.5, and kr = 3.
Figure 2 shows the r-dependence of the solution in ˆ Laplace space, p(r,s|r 0 ), for a particle that is initially unbound at r0 . The solution for r0 inside or outside the reaction sphere is given by Eqs. (3.5) and (3.6) or (3.9) and (3.10), ˆ respectively. For large s, the function p(r,s|r 0 ) peaks at r = r0 and spreads over r as s decreases. The derivative discontinuity at r = r0 is evident in the figure, whereas there are no discontinuities at r = σ . Figures 3–5 show the GF solutions in the time-domain, for the same three initial locations. The black lines are ˆ ˆ the (numerical) Laplace inverse of p(r,s|r 0 ) and q(r,s|r 0 ). The colored dotted lines are the exact numerical solution of the PDE in the time-domain using SSDP (version 2.66).60 There is an exact match between the two solutions. The GF for an unbound particle, p(r,t|r0 ) (red), looks initially like a widening Gaussian. Its peak moves to smaller
(3.12b) inverse LT SSDP p(r,t|r0)
t1 = 0.04663
The solution of Eq. (2.12) is the same as in Eqs. (3.10) above, with A and C multiplied by kd /(s + kd ). IV. DEMONSTRATIVE CALCULATIONS
In this section we demonstrate the above GF solution in both Laplace space and the time-domain. The latter is achieved by a numerical LT inversion using the algorithm of Ref. 59. It is compared with the numerical solution of Eq. (2.1) in the time-domain using the SSDP software.60 The agreement validates our analytical LTs (or, alternately, it validates the numerical procedures for LT inversion and PDE solution). The parameters used in the following figures are the same as in our previous publication:50 σ = 1, D = 1, kd = 0.5, and kr = 3D/σ 2 = 3.
t2 = 0.2489
p (r, t | r0 ), q (r, t | r0 )
C. Bound particle, initially inside: r0 < σ
0.02
1
SSDP q(r,t|r0)
t3 = 0.5749
0.01
2
3
3 2 0.00
0
1
2
3
4
r FIG. 3. The distance dependence of the GF in the time-domain, for an initially unbound particle starting outside the reaction sphere, at r0 = 2 > σ . Black lines demonstrate the Laplace inverse of Eqs. (3.5) and (3.6) for p(r, t|r0 ), as well as that of Eq. (2.10) for the probability density of the bound state, q(r, t|r0 ).
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-5
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012) 0.10
t4 = 4.659 t5 = 6.881 t6 = 11.7 1
q ( r, t | r0 )
p ( r, t | r0 )
t3 = 1.469
0.02
0.01
0.08
t2 = 0.5185
3
4
4
t1 = 0.02102
0.03 2
3 5
0.06
0.04
2 6
5
0.02
6 1 0.00 0.00 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.2 0.4 0.6 0.8 1.0 1.2
r FIG. 4. Same as Fig. 3 for r0 = 1 = σ .
values of r with increasing time due to reflection from r = 0. In contrast to these smooth profiles, the GF for the bound particle, q(r,t|r0 ) (green), has a cusp. This cusp is on the boundary of the reaction sphere (r = σ ) when r0 ≥ σ and at r0 when r0 < σ . Due to the recombination reaction, part of the initial δ-function density for the unbound particle becomes a δfunction for the bound particle. This population cannot diffuse in the bound state, only dissociate, giving rise to a derivative discontinuity in the time-domain. The time evolution of the GF for an initially bound particle (at r0 ≤ σ ) is shown in Fig. 6. The initial distribution is a δ-function q(r, 0|r0 ) = δ(r − r0 )/4π r02 (at r0 = 0.5) that persists for all times. Its amplitude, however, diminishes with time (not shown) as the particle dissociates. This creates the distribution p(r, t|r0 ) of the unbound particle, which is initially close to a δ-function (time t1 ), then widens (due to diffusion) and increases in amplitude (due to further dissociation – time t2 ), and finally (t3 onwards) decreases in amplitude due to escape of the particle to infinity. Recombination from the unbound state then creates the non-δ-function part of q(r, t|r0 ), its amplitude following that of p(r, t|r0 ) with a delay (due to the finite kr ). Thus the amplitude of q(r, t|r0 ) continues to in0.7
1 inverse LT SSDP p(r,t|r0)
0.6
2
p (r, t | r0 ), q ( r, t| r0 )
0.5
0.4
SSDP q(r,t|r0)
3
0.3
0.2
0.1
0.0
0.0
3 1 0.5
1.0
1.5
r FIG. 5. Same as Fig. 3 for r0 = 0.5 < σ , except that the times here are: t1 = 0.02019, t2 = 0.0307, and t3 = 0.05361. Black lines were calculated by Laplace inverting Eqs. (3.9) and (3.10).
r
FIG. 6. The time dependence of the GF for a particle that is initially bound ˆ s|r0 ) is given by at r0 = 0.5. The LT of the GF for the unbound state p(r, Eq. (3.9) multiplied by kd /(s + kd ) and its Laplace inverse is depicted by the ˆ s|r0 ), is calculated from blue lines. The LT of the GF for the bound state, q(r, ˆ s|r0 ) using Eq. (2.13) and its Laplace inverse is depicted by the red lines. p(r, The numerical solution of the PDE (2.1) for these initial conditions is shown by the dotted green lines.
crease up to about t4 , when p(r, t|r0 ) is already in its decreasing phase. Another interesting observation is that both parts of the GF show a derivative discontinuity at r0 . This contrasts with the initially unbound particle, for which only q(r, t|r0 ) shows such a discontinuity (Figs. 3–5). V. CONCLUSION
In this work we have obtained (analytically) the GF solution of the reversible FKE in Laplace space and (numerically) in the time-domain. The relevance of this equation to chemical kinetics is in providing a mathematical description for the “volume reactivity” model of geminate recombination, in which one particle is a reactive sphere and the other is a diffusing point particle that can bind (immobilize) and unbind (remobilize) within this sphere. This model is an alternative to the widely used “contact reactivity” model, in which the two particles can react only upon collision (at the “contact distance”). Interestingly, the GF solution for reversible binding of an initially unbound particle can be obtained from that of irreversible binding (kd = 0) by replacing kr /s by kr /(s + kd ). An initially bound particle is immobilized for all times in the irreversible case, yet its reversible GF can also be obtained from the solution of the irreversible case by a further multiplicative factor kd /(s + kd ). In Laplace space, all the components of the GF necessarily show a derivative discontinuity at the initial location of the particle, r0 . In the time-domain, p(r, t|r0 ) does not show such a discontinuity for an initially unbound particle, but all other GF components do show a derivative discontinuity at r0 , which is conceivably due to the lack of diffusion in the bound state. Appendices A and B show that the integrals of our GF’s (over r or r0 ) reproduce the solutions for the overall survival probability of an initially unbound particle (obtained in Ref. 55 for the irreversible case) and for an initially uniformly
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-6
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
distributed bound particle.50 Thus all previous solutions of the FKE are special cases of the GF obtained in the present work. The question is what is the practical difference between the contact and volume reactivity models and is a volume reactivity model really needed? As Fig. 5 in reference 50 shows, for an initially uniformly distributed particle within the reaction sphere the difference in the binding probability, Q(t|uni), between the two models may be rather subtle (see Fig. 5). However, the difference between the models may be larger for an initially unbound particle or for the full GF. This issue still needs to be studied. More interesting physically is the applicability of the volume reactivity model beyond the pair problem. Consider, for example, the case where reversible receptors are scattered within the reaction sphere.58 For example, a particle may be diffusing on the surface of a membrane,61 on which there is a circular patch containing many reversible receptors. (Such a situation is believed to occur, e.g., for protons on membranes.62, 63 ) Then, in a mean field sense, the kinetics of the particle is described by a two-dimensional version of the reversible FKE solved herein.
D −1/2 = βK1/2 (ασ )I3/2 (βσ ) + αK3/2 (ασ )I1/2 (βσ ). (A3) Bessel functions of the half-integer order can be written in terms of hyperbolic functions [Eqs. (10.2.13) and (10.2.17) in Ref. 65]. Then in Eq. (A3) reduces to Dβ 1/2 [β cosh(βσ ) + α sinh(βσ )]e−ασ . βσ = α (A4) The LTed binding probability in Eq. (A2) then becomes 1 sinh(βr0 ) kr 1 + ασ ˆ , ) = − Q(s|r 0 s s + kr β 2 r0 D β cosh(βσ ) + α sinh(βσ ) r0 ≤ σ, (A5a) kr e−α(r0 −σ ) βσ cosh(βσ ) − sinh(βσ ) ˆ Q(s|r , r0 ≥ σ. 0) = s β 2 r0 D β cosh(βσ ) + α sinh(βσ )
ACKNOWLEDGMENTS
Research supported by The Israel Science Foundation (Grant No. 122/08). The Fritz Haber Research Center is supported by the Minerva Gesellschaft für die Forschung, München, FRG. APPENDIX A: THE INTEGRATED FKE
The FKE usually encountered in the literature is for the LT of the survival probability, Pˆ (s|r0 ), obtained by integrating the GF as in Eq. (2.8). It is a solution of the integrated adjoint of Eq. (2.7), namely, skr H (σ − r0 ) s+ − L0 Pˆ (s|r0 ) = 1, (A1) s + kd with the adjoint diffusion operator L0 and the same boundary conditions at r0 = 0 and ∞ as above. This solution, for kd = 0, is given in the appendix of Ref. 55 for any dimensionality d. See also Eq. (1.5.1) in Chap. 4 of Part II of Ref. ˆ 64. The binding probability Q(s|r 0 ) can be obtained from the ˆ normalization condition Q(s|r0 ) = 1/s − Pˆ (s|r0 ). It is actually slightly simpler than the survival probability, so we focus on it here. For the irreversible case (kd = 0), Eqs. (59) and (60) in Ref. 55 in three dimensions become
αK3/2 (ασ ) σ 1/2 1 kr ˆ − 2 1/2 I1/2 (βr0 ) , Q(s|r0 ) = s s + kr β D r0 r0 ≤ σ, kr I3/2 (βσ ) ˆ Q(s|r 0) = s β D 1/2
Here Iν (z) and Kν (z) are the modified Bessel functions of order ν, and in the denominator is given by Eq. (61) of Ref. 55:
σ r0
(A2a)
1/2 K1/2 (αr0 ) , r0 ≥ σ. (A2b)
(A5b) The solution for reversible binding, kd = 0, is obtained by replacing kr /s by kr /(s + kd ) in the above equations, whereas β generalizes to its value in Eq. (3.4). This gives 1 + ασ kr ˆ 1− Q(s|r ) = 0 s(s + kd + kr ) r0 sinh(βr0 ) , r0 ≤ σ, × (A6a) β cosh(βσ ) + α sinh(βσ ) ˆ Q(s|r 0) =
e−α(r0 −σ ) βσ cosh(βσ ) − sinh(βσ ) kr , s(s +kd +kr ) r0 β cosh(βσ ) + α sinh(βσ ) r0 ≥ σ. (A6b)
Equation (A6) can also be obtained by direct integration of the GF over r. According to Eq. (2.10), obtaining the binding probability requires integration only inside the reaction sphere σ 4π kr ˆ ˆ s|r0 )r 2 dr, p(r, Q(s|r ) = (A7) 0 s + kd 0 and this is simpler than obtaining the survival probability that requires integration over all space. When r0 < σ , one needs ˆ s|r0 ) in regions 2 and 3, see Eqs. (3.9b) and to integrate p(r, (3.9c), but only in region 1, Eq. (3.5c), when r0 > σ . The surˆ vival probability can subsequently be obtained from Q(s|r 0) ˆ using the normalization condition Pˆ (s|r0 ) = 1/s − Q(s|r 0 ). Performing these integrals directly indeed gives the same result as above. Its Laplace inverse, P(t|r0 ), is demonstrated in Fig. 7. It starts from unity, then decreases to a minimum as the particle binds. But eventually it is released and escapes to infinity when P(t|r0 ) approaches unity again.
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-7
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012) 5 N.
1.0
P (t | r0 )
0.9
r0 = 1
0.8
0.7
SSDP P(t | r0)
^ |r ) inverse LT of P(s 0
r0 = 0.5
0.6
0
5
10
15
t FIG. 7. The time dependence of the survival probability, P(t|r0 ), obtained by Laplace inverting Eqs. (A5), as compared to numerically integrating the numerical solution of Eq. (2.1) using the SSDP software. Parameters are the same as in Fig. 2.
APPENDIX B: AN INITIAL UNIFORMLY BOUND STATE
Instead of a particle initially bound at r0 , we consider a bound particle that is initially uniformly distributed inside the reaction sphere. Denoting the survival probability in this case by P(r, t|uni), we define σ 4π P (r, t|uni) ≡ p(r, t|r0 ) r02 dr0 , (B1) V 0 where the division by V = 4π σ 3 /3 serves to normalize the initial distribution. This problem was considered in Ref. 50. Because the diffusion equation is symmetric to the exchange of r and r0 , we find from Eq. (A7) that s + kd ˆ Pˆ (r, s|uni) = Q(s|r). V kr
(B2)
Inserting the binding probability from Eq. (A6) with r0 replaced by r and multiplying by kd /(s + kd ) to account for an initially bound state, we obtain 1 + ασ kd /V 1− Pˆ (r, s|uni) = s(s + kd + kr ) r sinh(βr) × , r ≤ σ, β cosh(βσ ) + α sinh(βσ ) (B3a) Pˆ (r, s|uni) =
kd /V e−α(r−σ ) s(s + kd + kr ) r ×
βσ cosh(βσ ) − sinh(βσ ) , r ≥ σ. (B3b) β cosh(βσ ) + α sinh(βσ )
This is identical to Eq. (4.14) in Ref. 50. 1 S.
Agmon and A. Szabo, J. Chem. Phys. 92, 5270 (1990). Nadler and D. L. Stein, J. Chem. Phys. 104, 1918 (1996). 7 B. C. Lagerholm and N. L. Thompson, Biophys. J. 74, 1215 (1998). 8 I. V. Gopich, K. M. Solntsev, and N. Agmon, J. Chem. Phys. 110, 2164 (1999). 9 N. Agmon, J. Chem. Phys. 110, 2175 (1999). 10 I. V. Gopich and N. Agmon, J. Chem. Phys. 110, 10433 (1999). 11 H. Kim, K. J. Shin, and N. Agmon, J. Chem. Phys. 111, 3791 (1999). 12 H. Kim and K. J. Shin, Phys. Rev. Lett. 82, 1578 (1999). 13 N. Agmon and I. V. Gopich, J. Chem. Phys. 112, 2863 (2000). 14 H. Kim and K. J. Shin, J. Chem. Phys. 112, 8312 (2000). 15 H. Kim, K. J. Shin, and N. Agmon, J. Chem. Phys. 114, 3905 (2001). 16 H. Kim and K. J. Shin, J. Chem. Phys. 120, 9142 (2004). 17 N. Agmon, Chem. Phys. Lett. 417, 530 (2006). 18 S. Y. Reigh, K. J. Shin, and M. Tachiya, J. Chem. Phys. 129, 234501 (2008). 19 S. Park and N. Agmon, J. Chem. Phys. 130, 074507 (2009). 20 K. Park, K. J. Shin, and H. Kim, J. Chem. Phys. 131, 154105 (2009). 21 H. Kim, Chem. Phys. Lett. 507, 265 (2011). 22 I. Plante and F. A. Cucinotta, Phys. Rev. E 84, 051920 (2011). 23 Z. Kalay, J. Phys. A 45, 235001 (2012). 24 E. Pines, D. Huppert, and N. Agmon, J. Chem. Phys. 88, 5620 (1988). 25 D. Huppert, E. Pines, and N. Agmon, J. Opt. Soc. Am. B 7, 1545 (1990). 26 A. Masad and D. Huppert, Chem. Phys. Lett. 180, 409 (1991). 27 N. Agmon, D. Huppert, A. Masad, and E. Pines, J. Phys. Chem. 95, 10407 (1991); Erratum 96, 2020 (1992). 28 S. Y. Goldberg, E. Pines, and D. Huppert, Chem. Phys. Lett. 192, 77 (1992). 29 N. Agmon, S. Y. Goldberg, and D. Huppert, J. Mol. Liq. 64, 161 (1995). 30 I. Carmeli, D. Huppert, L. M. Tolbert, and J. E. Haubrich, Chem. Phys. Lett. 260, 109 (1996). 31 D. Huppert, L. M. Tolbert, and S. Linares-Samaniego, J. Phys. Chem. A 101, 4602 (1997). 32 K. M. Solntsev, D. Huppert, N. Agmon, and L. M. Tolbert, J. Phys. Chem. A 104, 4658 (2000). 33 K. M. Solntsev, D. Huppert, and N. Agmon, Phys. Rev. Lett. 86, 3427 (2001). 34 N. Koifman, B. Cohen, and D. Huppert, J. Phys. Chem. A 106, 4336 (2002). 35 C. Clower, K. M. Solntsev, J. Kowalik, L. M. Tolbert, and D. Huppert, J. Phys. Chem. A 106, 3114 (2002). 36 K. M. Solntsev, E. N. Sullivan, L. M. Tolbert, S. Ashkenazi, P. Leiderman, and D. Huppert, J. Am. Chem. Soc. 126, 12701 (2004). 37 L. Genosar, P. Leiderman, N. Koifman, and D. Huppert, J. Phys. Chem. A 108, 1779 (2004). 38 N. Agmon, J. Phys. Chem. A 109, 13 (2005). 39 P. Leiderman, R. Gepshtein, A. Uritski, L. Genosar, and D. Huppert, J. Phys. Chem. A 110, 9039 (2006). 40 M. R. di Nunzio, B. Cohen, and A. Douhal, J. Phys. Chem. A 115, 5094 (2011). 41 B. Cohen, D. Huppert, K. M. Solntsev, Y. Tsfadia, E. Nachliel, and M. Gutman, J. Am. Chem. Soc. 124, 7539 (2002). 42 R. Gepshtein, P. Leiderman, D. Huppert, E. Project, E. Nachliel, and M. Gutman, J. Phys. Chem. B 110, 26354 (2006). 43 N. Amdursky, R. Orbach, E. Gazit, and D. Huppert, J. Phys. Chem. C 113, 19500 (2009). 44 P. Leiderman, A. Uritski, and D. Huppert, J. Phys. Chem. A 111, 4998 (2007). 45 A. Uritski, I. Presiado, and D. Huppert, J. Phys. Chem. C 112, 11991 (2008). 46 P. Leiderman, D. Huppert, and N. Agmon, Biophys. J. 90, 1009 (2006). 47 N. Agmon, J. Phys. Chem. B 111, 7870 (2007). 48 O. Markovitch and N. Agmon, J. Chem. Phys. 129, 084505 (2008). 49 H. Chen, G. A. Voth, and N. Agmon, J. Phys. Chem. B 114, 333 (2010). 50 S. S. Khokhlova and N. Agmon, Bull. Korean Chem. Soc. 33, 1020 (2012). 51 J. Lipková, K. C. Zygalakis, S. J. Chapman, and R. Erban, SIAM J. Appl. Math. 71, 714 (2011). 52 M. C. M. Wright, Nat. Phys. 2, 646 (2006). 53 M. Kac, in Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, edited by J. Neyman (University of California, Berkeley, 1951), pp. 189–215. 54 S. N. Majumdar, Curr. Sci. 89, 2076 (2005), arXiv:cond-mat/0510064. 55 N. Agmon, J. Phys. Chem. A 115, 5838 (2011). 56 N. Agmon, J. Chem. Phys. 81, 3644 (1984). 57 S. Karlin and H. M. Taylor, A Second Course in Stochastic Processes (Academic, San-Diego, 1981). 58 A. L. Edelstein and N. Agmon, J. Comput. Phys. 132, 260 (1997). 6 W.
r0 = 2
A. Rice, in Diffusion-Limited Reactions, Comprehensive Chemical Kinetics Vol. 25 (Elsevier, Amsterdam, 1985). 2 F. C. Collins and G. E. Kimball, J. Colloid Sci. 4, 425 (1949). 3 N. Agmon, J. Chem. Phys. 81, 2811 (1984). 4 N. Agmon and G. H. Weiss, J. Chem. Phys. 91, 6937 (1989).
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-8
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
59 J.
63 N.
60 E.
64 A.
Abate and W. Whitt, ORSA J. Comput. 7, 36 (1995). B. Krissinel’ and N. Agmon, J. Comput. Chem. 17, 1085 (1996). 61 N. Agmon, Chem. Phys. 370, 232 (2010). 62 A. Springer, V. Hagen, D. A. Cherepanov, Y. N. Antonenko, and P. Pohl, Proc. Natl. Acad. Sci. U.S.A. 108, 14461 (2011).
Agmon and M. Gutman, Nat. Chem. 3, 840 (2011). N. Borodin and P. Salminen, in Handbook of Brownian Motion – Facts and Formulae, 2nd ed. (Birkhäser, Basel, 2002). 65 Handbook of Mathematical Functions, edited by M. Abramowitz and I. A. Stegun (Dover, New York, 1970).
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
Additional information on J. Chem. Phys. Journal Homepage: http://jcp.aip.org/ Journal Information: http://jcp.aip.org/about/about_the_journal Top downloads: http://jcp.aip.org/features/most_downloaded Information for Authors: http://jcp.aip.org/authors
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
THE JOURNAL OF CHEMICAL PHYSICS 137, 184103 (2012)
Green’s function for reversible geminate reaction with volume reactivity Svetlana S. Khokhlova and Noam Agmon The Fritz Haber Research Center, Institute of Chemistry, The Hebrew University of Jerusalem, Jerusalem 91904, Israel
(Received 27 September 2012; accepted 15 October 2012; published online 12 November 2012) The kinetics of a diffusing particle near a reversible trap may be described by an extension of the Feynman-Kac equation to the case of reversible binding, which can occur within a finite reaction sphere. We obtain the Green’s function solution for the Laplace transform of this equation when the particle is initially either bound or unbound. We study the solution in the time-domain by either inverting the Laplace transform numerically or propagating the partial differential equation in the time-domain. We show that integrals of this solution over the reaction sphere agree with previously obtained solutions. © 2012 American Institute of Physics. [http://dx.doi.org/10.1063/1.4764357] I. INTRODUCTION
Chemical reactions in solution are controlled by translational diffusion, which brings the reacting partners into a critical distance from which reaction may occur.1 A fundamental model for reaction involves “contact reactivity,” in which a pair of particles has first to collide in order to react. The reaction is then imposed as a boundary condition at the contact distance σ .2 This model can be generalized to a reversible diffusion-influenced reaction3 kr
A + B C,
(1.1)
kd
in which kr and kd are recombination and dissociation rate coefficients, respectively. In the case of a geminate pair, when we deal with a single C molecule or an isolated single pair of A and B, the solution can be obtained in Laplace space and inverted into the time-domain exactly, approximately, or numerically under a variety of conditions.4–23 This model has been applied extensively to (ultra-fast) excited-state proton transfer reactions in solution,24–40 within cavities,41–43 in ice,44, 45 and within the green fluorescent protein.46, 47 It has also been utilized for analyzing translational diffusion correlation functions in molecular simulations.48, 49 Despite this success, it is worthwhile to consider an alternate model for reversible geminate recombination. As opposed to the “contact reactivity” model, where reaction occurs at a specified A-B distance, r = σ , in the “volume reactivity” model reaction can occur throughout the spherical volume.50 In this model a particle diffuses freely in space, but when inside a “reaction sphere,” r ≤ σ , it can react reversibly with the rate coefficients kr and kd . The volume reactivity model can have different implementations with respect to the distance at which the pair A–B is released upon dissociation. Reference 51 considered the case when the release is at a specified distance, either smaller or larger than σ , and solved the ensuing diffusionreaction equations for the steady-state concentrations. Reference 50 considered the cases that the released particle is uniformly distributed within the reaction sphere or else when it appears at the same distance where it was trapped. This 0021-9606/2012/137(18)/184103/8/$30.00
immobilization/remobilization problem is the one treated here. The problem has been previously studied in Laplace space for a particle initially uniformly bound within the reaction sphere.50 Here we extend the solution to the Green’s function52 (GF) p(r,t|r0 ), namely, the probability that the particle transits from an initial distance r0 to r in threedimensional space during time t, in the presence of a reactive sphere of radius σ . Interestingly enough, the integral of the GF over r for kd = 0 obeys the celebrated Feynman-Kac equation (FKE).53–55 Its solution gives the “residence time” of the particle within the reaction sphere,56 called “occupation time” in random walk theory.57 The present work generalizes the solution, which can be found in the appendix of Ref. 55, to reversible binding and an initial δ-function distribution. This solution, which can only be obtained in Laplace space, is inverted numerically to the time-domain in order to exemplify its behavior for the various initial conditions. II. THE REVERSIBLE FEYNMAN-KAC EQUATION
In the “volume reactivity model,”50 a particle diffusing (with diffusion coefficient D) in three-dimensional space may bind and unbind only within a sphere of radius r = σ (centered at the origin), whose volume is V = 4π σ 3 /3. Thus we limit the discussion to spherical symmetry, where r is the norm of the radius vector. Let p(r, t) denote the probability density to find the particle unbound and at distance r from the origin at time t, while its probability density for being bound at some r < σ is denoted by q(r, t). Unbound and bound populations interchange at rates kr H(σ − r)p(r, t) and kd q(r, t) in the association and dissociation directions, respectively. H(x) is the Heaviside step-function, which is zero for x < 0 and 1 for x ≥ 0. It restricts the association reaction to occur only within the σ -sphere. No such restriction is needed in the dissociation direction because an initially unbound particle cannot bind at r > σ . The dissociating particle always resumes its diffusive motion at the same point to which it was bound. Considered as a model for reversible geminate recombination, here one particle (e.g., a ligand) can diffuse inside the
137, 184103-1
© 2012 American Institute of Physics
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-2
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
volume of the other (e.g., its binding site). For reaction between two small molecules this may seem artificial, yet of little quantitative consequence, because the volume is small. In other cases the present model may be more realistic than the contact reactivity model, for example, reversible binding of a diffusing particle to a circular receptor array.58 Using these notations, the following two coupled differential equations, a partial differential equation (PDE) plus an ordinary one, describe the reversible reaction in the volume reactivity model: ∂p(r, t) = L p(r, t) − kr H (σ − r) p(r, t) + kd q(r, t), ∂t (2.1a) ∂q(r, t) = kr H (σ − r) p(r, t) − kd q(r, t). ∂t
∂ 1 ∂ D(r)r 2 , 2 r ∂r ∂r
(2.2)
D(r) is the particle’s diffusion coefficient (for a pair it is their relative diffusion coefficient), which is hereafter assumed to be independent of r. In the irreversible limit there is no dissociation, so that kd = 0. Then the equations are no longer coupled: it suffices to solve Eq. (2.1a) for p(r, t) whereas q(r, t) could then be obtained by integrating the outgoing flux 4π r2 D∂p(r, t)/∂r with respect to t (for any r ≤ σ ). In this limit, Eq. (2.1a) is known as the Feynman-Kac equation. The boundary conditions (BCs) imposed on it are zero flux at the origin 2 ∂p(r, t) →0 (2.3) r ∂r r→0 and vanishing probability at infinite distances p(r, t) → 0 as
r → ∞.
(2.4)
There is no need to impose BCs on q(r, t) because in the bound state the particle does not diffuse. We solve these equations in Laplace space. The Laplace transformation (LT) of a function f(t), denoted by fˆ(s), is ∞ ˆ f (s) = f (t)e−st dt. 0
Taking the LT of Eqs. (2.1) we obtain ˆ s) = L p(r, ˆ s) − kr H (σ − r) p(r, ˆ s) s p(r, ˆ s) + p(r, 0), + kd q(r, ˆ s) = q(r,
For a particle that is initially unbound and located at distance r0 from the origin, we assume δ(r − r0 ) , 4π r02
p(r, 0|r0 ) =
q(r, 0|r0 ) = 0.
(2.5a)
1 ˆ s) + q(r, 0)]. [kr H (σ − r) p(r, s + kd (2.5b)
Here we still need to insert the initial conditions p(r, 0) and q(r, 0). For the GF, these involve a particle that is either unbound or bound at distance r0 .
(2.6)
By convention, the initial separation is denoted on the right of the vertical line, so that p(r, t|r0 ) depicts the GF in the timeˆ s|r0 ) is the GF in Laplace space. Insertdomain, whereas p(r, ing Eqs. (2.5b) and (2.6) into Eq. (2.5a) gives ˆ s|r0 ) − (s − L)p(r,
δ(r − r0 ) skr H (σ − r) ˆ s|r0 ), p(r, =− 2 s + kd 4π r0 (2.7)
(2.1b)
Here L is the spherically-symmetric diffusion operator in a three-dimensional space L≡
A. Initially unbound particle
which is an ordinary (because s is just a parameter) inhomogeneous differential equation in the variable r. This reversible FKE extends the conventional FKE in two respects. First, it is reversible, though this only amounts to replacing kr by skr /(s + kd ). Second, it applies to a δ-function initial condition, from which the solution for any other initial distribution can be obtained by averaging the GF over this distribution. The spatial integral of p(r, t|r0 ) over r is the “survival probability” of the unbound particle ∞ P (t|r0 ) ≡ 4π p(r, t|r0 )r 2 dr. (2.8) 0
The customary FKE deals with this function,54, 55 which is discussed in Appendix A, whereas Appendix B deals with a uniform initial distribution within the reaction sphere.50 It is convenient to solve Eq. (2.7) for r = r0 when it becomes homogeneous and then match the solutions at r0 . Denote the solutions for r < r0 and r > r0 by pˆ − (r, s|r0 ) and pˆ + (r, s|r0 ), respectively. They must be continuous at r0 , but their derivative is discontinuous there, pˆ − (r, s|r0 ) = pˆ + (r, s|r0 ),
(2.9a)
1 ∂ pˆ − (r, s|r0 ) ∂ pˆ + (r, s|r0 ) = − . (2.9b) ∂r ∂r 4π Dr02 r=r0 − r=r0 + Because Eq. (2.7) is of second order, each of these solutions is a linear combination of two fundamental solutions (see below) with coefficients determinable from the above condition. Due to the discontinuous Heaviside function, there is another matching sphere at r = σ , but there one demands that both the function and its derivative are continuous. Thus we solve the equation separately in three regions of space bounded by the spheres r = r0 and r = σ and match the solutions at these two interfaces. Finally, after solving for ˆ s|r0 ), one can obtain the probability density of the bound p(r, particle from Eq. (2.5b), namely ˆ s|r0 ) = q(r,
kr ˆ s|r0 ) , r ≤ σ. p(r, s + kd
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
(2.10)
184103-3
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
B. Initially bound particle
A. Unbound particle, initially outside: r0 > σ
For a particle that is initially bound and at r0 ≤ σ we impose
For r0 > σ , we consider the solution in the three regions: (1) r < σ , (2) σ < r < r0 , and (3) r > r0 , which are depicted in Fig. 1. In region 3 the BC at infinity eliminates the rising exponential, so that
p(r, 0|r0 ) = 0 ,
q(r, 0|r0 ) =
δ(r − r0 ) . 4π r02
(2.11)
ˆ s|r0 ) = p(r,
Instead of Eq. (2.7), we now obtain ˆ s|r0 ) − (s − L)p(r,
kd δ(r − r0 ) s + kd 4π r02
which has a modified δ-function term. Hence its solution is simply that of Eq. (2.7) times kd /(s + kd ). The probability density of the bound particle is subsequently obtained from Eq. (2.5b), namely, 1 δ(r − r0 ) ˆ s|r0 ) + ˆ s|r0 ) = kr p(r, . (2.13) q(r, s + kd 4π r02
ˆ s|r0 ) = p(r,
ˆ s|r0 ) = p(r,
Given the differential operator in Eq. (2.2), the two fundamental solutions of this equation are exp (± αr)/r, where (3.2) α ≡ s/D.
A sinh(βr), r < σ. r
(3.5c)
1 e−α(r0 −σ ) [α sinh(βσ ) + β cosh(βσ )]−1 , 4π Dr0 (3.6a)
B1 =
1 e−αr0 , 8π Dαr0
B2 =
α sinh(βσ ) − β cosh(βσ ) −α(r0 −2σ ) 1 , e 8π Dαr0 α sinh(βσ ) + β cosh(βσ )
(3.6b)
(3.6c) (3.3)
and the two fundamental solutions are exp (± βr)/r, where β is defined by β = α[1 + kr /(s + kd )]1/2 .
(3.5b)
The coefficients A, B1 , B2 , and C are functions of s and r0 (but not of r), to be determined from the matching conditions at r = r0 and r = σ . Thus we find that
Let us consider Eqs. (2.7) and (2.12) for r = r0 and either r > σ or r < σ . For r > σ they simplify to (3.1)
B1 αr B2 −αr e + e , σ < r < r0 . r r
Finally, in region 1 the vanishing flux at the origin dictates that
A=
III. THE GREEN’S FUNCTION SOLUTION
For r < σ they become kr ˆ s|r0 ) = 0 p(r, L−s 1+ s + kd
(3.5a)
In region 2 both exponentials contribute,
skr H (σ − r) ˆ s|r0 ), (2.12) p(r, =− s + kd
ˆ s|r0 ) = 0. (L − s)p(r,
C exp(−αr), r > r0 . r
(3.4)
When kr = 0, β = α and Eq. (2.7) reduces to free diffusion. Let us now solve these equations separately for r0 > σ or r0 < σ (Figure 1).
C = B2 +
1 eαr0 . 8π Dαr0
(3.6d)
These coefficients obey the following four identities: (i) From continuity of the solution and its derivative jump at r = r0 : Ce−αr0 = B1 eαr0 + B2 e−αr0 , Ce−αr0 =
1 − B1 eαr0 + B2 e−αr0 . 4π αDr0
(3.7a) (3.7b)
(ii) From continuity of the solution and its first derivative at r = σ:
FIG. 1. The spatial arrangement of elements pertaining to the GF solution, which depends on the two distances r and r0 . The initial location of the particle, r = r0 , can be either inside or outside the reactive sphere, r = σ . This divides space, in two different ways, into three regions (indicated by 1, 2, and 3). The particle can react only within the colored sphere, where r ≤ σ .
A sinh(βσ ) = B1 eασ + B2 e−ασ ,
(3.8a)
Aβ cosh(βσ ) = B1 αeασ − B2 αe−ασ .
(3.8b)
B. Unbound particle, initially inside: r0 < σ
For r0 < σ we consider the solution in the three regions: (1) r < r0 , (2) r0 < r < σ , (3) r > σ (see Fig. 1). The solutions are the same as in Eqs. (3.5), except for region 2, where α is
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-4
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
replaced by β, ˆ s|r0 ) = p(r,
C exp(−αr), r
r > σ,
(3.9a)
ˆ s|r0 ) = p(r,
B1 −β(σ −r) B2 β(σ −r) e e + , r r r0 < r < σ,
(3.9b)
ˆ s|r0 ) = p(r,
A sinh(βr), r
r < r0 .
(3.9c)
The matching conditions now give A=
1 α sinh[β(σ − r0 )] + β cosh[β(σ − r0 )] , 4π Dβr0 α sinh(βσ ) + β cosh(βσ ) (3.10a)
B1 =
β − α −ασ e C, 2β
(3.10b)
B2 =
β + α −ασ e C, 2β
(3.10c)
eασ sinh(βr0 ) . 4π Dr0 α sinh(βσ ) + β cosh(βσ )
(3.10d)
C=
These coefficients obey the following four identities: (i) From continuity of the solution and its first derivative at r = σ: Ce−ασ = B1 + B2 ,
(3.11a)
αCe−ασ = β(B2 − B1 ).
(3.11b)
(ii) From continuity of the solution and its derivative jump at r = r0 : A sinh(βr0 ) = B1 e−β(σ −r0 ) + B2 eβ(σ −r0 ) , A cosh(βr0 ) =
(3.12a)
1 + B1 e−β(σ −r0 ) − B2 eβ(σ −r0 ) . 4πβDr0
FIG. 2. The distance dependence of the Laplace transformed GF for an initially unbound particle at three values of r0 : inside (r0 = 0.5), on the surface (r0 = 1), or outside (r0 = 2) the reaction sphere. Lines demonstrate Eqs. (3.5) and (3.6), or (3.9) and (3.10), for s = 1 (full line), 5 (dashed line), and 10 (dashed–dot line). All the figures in this publication use the same parameters set: σ = 1, D = 1, kd = 0.5, and kr = 3.
Figure 2 shows the r-dependence of the solution in ˆ Laplace space, p(r,s|r 0 ), for a particle that is initially unbound at r0 . The solution for r0 inside or outside the reaction sphere is given by Eqs. (3.5) and (3.6) or (3.9) and (3.10), ˆ respectively. For large s, the function p(r,s|r 0 ) peaks at r = r0 and spreads over r as s decreases. The derivative discontinuity at r = r0 is evident in the figure, whereas there are no discontinuities at r = σ . Figures 3–5 show the GF solutions in the time-domain, for the same three initial locations. The black lines are ˆ ˆ the (numerical) Laplace inverse of p(r,s|r 0 ) and q(r,s|r 0 ). The colored dotted lines are the exact numerical solution of the PDE in the time-domain using SSDP (version 2.66).60 There is an exact match between the two solutions. The GF for an unbound particle, p(r,t|r0 ) (red), looks initially like a widening Gaussian. Its peak moves to smaller
(3.12b) inverse LT SSDP p(r,t|r0)
t1 = 0.04663
The solution of Eq. (2.12) is the same as in Eqs. (3.10) above, with A and C multiplied by kd /(s + kd ). IV. DEMONSTRATIVE CALCULATIONS
In this section we demonstrate the above GF solution in both Laplace space and the time-domain. The latter is achieved by a numerical LT inversion using the algorithm of Ref. 59. It is compared with the numerical solution of Eq. (2.1) in the time-domain using the SSDP software.60 The agreement validates our analytical LTs (or, alternately, it validates the numerical procedures for LT inversion and PDE solution). The parameters used in the following figures are the same as in our previous publication:50 σ = 1, D = 1, kd = 0.5, and kr = 3D/σ 2 = 3.
t2 = 0.2489
p (r, t | r0 ), q (r, t | r0 )
C. Bound particle, initially inside: r0 < σ
0.02
1
SSDP q(r,t|r0)
t3 = 0.5749
0.01
2
3
3 2 0.00
0
1
2
3
4
r FIG. 3. The distance dependence of the GF in the time-domain, for an initially unbound particle starting outside the reaction sphere, at r0 = 2 > σ . Black lines demonstrate the Laplace inverse of Eqs. (3.5) and (3.6) for p(r, t|r0 ), as well as that of Eq. (2.10) for the probability density of the bound state, q(r, t|r0 ).
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-5
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012) 0.10
t4 = 4.659 t5 = 6.881 t6 = 11.7 1
q ( r, t | r0 )
p ( r, t | r0 )
t3 = 1.469
0.02
0.01
0.08
t2 = 0.5185
3
4
4
t1 = 0.02102
0.03 2
3 5
0.06
0.04
2 6
5
0.02
6 1 0.00 0.00 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.2 0.4 0.6 0.8 1.0 1.2
r FIG. 4. Same as Fig. 3 for r0 = 1 = σ .
values of r with increasing time due to reflection from r = 0. In contrast to these smooth profiles, the GF for the bound particle, q(r,t|r0 ) (green), has a cusp. This cusp is on the boundary of the reaction sphere (r = σ ) when r0 ≥ σ and at r0 when r0 < σ . Due to the recombination reaction, part of the initial δ-function density for the unbound particle becomes a δfunction for the bound particle. This population cannot diffuse in the bound state, only dissociate, giving rise to a derivative discontinuity in the time-domain. The time evolution of the GF for an initially bound particle (at r0 ≤ σ ) is shown in Fig. 6. The initial distribution is a δ-function q(r, 0|r0 ) = δ(r − r0 )/4π r02 (at r0 = 0.5) that persists for all times. Its amplitude, however, diminishes with time (not shown) as the particle dissociates. This creates the distribution p(r, t|r0 ) of the unbound particle, which is initially close to a δ-function (time t1 ), then widens (due to diffusion) and increases in amplitude (due to further dissociation – time t2 ), and finally (t3 onwards) decreases in amplitude due to escape of the particle to infinity. Recombination from the unbound state then creates the non-δ-function part of q(r, t|r0 ), its amplitude following that of p(r, t|r0 ) with a delay (due to the finite kr ). Thus the amplitude of q(r, t|r0 ) continues to in0.7
1 inverse LT SSDP p(r,t|r0)
0.6
2
p (r, t | r0 ), q ( r, t| r0 )
0.5
0.4
SSDP q(r,t|r0)
3
0.3
0.2
0.1
0.0
0.0
3 1 0.5
1.0
1.5
r FIG. 5. Same as Fig. 3 for r0 = 0.5 < σ , except that the times here are: t1 = 0.02019, t2 = 0.0307, and t3 = 0.05361. Black lines were calculated by Laplace inverting Eqs. (3.9) and (3.10).
r
FIG. 6. The time dependence of the GF for a particle that is initially bound ˆ s|r0 ) is given by at r0 = 0.5. The LT of the GF for the unbound state p(r, Eq. (3.9) multiplied by kd /(s + kd ) and its Laplace inverse is depicted by the ˆ s|r0 ), is calculated from blue lines. The LT of the GF for the bound state, q(r, ˆ s|r0 ) using Eq. (2.13) and its Laplace inverse is depicted by the red lines. p(r, The numerical solution of the PDE (2.1) for these initial conditions is shown by the dotted green lines.
crease up to about t4 , when p(r, t|r0 ) is already in its decreasing phase. Another interesting observation is that both parts of the GF show a derivative discontinuity at r0 . This contrasts with the initially unbound particle, for which only q(r, t|r0 ) shows such a discontinuity (Figs. 3–5). V. CONCLUSION
In this work we have obtained (analytically) the GF solution of the reversible FKE in Laplace space and (numerically) in the time-domain. The relevance of this equation to chemical kinetics is in providing a mathematical description for the “volume reactivity” model of geminate recombination, in which one particle is a reactive sphere and the other is a diffusing point particle that can bind (immobilize) and unbind (remobilize) within this sphere. This model is an alternative to the widely used “contact reactivity” model, in which the two particles can react only upon collision (at the “contact distance”). Interestingly, the GF solution for reversible binding of an initially unbound particle can be obtained from that of irreversible binding (kd = 0) by replacing kr /s by kr /(s + kd ). An initially bound particle is immobilized for all times in the irreversible case, yet its reversible GF can also be obtained from the solution of the irreversible case by a further multiplicative factor kd /(s + kd ). In Laplace space, all the components of the GF necessarily show a derivative discontinuity at the initial location of the particle, r0 . In the time-domain, p(r, t|r0 ) does not show such a discontinuity for an initially unbound particle, but all other GF components do show a derivative discontinuity at r0 , which is conceivably due to the lack of diffusion in the bound state. Appendices A and B show that the integrals of our GF’s (over r or r0 ) reproduce the solutions for the overall survival probability of an initially unbound particle (obtained in Ref. 55 for the irreversible case) and for an initially uniformly
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-6
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
distributed bound particle.50 Thus all previous solutions of the FKE are special cases of the GF obtained in the present work. The question is what is the practical difference between the contact and volume reactivity models and is a volume reactivity model really needed? As Fig. 5 in reference 50 shows, for an initially uniformly distributed particle within the reaction sphere the difference in the binding probability, Q(t|uni), between the two models may be rather subtle (see Fig. 5). However, the difference between the models may be larger for an initially unbound particle or for the full GF. This issue still needs to be studied. More interesting physically is the applicability of the volume reactivity model beyond the pair problem. Consider, for example, the case where reversible receptors are scattered within the reaction sphere.58 For example, a particle may be diffusing on the surface of a membrane,61 on which there is a circular patch containing many reversible receptors. (Such a situation is believed to occur, e.g., for protons on membranes.62, 63 ) Then, in a mean field sense, the kinetics of the particle is described by a two-dimensional version of the reversible FKE solved herein.
D −1/2 = βK1/2 (ασ )I3/2 (βσ ) + αK3/2 (ασ )I1/2 (βσ ). (A3) Bessel functions of the half-integer order can be written in terms of hyperbolic functions [Eqs. (10.2.13) and (10.2.17) in Ref. 65]. Then in Eq. (A3) reduces to Dβ 1/2 [β cosh(βσ ) + α sinh(βσ )]e−ασ . βσ = α (A4) The LTed binding probability in Eq. (A2) then becomes 1 sinh(βr0 ) kr 1 + ασ ˆ , ) = − Q(s|r 0 s s + kr β 2 r0 D β cosh(βσ ) + α sinh(βσ ) r0 ≤ σ, (A5a) kr e−α(r0 −σ ) βσ cosh(βσ ) − sinh(βσ ) ˆ Q(s|r , r0 ≥ σ. 0) = s β 2 r0 D β cosh(βσ ) + α sinh(βσ )
ACKNOWLEDGMENTS
Research supported by The Israel Science Foundation (Grant No. 122/08). The Fritz Haber Research Center is supported by the Minerva Gesellschaft für die Forschung, München, FRG. APPENDIX A: THE INTEGRATED FKE
The FKE usually encountered in the literature is for the LT of the survival probability, Pˆ (s|r0 ), obtained by integrating the GF as in Eq. (2.8). It is a solution of the integrated adjoint of Eq. (2.7), namely, skr H (σ − r0 ) s+ − L0 Pˆ (s|r0 ) = 1, (A1) s + kd with the adjoint diffusion operator L0 and the same boundary conditions at r0 = 0 and ∞ as above. This solution, for kd = 0, is given in the appendix of Ref. 55 for any dimensionality d. See also Eq. (1.5.1) in Chap. 4 of Part II of Ref. ˆ 64. The binding probability Q(s|r 0 ) can be obtained from the ˆ normalization condition Q(s|r0 ) = 1/s − Pˆ (s|r0 ). It is actually slightly simpler than the survival probability, so we focus on it here. For the irreversible case (kd = 0), Eqs. (59) and (60) in Ref. 55 in three dimensions become
αK3/2 (ασ ) σ 1/2 1 kr ˆ − 2 1/2 I1/2 (βr0 ) , Q(s|r0 ) = s s + kr β D r0 r0 ≤ σ, kr I3/2 (βσ ) ˆ Q(s|r 0) = s β D 1/2
Here Iν (z) and Kν (z) are the modified Bessel functions of order ν, and in the denominator is given by Eq. (61) of Ref. 55:
σ r0
(A2a)
1/2 K1/2 (αr0 ) , r0 ≥ σ. (A2b)
(A5b) The solution for reversible binding, kd = 0, is obtained by replacing kr /s by kr /(s + kd ) in the above equations, whereas β generalizes to its value in Eq. (3.4). This gives 1 + ασ kr ˆ 1− Q(s|r ) = 0 s(s + kd + kr ) r0 sinh(βr0 ) , r0 ≤ σ, × (A6a) β cosh(βσ ) + α sinh(βσ ) ˆ Q(s|r 0) =
e−α(r0 −σ ) βσ cosh(βσ ) − sinh(βσ ) kr , s(s +kd +kr ) r0 β cosh(βσ ) + α sinh(βσ ) r0 ≥ σ. (A6b)
Equation (A6) can also be obtained by direct integration of the GF over r. According to Eq. (2.10), obtaining the binding probability requires integration only inside the reaction sphere σ 4π kr ˆ ˆ s|r0 )r 2 dr, p(r, Q(s|r ) = (A7) 0 s + kd 0 and this is simpler than obtaining the survival probability that requires integration over all space. When r0 < σ , one needs ˆ s|r0 ) in regions 2 and 3, see Eqs. (3.9b) and to integrate p(r, (3.9c), but only in region 1, Eq. (3.5c), when r0 > σ . The surˆ vival probability can subsequently be obtained from Q(s|r 0) ˆ using the normalization condition Pˆ (s|r0 ) = 1/s − Q(s|r 0 ). Performing these integrals directly indeed gives the same result as above. Its Laplace inverse, P(t|r0 ), is demonstrated in Fig. 7. It starts from unity, then decreases to a minimum as the particle binds. But eventually it is released and escapes to infinity when P(t|r0 ) approaches unity again.
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-7
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012) 5 N.
1.0
P (t | r0 )
0.9
r0 = 1
0.8
0.7
SSDP P(t | r0)
^ |r ) inverse LT of P(s 0
r0 = 0.5
0.6
0
5
10
15
t FIG. 7. The time dependence of the survival probability, P(t|r0 ), obtained by Laplace inverting Eqs. (A5), as compared to numerically integrating the numerical solution of Eq. (2.1) using the SSDP software. Parameters are the same as in Fig. 2.
APPENDIX B: AN INITIAL UNIFORMLY BOUND STATE
Instead of a particle initially bound at r0 , we consider a bound particle that is initially uniformly distributed inside the reaction sphere. Denoting the survival probability in this case by P(r, t|uni), we define σ 4π P (r, t|uni) ≡ p(r, t|r0 ) r02 dr0 , (B1) V 0 where the division by V = 4π σ 3 /3 serves to normalize the initial distribution. This problem was considered in Ref. 50. Because the diffusion equation is symmetric to the exchange of r and r0 , we find from Eq. (A7) that s + kd ˆ Pˆ (r, s|uni) = Q(s|r). V kr
(B2)
Inserting the binding probability from Eq. (A6) with r0 replaced by r and multiplying by kd /(s + kd ) to account for an initially bound state, we obtain 1 + ασ kd /V 1− Pˆ (r, s|uni) = s(s + kd + kr ) r sinh(βr) × , r ≤ σ, β cosh(βσ ) + α sinh(βσ ) (B3a) Pˆ (r, s|uni) =
kd /V e−α(r−σ ) s(s + kd + kr ) r ×
βσ cosh(βσ ) − sinh(βσ ) , r ≥ σ. (B3b) β cosh(βσ ) + α sinh(βσ )
This is identical to Eq. (4.14) in Ref. 50. 1 S.
Agmon and A. Szabo, J. Chem. Phys. 92, 5270 (1990). Nadler and D. L. Stein, J. Chem. Phys. 104, 1918 (1996). 7 B. C. Lagerholm and N. L. Thompson, Biophys. J. 74, 1215 (1998). 8 I. V. Gopich, K. M. Solntsev, and N. Agmon, J. Chem. Phys. 110, 2164 (1999). 9 N. Agmon, J. Chem. Phys. 110, 2175 (1999). 10 I. V. Gopich and N. Agmon, J. Chem. Phys. 110, 10433 (1999). 11 H. Kim, K. J. Shin, and N. Agmon, J. Chem. Phys. 111, 3791 (1999). 12 H. Kim and K. J. Shin, Phys. Rev. Lett. 82, 1578 (1999). 13 N. Agmon and I. V. Gopich, J. Chem. Phys. 112, 2863 (2000). 14 H. Kim and K. J. Shin, J. Chem. Phys. 112, 8312 (2000). 15 H. Kim, K. J. Shin, and N. Agmon, J. Chem. Phys. 114, 3905 (2001). 16 H. Kim and K. J. Shin, J. Chem. Phys. 120, 9142 (2004). 17 N. Agmon, Chem. Phys. Lett. 417, 530 (2006). 18 S. Y. Reigh, K. J. Shin, and M. Tachiya, J. Chem. Phys. 129, 234501 (2008). 19 S. Park and N. Agmon, J. Chem. Phys. 130, 074507 (2009). 20 K. Park, K. J. Shin, and H. Kim, J. Chem. Phys. 131, 154105 (2009). 21 H. Kim, Chem. Phys. Lett. 507, 265 (2011). 22 I. Plante and F. A. Cucinotta, Phys. Rev. E 84, 051920 (2011). 23 Z. Kalay, J. Phys. A 45, 235001 (2012). 24 E. Pines, D. Huppert, and N. Agmon, J. Chem. Phys. 88, 5620 (1988). 25 D. Huppert, E. Pines, and N. Agmon, J. Opt. Soc. Am. B 7, 1545 (1990). 26 A. Masad and D. Huppert, Chem. Phys. Lett. 180, 409 (1991). 27 N. Agmon, D. Huppert, A. Masad, and E. Pines, J. Phys. Chem. 95, 10407 (1991); Erratum 96, 2020 (1992). 28 S. Y. Goldberg, E. Pines, and D. Huppert, Chem. Phys. Lett. 192, 77 (1992). 29 N. Agmon, S. Y. Goldberg, and D. Huppert, J. Mol. Liq. 64, 161 (1995). 30 I. Carmeli, D. Huppert, L. M. Tolbert, and J. E. Haubrich, Chem. Phys. Lett. 260, 109 (1996). 31 D. Huppert, L. M. Tolbert, and S. Linares-Samaniego, J. Phys. Chem. A 101, 4602 (1997). 32 K. M. Solntsev, D. Huppert, N. Agmon, and L. M. Tolbert, J. Phys. Chem. A 104, 4658 (2000). 33 K. M. Solntsev, D. Huppert, and N. Agmon, Phys. Rev. Lett. 86, 3427 (2001). 34 N. Koifman, B. Cohen, and D. Huppert, J. Phys. Chem. A 106, 4336 (2002). 35 C. Clower, K. M. Solntsev, J. Kowalik, L. M. Tolbert, and D. Huppert, J. Phys. Chem. A 106, 3114 (2002). 36 K. M. Solntsev, E. N. Sullivan, L. M. Tolbert, S. Ashkenazi, P. Leiderman, and D. Huppert, J. Am. Chem. Soc. 126, 12701 (2004). 37 L. Genosar, P. Leiderman, N. Koifman, and D. Huppert, J. Phys. Chem. A 108, 1779 (2004). 38 N. Agmon, J. Phys. Chem. A 109, 13 (2005). 39 P. Leiderman, R. Gepshtein, A. Uritski, L. Genosar, and D. Huppert, J. Phys. Chem. A 110, 9039 (2006). 40 M. R. di Nunzio, B. Cohen, and A. Douhal, J. Phys. Chem. A 115, 5094 (2011). 41 B. Cohen, D. Huppert, K. M. Solntsev, Y. Tsfadia, E. Nachliel, and M. Gutman, J. Am. Chem. Soc. 124, 7539 (2002). 42 R. Gepshtein, P. Leiderman, D. Huppert, E. Project, E. Nachliel, and M. Gutman, J. Phys. Chem. B 110, 26354 (2006). 43 N. Amdursky, R. Orbach, E. Gazit, and D. Huppert, J. Phys. Chem. C 113, 19500 (2009). 44 P. Leiderman, A. Uritski, and D. Huppert, J. Phys. Chem. A 111, 4998 (2007). 45 A. Uritski, I. Presiado, and D. Huppert, J. Phys. Chem. C 112, 11991 (2008). 46 P. Leiderman, D. Huppert, and N. Agmon, Biophys. J. 90, 1009 (2006). 47 N. Agmon, J. Phys. Chem. B 111, 7870 (2007). 48 O. Markovitch and N. Agmon, J. Chem. Phys. 129, 084505 (2008). 49 H. Chen, G. A. Voth, and N. Agmon, J. Phys. Chem. B 114, 333 (2010). 50 S. S. Khokhlova and N. Agmon, Bull. Korean Chem. Soc. 33, 1020 (2012). 51 J. Lipková, K. C. Zygalakis, S. J. Chapman, and R. Erban, SIAM J. Appl. Math. 71, 714 (2011). 52 M. C. M. Wright, Nat. Phys. 2, 646 (2006). 53 M. Kac, in Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, edited by J. Neyman (University of California, Berkeley, 1951), pp. 189–215. 54 S. N. Majumdar, Curr. Sci. 89, 2076 (2005), arXiv:cond-mat/0510064. 55 N. Agmon, J. Phys. Chem. A 115, 5838 (2011). 56 N. Agmon, J. Chem. Phys. 81, 3644 (1984). 57 S. Karlin and H. M. Taylor, A Second Course in Stochastic Processes (Academic, San-Diego, 1981). 58 A. L. Edelstein and N. Agmon, J. Comput. Phys. 132, 260 (1997). 6 W.
r0 = 2
A. Rice, in Diffusion-Limited Reactions, Comprehensive Chemical Kinetics Vol. 25 (Elsevier, Amsterdam, 1985). 2 F. C. Collins and G. E. Kimball, J. Colloid Sci. 4, 425 (1949). 3 N. Agmon, J. Chem. Phys. 81, 2811 (1984). 4 N. Agmon and G. H. Weiss, J. Chem. Phys. 91, 6937 (1989).
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
184103-8
S. S. Khokhlova and N. Agmon
J. Chem. Phys. 137, 184103 (2012)
59 J.
63 N.
60 E.
64 A.
Abate and W. Whitt, ORSA J. Comput. 7, 36 (1995). B. Krissinel’ and N. Agmon, J. Comput. Chem. 17, 1085 (1996). 61 N. Agmon, Chem. Phys. 370, 232 (2010). 62 A. Springer, V. Hagen, D. A. Cherepanov, Y. N. Antonenko, and P. Pohl, Proc. Natl. Acad. Sci. U.S.A. 108, 14461 (2011).
Agmon and M. Gutman, Nat. Chem. 3, 840 (2011). N. Borodin and P. Salminen, in Handbook of Brownian Motion – Facts and Formulae, 2nd ed. (Birkhäser, Basel, 2002). 65 Handbook of Mathematical Functions, edited by M. Abramowitz and I. A. Stegun (Dover, New York, 1970).
Downloaded 25 Nov 2012 to 132.64.96.165. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
