Theory of diffusion-influenced reactions in complex geometries
aa r X i v : . [ phy s i c s . c h e m - ph ] O c t Theory of diffusion-influenced reactions in complex geometries
Marta Galanti , , , , Duccio Fanelli , , Francesco Piazza Universit`a degli Studi di Firenze, Dipartimento di Fisica e Astronomia and CSDC,via G. Sansone 1, IT-50019 Sesto Fiorentino, Firenze, Italia, Dipartimento di Sistemi e Informatica, Universit`a di Firenze,Via S. Marta 3, IT-50139 Florence, Italy, INFN, Sezione di Firenze,Italia, Universit´e d’Orl´eans, Centre de Biophysique Mol´eculaire,CNRS-UPR4301, Rue C. Sadron, 45071, Orl´eans, France.
Chemical reactions involving diffusion of reactants and subsequent chemical fixation steps are gen-erally termed “diffusion-influenced” (DI). Virtually all biochemical processes in living media canbe counted among them, together with those occurring in an ever-growing number of emergingnano-technologies. The role of the environment’s geometry (obstacles, compartmentalization) anddistributed reactivity (competitive reactants, traps) is key in modulating the rate constants of DIreactions, and is therefore a prime design parameter. Yet, it is a formidable challenge to build acomprehensive theory able to describe the environment’s “reactive geometry”. Here we show thatsuch a theory can be built by unfolding this many-body problem through addition theorems forspecial functions. Our method is powerful and general and allows one to study a given DI reactionoccurring in arbitrary “reactive landscapes”, made of multiple spherical boundaries of given sizeand reactivity. Importantly, ready-to-use analytical formulas can be derived easily in most cases.
PACS numbers: 82.40.Qt, 82.39.Rt, 87.10.Ed, 02.30.Jr
Diffusion-influenced reactions (DIR) are ubiquitous inmany contexts in physics, chemistry and biology [1, 2]and keep on sparking intense theoretical and computa-tional activity in many fields [3–11]. Modern examples ofemerging nanotechnologies that rely on controlled alter-ations of diffusion and reaction pathways in DIRs includedifferent sorts of chemical and biochemical catalysis in-volving complex nano-reactors [12, 13], nanopore-basedsequencing engines [14] and morphology control and sur-face functionalization of inorganic-based delivery vehiclesfor controlled intracellular drug release [15, 16].However, while the mathematical foundations for thedescription of such problems have been laid nearly a cen-tury ago [17], many present-day problems of the utmostimportance at both the fundamental and applied levelare still challenging. Notably, arduous difficulties arisein the quantification of the important role played bythe environment’s geometry (obstacles, compartmental-ization) [18] and distributed reactivity (patterns of com-petitive reaction targets or traps) in coupling transportand reaction pathways in many natural and artificial(bio)chemical reactions [1, 19, 20].A formidable challenge in modeling environment-related effects on chemical reactions is represented bythe intrinsic many-body nature of the problem. This isbrought about essentially by two basic features, commonto virtually all realistic situations, namely ( i ) finite den-sity of reactants and other inert species (in biology alsoreferred to as macromolecular crowding [21, 22]) and ( ii )confining geometry of natural or artificial reaction do-mains in 3 D space. In general, the presence of multiplereactive and non-reactive particles/boundaries cannot beneglected in the study of (bio)chemical reactions occur-ring in real milieux , where the geometrical compactness of the environment may have profound effects, such as first-passage times that are non-trivially influenced bythe starting point [23]. Relevant complex media includethe cell cytoplasm [8, 24–26], porous or other artificialconfining media [23, 27–32], which be considered as of-fering important tunable features for technological appli-cations [12, 14, 16].In this paper, we take a major step forward by solvingthe general problem of computing the steady-state reac-tion rate constant for a diffusion-influenced chemical re-action between a mobile ligand and an explicit arbitrary,static 3 D configuration of spherical reactive boundariesof arbitrary sizes and intrinsic reactivities.To set the stage for the forthcoming discussion, letus first consider the simple problem of two molecules A and B of size R and a , respectively, diffusing in solution.Upon encountering, the two species can form a complex,which catalyzes the transformation of species B into someproduct P with rate constant k , A + B k −→ A + P (1)Under the hypotheses that ( i ) A molecules diffuse muchmore slowly than B molecules, ( ii ) both species arehighly diluted and ( iii ) the bulk concentration of A molecules ρ A is much smaller than the bulk concentra-tion of B molecules ρ B [2, 33], the rate constant k can becomputed by solving the following stationary two-bodyboundary problem [1] ▽ u = 0 with u | ∂ Ω = 0 , lim r →∞ u = 1 (2)where ∂ Ω is a spherical sink of radius σ = R + a (theencounter distance) and u ( r ) = ρ ( r ) /ρ B is the stationarynormalized concentration of B molecules around the sink.The rate constant is simply the total flux into the sink, i.e. k = D Z ∂ Ω ∂u∂r (cid:12)(cid:12)(cid:12)(cid:12) r = σ dS (3)where D = D A + D B is the relative diffusion constant.The solution to the boundary problem (2) is u ( r ) =1 − σ/r , which yields the so-called Smoluchowski rate con-stant for an isolated spherical sink, namely k S = 4 πDσ .These simple ideas, originally developed to describe co-agulation in colloidal systems [17, 34], together with therelated subsequent major advances by Debye [35] andCollins & Kimball [36] represent the basic building blockof many modern theoretical approaches in chemicalandsoft-matter [37, 38] kinetics.In many realistic situations in chemical and biochem-ical kinetics, a single ligand ( B ) molecule has to diffuseamong many competing reactive particles A . In addi-tion, it might be forced to find its target within a spe-cific confining geometry, which in principle can be mod-eled through a collection of reflecting boundaries. Suchsettings define a genuinely many-body problem, as theoverall flux of ligands is shaped by the mutual screeningamong all the different reactive boundaries (the reactiveenvironment ), known as diffusive interaction [39–41]. Inthe following we show how this kind of problems can beformulated and solved in a rather general form.Let us imagine a reaction of the kind (1) to be cat-alyzed at N + 1 spherical boundaries ∂ Ω α of radius (en-counter distance) σ α = R α + a , α = 0 , , . . . , N arrangedin space at positions X α . With reference to the Smolu-chowski problem, this means that we are explicitly relax-ing the assumption of vanishing density of the reactivecenters A . In the most general setting, each sphere canbe endowed with an intrinsic reaction rate constant k ∗ α ,that specifies the conversion rate from encounter complexto product at its surface. Then, the stationary densityof B molecules is the solution of the following boundaryvalue problem ▽ u = 0 (4a) (cid:18) σ α ∂u∂r α − h α u (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ∂ Ω α = 0 α = 0 , , . . . , N (4b)lim r →∞ u = 1 (4c)where h α = k ∗ α /k S α with k S α = 4 πDσ α . The boundaryconditions (BC) (4b) are called radiative or Robin bound-ary conditions. The limits h α → ∞ and h α = 0 corre-spond to perfectly absorbing (sink) and reflecting (obsta-cle) boundaries, respectively, while values 0 < h α < ∞ correspond to finite surface reactivity [36].The boundary problem (4a), (4b), (4c) provides a rig-orous mathematical description of a wide assortmentof physical situations, ranging from one or many sinksscreened by neighboring competing reactive boundariesto hindered diffusion to a sink located among a collectionof static reflecting obstacles placed at given positions inspace. In order to solve the problem, it is expedient to consider as many sets of spherical coordinate systems asthere are boundaries, r α ≡ ( r α , θ α , φ α ). The solutioncan then be written formally as an expansion in series ofirregular solid harmonics, namely u = 1 + N X α =0 u − α ( r α ) , u − α = ∞ X ℓ =0 ℓ X m = − ℓ B αmℓ r ℓ +1 α Y mℓ ( r α )(5)where Y mℓ ( r α , θ α φ α ) are spherical harmonics expressedin the local frame centered at the α -th sphere [42].The coefficients B αmℓ should be determined by impos-ing the BCs (4b). In order to do so, we use known addi-tion theorems for spherical harmonics [43] to express thesolution (5) in all the N +1 different reference frames cen-tered at each sphere. The result is the following infinite-dimensional system of linear equations B αgq − ( q − h α )( h α + q + 1) (cid:20) δ g δ q + ∞ X ℓ =0 ℓ X m = − ℓ N X β =0 β = α B βmℓ W αβgqmℓ (cid:21) = 0(6)for α = 0 , , . . . , N , q = 1 , , . . . , ∞ and g = − q, ..., q (seesupplementary material for a detailed derivation and theexplicit expression of the matrix elements W αβgqmℓ ). Tosolve the problem one simply needs to truncate the sumon ℓ in eq. (6), by including a finite number of multipolesso as to attain the desired accuracy on the overall rateconstant. In analogy to eq. (3), and taking into accountthe definition (5), the rate constant corresponding to agiven subset of reactive boundaries S can be computedas k = D X α ∈S Z ∂ Ω α ∂u∂r (cid:12)(cid:12)(cid:12)(cid:12) r = σ α dS = − X α ∈S k S α B α (7)The theoretical framework that culminates in for-mula (7) provides an extremely powerful tool to inves-tigate how specific geometries of obstacles and/or com-petitive reactive boundaries modulate the rate constantof a given diffusion-influenced reaction.A clear and instructive illustration of our general ap-proach can be outlined by focussing on the simplestmodel of diffusion-influenced reaction, namely diffusionof a ligand to a perfect sink. Even if our method couldbe employed to examine far more complex reactive ge-ometries , realized by assembling a large number of spher-ical boundaries of arbitrary sizes and reactivities, for thesake of clarity we shall specialize here to the case of asink of radius σ surrounded by N identical spheres of ra-dius σ = λσ arranged randomly at a fixed distance d .This problem has been tackled recently for λ = 1 and N ≤ ✥(cid:0)✁✥(cid:0)✂✥(cid:0)✄☎ ☎ ✆ ✝ ✁ ✞◆✟✠✡☛☞✌✍✎✏✠☛✑✎✒✟✓✔✑☛✓✑ ❞✴✕q ✖ ✗✘♦✱ ❘✙✚✛ ✜✢✣✤q ✖ ✦✧✘♦✱ ❘✙✚✛ ✜✢✣✤ FIG. 1.
The approximate finite-element calculationscompared with the exact results.
Total flux into a sinkof diameter σ normalized to k S = 4 πDσ (flux into an isolatedsink) in the presence of two spherical screening boundaries ofthe same diameter placed at a fixed distance d from the sinkand forming an angle θ . Light blue and dark orange denotereflecting and absorbing particles, respectively. Symbols arenumerical results of finite-element calculations from Ref. [41].Solid lines are the exact results, obtained by solving eqs. (6)with a relative accuracy of 10 − . N o r m a li z ed r a t e c on s t an t d / 2σ N = 5N = 10N = 50
FIG. 2.
Competitive screening greatly reduces the rateconstant compared to inert obstacles, and is stronglymodulated by the configuration.
Total flux into a sinkof radius σ surrounded by N spherical boundaries of thesame size arranged randomly at distance d (normalized to k S = 4 πDσ ). Symbols denote the exact results (solution ofeqs. (6)), averaged over 100 independent configurations foreach value of d . The shaded bands highlight the regions com-prised between the minimum and maximum rates. For re-flecting screening boundaries, these regions are as small asthe truncation error. The light blue and orange lines areplots of eq. (8) and eq. (10), respectively, with λ = 1. Thearrows flag values of d corresponding to the two configura-tions shown ( N = 50) with the screening spheres made allabsorbing (bottom) and all reflecting (top). In Fig. 1 we compare the FE numerical results withthe exact solution for N = 2. It appears clear that thescreening effect is harder to capture via a FE scheme inthe case of reflecting obstacles than in the presence ofcompetitive sinks. However, our exact approach allowsone to dig much further into this problem and investigateanalytically the screening effect of configurations com-prising a large number of spheres. For example, one canexpand the system (6) in powers of ε ≡ σ/d to derivesimple analytical estimates of the rate constant to thesink (see supplementary material for the detailed calcu-lations). In the case of reflecting obstacles, one gets kk S = 1 − (cid:18) λ N (cid:19) ε − (cid:18) λ N (cid:19) ε + . . . (8)which is independent of the screeninig configuration andlinear in N , as suggested in Ref. [41]. However, we findthat this only holds up to sixth order in ε – it can be seenfrom the expansion that the configuration enters explic-itly successive powers of ε (see supplementary material).On the other hand, a similar procedure in the case of N screening sinks yields kk s = 1 − λN ε + (cid:20) λN + λ N X α,β =1 β = α αβ (cid:21) ε − (cid:20) λ N + λ N X α,β =1 β = α αβ + λ N X α,β,δ =1 β,δ = α αβ Γ αδ (cid:21) ε + . . . (9)where Γ αβ = 2 sin( ω αβ / ω αβ being the angle formedby the sinks α and β with respect to the origin. Eq. (9)makes it very clear that the configuration of competitivereactive boundaries does influence the screening effect onthe central sink. A clear signature of this is also that thecorrections in Eq. (9) alternate in sign. This observationsheds considerable light onto the many-body characterof the rate constant, whose perturbative series is alterna-tively reduced by the diffusive interactions between thescreening boundaries and the sink (shielding ligand fluxfrom it) and increased by the diffusive interactions amongthe screening particles (shielding flux from each other).On the contrary, the screening action of inert obstaclesis largely dominated by the excluded-volume effect, andthus can only yield negative corrections at all orders.Due to its perturbative nature, eq. (9) can be usedto quantify the shielding action of specific 3 D arrange-ments of sinks only for N ε ∝ N/d ≪
1. However, itstill provides a powerful analytical tool to compare dif-ferent geometries, as the perturbative rate is always pro-portional to the true rate (see supplementary material).For example, eq. (9) could be used to design the specialconfigurations that minimize or maximize the screeningeffect on the central sink for given values of N and d .The average shielding action exerted by N equidis-tant sinks can be easily obtained analytically in themonopole approximation, i.e. by keeping only the ℓ = 0and q = 0 terms in eqs. (6). The ensuing reducedsystem can be averaged over different configurations inthe hypothesis of vanishing many-body spatial correla-tions, i.e. by integrating over the probability density P N = Q α = β (sin ω αβ ) / π , with 2 arcsin( σ /d ) ≤ ω αβ ≤ π (excluded-volume constraint between screening sinks).The result is (see supplementary material for the details) (cid:28) kk S (cid:29) = 1 − λε [ N − ( N − − λε )]1 − λε [ N ε − ( N − − λε )] (10)Fig. 2 shows that for λ = 1 eq. (10) provides an ex-tremely good estimate of the configurational averages ofthe exact results at separations greater than a few di-ameters, highlighting the dramatic screening action ofcompetitive reactive boundaries with respect to inert ob-stacles. Furthermore, a simple analysis of the rate fluctu-ations over the configuration ensembles at fixed d allowsone to gauge how sensitive competitive screening is tothe specific 3 D arrangement of the sinks. Remarkably,this analysis reveals stretches between the minimum andmaximum rates for a given value of d as high as 50 % ofthe average (see shaded bands in Fig. 2). More precisely,we remark that the variability associated with differentgeometries is greater ( i ) at short distances and ( ii ) forfew screening particles.Eqs. (8), (9) and (10) and the ensuing argumentsare rather exemplary illustrations of the powerful ana-lytical insight afforded by our general approach. Thecase of small screening sinks, λ <
1, provides a furtherdemonstration of non-trivial effects that are captured byour analytics. It turns out that the function (10) dis-plays a minimum for certain choices of the parameters
N, λ . More precisely, a minimum exists at fixed N for λ ≤ λ ∗ ( N ) ≡ ( √ N + 1 − / (2 N ) <
1, or, alternatively,at fixed λ for N ≤ N ∗ ( λ ) ≡ (1 − λ ) /λ . Fig. 3 providesa clear illustration of this subtle effect.The flux into a large sink features a minimum forscreening configurations of tiny absorbing particles closeto its surface. This is the result of the competition be-tween two effects. When the small particles lie very closeto the surface of the large sink S σ , the latter behaves asan effective isolated sink of size σ eff < σ , absorbing a fluxΦ eff = 4 πDσ eff . k S . Upon increasing the distance d ,the total flux to the screening sinks will increase (their active surfaces get larger and they also get farther apartfrom each other). Now it is clear that the effect of thison the flux to S σ will depend on the size of the screeningparticles. If σ is small enough, σ eff is not much smallerthan σ , so that Φ d = σ + σ ≡ Φ eff is not much smaller thanΦ ∞ = k S . Under these conditions, the flux into the largesink starts decreasing , as the screening ensemble effec-tively steal more and more flux from it. However, uponincreasing d past a critical distance, the small particlescan no longer catch enough ligand flux, so that the fluxto S σ starts increasing, as it should, towards k S .Summarizing, in this paper we have introduced a gen-eral theoretical framework to quantify how the geometry d / σ λ = 0.01 N = 500 N = 1000 λ = 0.1 N o r m a li z ed r a t e c on s t an t d / σ ExactEq. (10), monop. + conf. aver.
FIG. 3.
Making the screening boundaries smallermakes the flux into the sink non-monotonic . Totalflux into a sink of radius σ surrounded by N = 50 smallersinks of radius σ = σ/
10 arranged randomly at a distance d (normalized to k S = 4 πDσ ). The left-most and right-mostcartoons depict two configurations that screen exactly thesame amount of flux, despite being at considerably differ-ent distances ( d/σ = 1 . d/σ = 8). The configurationshown in the middle corresponds to the predicted minimumat d/σ = 1 + p − λ [1 + ( N − λ ] ≈ .
64. The solid line is aplot of formula (10). Each symbol is the average over 250 in-dependent configurations, while the filled band comprises theregion between the minimum and maximum rates. The toppanel illustrates the case of screening by a large number oftiny particles, highlighting the sizeable non-monotonic effect.The curves are plots of eq. (10). and distributed reactivity patterns of the environmentmodulate the rate constant of diffusion-influenced chem-ical reactions. Our method can be used to examine arbi-trary reactive landscapes , made by assembling sphericalboundaries of selected size at given locations in space andendowed with arbitrary surface reactivity. We note thatour method is utterly general, as it can be easily extendedto accommodate for reactive environments realized withmore complex, non-spherical boundaries. The only re-quirement is that one of the 13 coordinate systems inwhich Laplace’s equation is separable be used [43], andthat addition theorems exist for the corresponding ele-mentary solutions [44]. Moreover, our method can beextended to Laplace space [45], so as to work out exactlythe effect of the environment on time-dependent prob-lems. This technique could be employed to shed furtherlight on the intriguing sensitivity of time-dependent ef-fects on initial conditions, which seems to constitute arather generic feature of complex media [23].Finally, we stress that our method can be easily em-ployed to derive approximate closed analytic formulas,that can be used to investigate naturally occurring re-active geometries and to assist in the design of effectiveartificial nano-reactors for different technological appli-cations.The authors wish to thank Sergey Traytak, P. De Los Rios and G. Foffi for insightful comments. F. P. andD. F. acknowledge joint funding from the French CNRS(PICS). [1] S. A. Rice, ed.,
Diffusion-limited reactions , Comprehen-sive Chemical Kinetics, Vol. 25 (Elsevier, Amsterdam,1985).[2] A. Szabo, The Journal of Physical Chemistry , 6929(1989).[3] G. Foffi, A. Pastore, F. Piazza, and P. A. Temussi, Phys-ical Biology , 040301 (2013).[4] J. Sch¨oneberg and F. No´e, PLoS ONE , e74261 EP(2013).[5] K. Seki, M. Wojcik, and M. Tachiya, Physical review. E , 011131 (2012).[6] N. Dorsaz, C. De Michele, F. Piazza, P. De Los Rios, andG. Foffi, Physical Review Letters , 120601 (2010).[7] J. D. Schmit, E. Kamber, and J. Kondev,Physical Review Letters , 218302 (2009).[8] D. Ridgway, G. Broderick, A. Lopez-Campistrous,M. Ru’aini, P. Winter, M. Hamilton, P. Boulanger,A. Kovalenko, and M. J. Ellison, Biophysical Journal , 3748 (2008).[9] D. ben Avraham and S. Havlin, Diffusion and Reactionsin Fractals and Disordered Systems (Cambridge Univer-sity Press, 2000).[10] S. D. Traytak, Chemical Physics , 1 (1995).[11] P. P. Mitra, P. N. Sen, L. M. Schwartz, and P. Le Dous-sal, Phys. Rev. Lett. , 3555 (1992).[12] Y. Lu and M. Ballauff, Progress in Polymer Science ,767 (2011).[13] N. Welsch, A. Wittemann, and M. Ballauff,Journal of Physical Chemistry B , 16039 (2009).[14] K. T. Brady and J. E. Reiner, The Journal of ChemicalPhysics , 074904 (2015).[15] Y. Gao, Y. Chen, X. Ji, X. He, Q. Yin, Z. Zhang, J. Shi,and Y. Li, ACS Nano , ACS Nano , 9788 (2011).[16] J. L. Vivero-Escoto, I. I. Slowing, B. G. Trewyn, andV. S.-Y. Lin, Small , 1952 (2010).[17] M. von Smoluchowski, Physik Z , 557 (1916).[18] O. B´enichou, C. Chevalier, J. Klafter, B. Meyer, andR. Voituriez, Nature Chemistry , 472 (2010).[19] M. F. Shlesinger, G. M. Zaslavsky, and J. Klafter, Nature , 31 (1993).[20] R. Kopelman, Science , 1620 (1988).[21] R. J. Ellis, Current Opinion in Structural Biology ,114 (2001).[22] H. X. Zhou, G. Rivas, and A. P. Minton, Annual reviewof biophysics , 375 (2008).[23] O. B´enichou, C. Chevalier, J. Klafter, B. Meyer, andR. Voituriez, Nature Chemistry , 472 (2010). [24] J. S. Kim and A. Yethiraj, Biophysical journal , 1333(2009).[25] J. A. Dix and A. Verkman, Annual Review of Biophysics , 247 (2008).[26] K. Luby-Phelps, International review of cytology ,189 (2000).[27] J. Kurzidim, D. Coslovich, and G. Kahl, Journal ofPhysics: Condensed Matter , 234122 (2011).[28] B. Nguyen and D. Grebenkov, Journal of StatisticalPhysics , 532 (2010).[29] J. Kurzidim, D. Coslovich, and G. Kahl, Physical ReviewLetters , 138303 (2009).[30] I. L. Novak, P. Kraikivski, and B. M. Slepchenko, Bio-physical Journal , 758 (2009).[31] P. S. Burada, P. H¨anggi, F. Marchesoni, G. Schmid, andP. Talkner, ChemPhysChem , 45 (2009).[32] I. C. Kim and S. Torquato, The Journal of ChemicalPhysics , 1498 (1992).[33] F. Piazza, G. Foffi, and C. De Michele, Journal ofPhysics: Condensed Matter , 245101 (2013).[34] M. von Smoluchowski, Z Phys. Chem. , 129 (1917).[35] P. Debye, Transactions of The Electrochemical Society , 265 (1942).[36] F. C. Collins and G. E. Kimball, Journal of Colloid Sci-ence , 425 (1949).[37] M. M¨uller and N. Lebovka, “Advances in polymer sci-ence,” in Polyelectrolyte Complexes in the Dispersed andSolid State I , Vol. 255 (Springer Berlin Heidelberg, 2014)pp. 57–96.[38] G. Oshanin, M. Moreau, and S. Burlatsky, Advances inColloid and Interface Science , 1 (1994).[39] S. D. Traytak, Chemical Physics Letters , 247 (1992).[40] S. D. Traytak, Chemical Physics , 351 (1995).[41] C. Eun, P. M. Kekenes-Huskey, and J. A. McCammon,Journal of Chemical Physics , 044117 (2013).[42] Here we use the definition Y mℓ ( θ, φ ) = P mℓ (cos θ ) e imφ where P mℓ (cos θ ) are associated Legendre polynomi-als [43].[43] P. M. Morse and H. Feshbach, Methods of theoreticalphysics , Vol. 2 (McGraw-Hill Science/Engineering/Math,1953) pp. 409–415.[44] E. W. Hobson,
The theory of Spherical and EllipsoidalHarmonics (Cambridge University Press, 1931).[45] E. Gordeliy, S. L. Crouch, and S. G. Mogilevskaya, Inter-national Journal for Numerical Methods in Engineering77