Abstract
Textbook concepts of diffusionversus kineticcontrol are welldefined for reactionkinetics involving macroscopic concentrations of diffusive reactants that are adequately described by rateconstants—the inverse of the meanfirstpassagetime to the reactionevent. In contradiction, an open important question is whether the meanfirstpassagetime alone is a sufficient measure for biochemical reactions that involve nanomolar reactant concentrations. Here, using a simple yet generic, exactly solvable model we study the effect of diffusion and chemical reactionlimitations on the full reactiontime distribution. We show that it has a complex structure with four distinct regimes delineated by three characteristic time scales spanning a window of several decades. Consequently, the reactiontimes are defocused: no unique timescale characterises the reactionprocess, diffusion and kineticcontrol can no longer be disentangled, and it is imperative to know the full reactiontime distribution. We introduce the concepts of geometry and reactioncontrol, and also quantify each regime by calculating the corresponding reaction depth.
Introduction
Reactions between chemically active molecules in condensed matter systems are typically controlled by two factors: the diffusive search of the species for each other^{1,2,3,4} and the intrinsic reactivity κ associated with the probability that a reaction indeed occurs when the particles collide with each other^{5}. For chemical reactions involving sufficiently high concentrations of particles, which are initially uniformly distributed in the container or reactor such that encounters between reactive species occur more or less uniformly in time, theories based on mean effective reaction rates provide an adequate description of the reaction kinetics^{1,2,3}—apart from some singular and wellknown reaction schemes which exhibit anomalous, fluctuationinduced kinetics under special physical conditions (see, for instance refs. ^{4,6,7,8,9,10,11}). Since the seminal works by Smoluchowski^{12} and Collins and Kimball^{13} a vast number of theoretical advances have scrutinised a combined effect of both ratecontrolling factors on the mean effective rates providing a comprehensive understanding of this effect^{1,2,3,4,14,15,16,17}. In particular, the mean reaction time is the sum of two time scales corresponding to the inverse diffusion coefficient and the inverse intrinsic reactivity (see Eq. (5)), such that the influence of diffusion control and (chemical) rate control are separable^{13}.
For many biochemical reactions, however, the reactive species do not exist in sufficiently abundant amounts to give rise to smooth concentration levels. In contrast, only small numbers of biomolecules, released at certain prescribed positions, are often involved in the reaction process. Indeed, in systems such as the wellstudied Lac and phage lambda repressor proteins only few to few tens of molecules are typically present in a living biological cell, corresponding to nanomolar concentrations. The starting positions of biomolecules can either be rather close to the target or relatively far away. Particularly in the context of the rapid search hypothesis of gene expression, it was shown that the geometric distance between two genes, communicating with each other via signalling proteins, is typically kept short by design in biological cells^{18}, guaranteeing higherthanaverage concentrations of proteins around the target in conjunction with fast and reliable signalling^{19}. Quite generically, many intracellular processes of signalling, regulation, infection, immune reactions, metabolism or transmitter release in neurons are triggered by the arrival of one or few biomolecules to a small spatially localised region^{20,21}. In such cases, it becomes inappropriate to rely on mean rates, and one needs to know the whole distribution of random reaction times, also called the first passage times to a reaction event. Lacking a large number of molecules, reaction times become strongly defocused such that the mean reaction time is no longer representative and the most probable reaction time becomes relevant. We note that even for perfect reactions that occur immediately upon the first encounter between two particles and have thus infinitely large intrinsic reactivity, the mean and the most probable first passage times can differ by orders of magnitude^{22,23} and two first passage events in the same system may be dramatically disparate^{24,25,26}.
For such effectively fewbody reactions, most of the available theoretical effort has been concentrated on the analysis of perfect reactions and hence, on the impact of diffusion control only^{27,28,29}. In particular, in ref. ^{28}, it was argued that for perfect reactions the reaction time density (RTD) can be accurately modelled as
where t_{mean} is the MFPT and q is the contribution of trajectories that arrive to the target site immediately. Conversely, fluctuations of the cycle completion time for enzymatic reactions, in absence of any diffusion stage, have been quantified through the coefficient of variation, γ, of the corresponding distribution function of these times^{30}. Few other works^{31,32,33,34,35} analysed the combined effect of both ratecontrolling factors, but solely for the mean reaction time. These works have shown that the effect of the intrinsic reactivity is certainly significant and even most likely is the dominant factor. The question of the combined influence of both factors on the full distribution of reaction times has been only addressed most recently^{36}, with the focus on the target search kinetics in cylindrical geometries. However, the results of Grebenkov et al.^{36} rely on the socalled selfconsistent approximation^{31} and moreover, have a somewhat cumbersome and thus less practical form. Hence, it is highly desirable to consider particular yet generic examples for which the RTD can be calculated exactly and the results can be presented in a lucid, compact and easy to use form revealing numerous insightful features well beyond the simple approximation in Eq. (1). This is clearly an appealing problem of utmost significance for a conceptual understanding of the kinetics of biochemical reactions.
We here focus on the conceptually and practically relevant question of the influence of the intrinsic chemical reactivity and the initial position of the reacting particles onto the form of the full distribution of reaction times. We demonstrate that when the reactivity is finite and no longer guarantees immediate reaction on mutual encounter, the defocusing of reaction times is strongly enhanced. Remarkably, an extended plateau of the reaction time distribution emerges due to this reaction control, such that the reaction times turn out to be equally probable over several orders of magnitude. A direct consequence of the defocusing is that the contributions of diffusion and rate effects are no longer separable. To distinguish them from the classical concepts of diffusion and kinetic control, we will talk about geometry (initial distance) control and reaction (intrinsic reaction rate) control, keeping in mind that the latter not only specifies the dominant ratecontrolling factor for the MFPT, but affects the shape of the full RTD. An exact solution for the RTD provides us a unique opportunity to derive explicit formulae, for arbitrary initial conditions and arbitrary values of the intrinsic reaction constant κ for several characteristic properties of the distribution such as, e.g., its precise functional forms in different asymptotic regimes, the corresponding crossover times between these kinetic regimes, and also the reaction depths corresponding to these time scales.
Results
Mathematical model
We consider a model involving a pair of reactive molecules: a partially absorbing, immobile target site of radius ρ within a bounded domain of radius R limited by an impenetrable boundary, and a molecule, initially placed at some prescribed position and diffusing with diffusivity D. Once the diffusing particle hits the surface of the target site, it reacts with (binds to) the latter with a finite, intrinsic reaction rate κ. The reflecting outer boundary can mimic an impenetrable cell membrane, the reaction container’s surface or be an effective virtual frontier of the “zone of influence” of the target molecule, separating it from other remotely located target molecules.
Assuming that the domain has a spherical shape and placing the target at the origin of this domain renders the model exactly solvable. We note that although such a geometrical setup is simplified as compared with realistic situations (e.g., the target site is not necessarily located at the centre of the domain^{28,29} or may be attached to some structure which partially screens it^{35,36}), this model captures explicitly two essential ingredients of the reaction process: the diffusive search for the target site and its finite intrinsic reactivity. Importantly, the fact that the model is exactly solvable, permits us to unveil some generic features of the full RTD without resorting to any approximation.
The probability density function H(r, t) of the reaction time t for a particle released a radial distance r−ρ away from the spherical target of radius ρ is calculated using standard tools:^{27,37,38} one first finds the survival probability S(r, t) of a diffusing particle in a radially symmetric situation subject to the zerocurrent boundary condition on the outer boundary of the domain, and the “radiation”, or partiallyreflecting boundary condition^{1,2,3,39}
imposed on the surface of the target site. The proportionality factor κ in Eq. (2) is an intrinsic rate constant (of dimension length/time) whose value shows how readily the particle reacts with the target site upon encounter. When κ = 0 no reaction occurs, while the limit κ = ∞ corresponds to a perfect reaction, when a particle reacts with the target site upon a first encounter. These limiting cases therefore correspond to perfectly reflecting or absorbing boundaries, respectively. The RTD H(r, t) is obtained as the negative derivative of S(r, t) and is valid for arbitrary values of the system parameters. Details of these calculations are presented in the beginning of the Methods section.
Structure of the full distribution of reaction times
The typical shapes of the reaction time density H(r, t) are shown in Fig. 1 for two different release radii r and different values of the dimensionless reactivity κR/D. Note that the parameter κR/D represents a combined effect of two factors: based on the definition of the standard chemical constant K_{on} = 4πρ^{2}κ for a forward reaction and the definition of the socalled Smoluchowski constant K_{S} = 4πDρ we see that κR/D = (K_{on}/K_{S})(R/ρ) and, hence, this is the ratio of the chemical rate and the Smoluchowski rate constant, multiplied by the ratio of the sizes of the domain and of the target site.
We notice that H(r, t) has a much richer structure than the previously proposed simple form in Eq. (1). The RTD consists of four distinct time domains seen in Figs. 1–3: first, a sharp exponential cutoff at short reaction times terminating at the most probable time t_{mp}; second, a region spanning from the most probable reaction time to the crossover time t_{c} in which H(r, t) shows a slow powerlaw decrease; third, an extended plateau region beyond t_{c} which stretches up to the mean reaction time t_{mean}; and fourth, an ultimate longtime exponential cutoff. The shape of the RTD for varying reactivities highlighting the geometrycontrolled Lévy–Smirnov hump and the reactioncontrolled plateau region is our central result. In order to get a deeper understanding of the time scales involved in the reaction process, we also introduce and analyse in the Methods section the forms of two complementary characteristic times: the harmonic mean reaction time t_{harm} = 1/〈1/t〉 and the typical reaction time t_{typ} = t_{0}exp(〈ln(t/t_{0})〉), where the angular brackets denote averaging with respect to the RTD depicted in Figs. 1 and 2, and t_{0} is an arbitrary time scale. Since the logarithm is a slowly varying function, its average value is dominated by the most frequent values of t, while anomalously large/small values corresponding to rare events provide a negligible contribution. Such an averaged value is widely used to estimate a typical behaviour in diverse situations^{40,41}.
Three characteristic time scales
The most probable reaction time, corresponding to the very pronounced maximum, can be calculated explicitly (see the Methods section) and has the approximate form
Interestingly, this simple estimate, which depends only on the diffusion coefficient and the initial distance to the target site, appears to be very robust: t_{mp} indeed shows very little variation with the reactivity κ, as one may infer from Figs. 1 and 3. In the Methods section, we show that when κ decreases from infinity to zero, the value of t_{mp} varies only by a factor of 3. This characteristic time is always strongly skewed towards the left tail of the distribution, that is, to short reaction times: t_{mp} in fact corresponds to particles moving relatively directly from their starting point to the target followed by an immediate reaction and thus generalises the concept of direct, purely geometrycontrolled trajectories^{22} to systems with reaction control. Note that expression (3) is different from the diffusioncontrolled additive contribution proportional to 1/D in the mean reaction time (5).
The second characteristic time scale is the crossover time t_{c} from the humplike Lévy–Smirnov region specific to an unbounded system, to the plateau region. Hence, t_{c} can be interpreted as the time at which a molecule starts to feel the confinement. This can be nicely discerned from comparison with the density H_{∞}(r, t) for the unbounded case (Fig. 2). Thus, reaction times beyond t_{c} correspond to indirect trajectories^{22}. From the result
obtained in the Methods section, we see that t_{c} is independent of the starting point and of the reactivity κ, being entirely dominated by the diffusivity and the difference between the sizes of the domain R and of the target. Writing t_{mp}/t_{c} = π^{2}(r−ρ)^{2}/[12(R−ρ)^{2}], one realises that the crossover time can be comparable with the most probable time (such that the humplike region shrinks), but may also become much larger than the latter when r is close to ρ, as it happens, e.g., when proteins are produced in a close vicinity of a first gene activated at t = 0. In this case, of course, the humplike region will be most pronounced (Fig. 2).
Finally, the onset of the right exponential shoulder at long reaction times coincides with the mean reaction time, as indicated by the arrows in Figs. 1 and 2. The latter is obtained from the Laplace transformed distribution (see the Methods section) and is given by the exact formula
which can be thought of as an analogue of the celebrated Collins–Kimball relation for the apparent reaction rate^{13}. The first term in Eq. (5) is the standard MFPT to a perfectly reactive target and corresponds to the classical notion of diffusioncontrolled rate. The additional contribution to t_{mean} proportional to κ^{−1} accounts for the imperfect reaction with finite reactivity, independent of the particle’s starting point. When t_{mean} is a unique time scale characterising exhaustively well the reaction kinetics, as it happens for reactions with sufficiently high concentrations of reactants, one can indeed distinguish between diffusion or kinetic control. In contradiction, for reactions with nanomolar concentrations of reactive species, the other time scales t_{mp} and t_{c} are equally important and no clearcut separation between diffusion and kinetic control can be made. In the Methods section, we also present an explicit exact expression for the variance of the first reaction time, which permits us to determine the coefficient of variation of the RTD and hence, to quantify its broadness.
Geometry versus reaction control
We emphasise that even for perfect reactions, for which κ = ∞, the mean reaction time is orders of magnitude longer than the most probable reaction time. For imperfect reactions (finite κ values) the mean reaction time becomes even longer, and diverges as 1/κ when κ → 0. The fact that the most probable reaction time is very weakly dependent on κ renders the difference between the most probable and the mean reaction times so much more severe for finite κ. Another remarkable and so far unnoticed feature is that a pronounced plateau develops beyond t_{c}, reflecting an emergent regime of reactioncontrol. This plateau exists even for κ = ∞ (Fig. 1) and becomes increasingly longer with decreasing reactivity κ, implying that over several decades the values of the reaction time become equally probable. Mathematically speaking this plateau appears due to the fact that the smallest eigenvalue of the boundary value problem—the only eigenvalue with an appreciable dependence on κ—disentangles from the remaining eigenvalues. This point is discussed in more detail in the Methods section. Physically, the emergence of the plateau implies that the first passage process to the reaction event becomes even more defocused with decreasing κ, i.e., that the spread of possible reaction times increases significantly. The long spread of reaction times within this plateau region is a consequence of geometrically defocused trajectories exploring the boundary of the reaction volume reinforced by the necessary multiple collisions with the target before a final reaction event due to the reactioncontrol with finite reactivity. An important consequence of the existence of the extended plateau region is that all positive moments of H(r, t), not only the mean reaction time, will be dominated by integration over this region. In other words, the resulting RTD is a concerted effect of geometrycontrol and reactioncontrol.
In Fig. 2, we analyse the effect of the initial distance to the surface of the target site for both perfect and imperfect reactions. The exponential shoulder at long reaction times almost coincides for all cases, especially when the reactivity is finite. This part of the reaction time distribution is dominated by trajectories that equilibrate in the volume before eventual reaction (indirect trajectories^{22}). In contrast, we see a strong variation of the most probable reaction time. The exponential cutoff at short reaction times and the position of the maximum of the distribution is geometrycontrolled, as can be anticipated from the Lévy–Smirnov form for the unbounded problem (see the Methods section): direct trajectories from the initial position to the target need a minimum travel time. For increasing initial distance, the most probable reaction time thus moves to longer times and the relative contribution of the geometrycontrolled fraction of direct trajectories becomes less relevant: instead the particles almost fully equilibrate in the confined volume until they finally react with the target. This reactioncontrol effect is accentuated for decreasing reactivity. We stress that for biological applications, both cases are relevant: shorter initial distances, for instance, are involved when proteins are produced around a first gene activated at time t = 0 and these proteins then need to move to a closeby second gene, here represented by the inner target. This scenario is very similar to the one discussed in reference^{19} as an example for the rapid search hypothesis^{18}. Longer initial distances are relevant when a molecular signal passes the cellular membrane or is produced around a cytoplasmic plasmid, and when these molecules then need to diffuse to the nucleoid region in a bacterial cell or pass the nuclear membrane in an eukaryotic cell. Figure 3 summarises the effects of the finite reactivity and of the distance to the target onto the reaction time distribution in the form of a “heat map”.
Short and longtime behaviour
We now turn to the discussion of the short and longtime tails of H(r, t). The longtime behaviour of the density H(r, t) is determined by the smallest eigenvalue λ_{0} of the Laplace operator. For the spherical domain, one can accurately compute this eigenvalue by solving a trigonometric equation (see the Methods section). When both the target and its reactivity are small one gets λ_{0} ≈ κS_{ρ}/(DV), where the surface area S_{ρ} = 4πρ^{2} of the target and the volume of the domain V ≈ 4πR^{3}/3 are introduced. According to Eq. (5), in this limit t_{mean} ≈ 1/(Dλ_{0}), i.e., the mean reaction time is dominated by multiple returns to the target until the reaction occurs. As the target shrinks (ρ vanishes), the smallest eigenvalue tends to zero. In turn, the other eigenvalues λ_{n}, corresponding to rotationinvariant eigenfunctions of the Laplace operator in the spherical domain, are bounded from below: λ_{n} > π^{2}n^{2}/R^{2} for n = 1,2,…. As a consequence, there is an intermediate range of times, \(1/(D\lambda _1) \ll t \ll 1/(D\lambda _0)\), for which the contribution of all higherorder eigenmodes vanishes, that is, \(e^{  Dt\lambda _n} \ll 1\), whereas the contribution of the lowest eigenmode is almost constant in time, \(e^{  Dt\lambda _0} \approx 1\). This is precisely the reason why the intermediate, plateaulike region emerges, see Fig. 1. Note that this region protrudes over an increasing range of time scales when either the reactivity κ or the target radius ρ decrease, or both. Note also that this intermediate regime corresponds approximately to an exponential law which is often evoked in the context of the first passage statistics to small targets, see, for instance, references^{28,42,43}.
While the smallest eigenvalue determines the plateau and the ultimate exponential cutoff, the shorttime behaviour of the reaction time density H(r, t) is determined by other eigenmodes. Since the limit of a small target \(\left( {\rho \ll R} \right)\) can alternatively be seen as the limit of large domain size, one can use the density H_{∞}(r, t) for diffusion in the exterior of a target, which was first derived by Collins and Kimball^{13},
where \({\mathrm{erfcx}}(x) = e^{x^2}{\mathrm{erfc}}(x)\) is the scaled complementary error function (its derivation is reproduced in the Methods section). As demonstrated in Fig. 1, Eq. (6) fully captures the geometrycontrolled part of the reaction time distribution. In the limit of a perfectly absorbing target, κ → ∞, this expression reduces to
whose normalisation ρ/r ≤ 1 reflects the transient nature of diffusion in three dimensions. One can easily check that the maximum of this Lévy–Smirnovtype density is given precisely by Eq. (3), as intuitively expected.
Approximate form of the full distribution
Combining the short and longtime contributions, we arrive at the following approximate formula for the reaction time density
where t_{mean} ≈ 1/(Dλ_{0}) and
is the hitting probability of the target. The correct normalisation of H(r, t) is ensured by the prefactor in front of the second term. Result (8) is substantially more general than the simple form (1) suggested in ref. ^{28}. The form (8) not only extends expression (1) to the partiallyreactive case, i.e., for arbitrary finite values of κ, but also emphasises and provides an explicit form for the contribution from the humplike region around t_{mp}, which is most relevant for reactions in which the molecule starts close to the target.
Figure 4 illustrates the quality of this approximation, showing that it becomes most accurate when the target radius ρ or reactivity κ are small. One observes that it accurately captures both the maximum, the plateau, and the exponential cutoff of the reaction time distribution. In turn, the transition between the maximum and the plateau region is less sharp than in the exact form. A minor inaccuracy of the approximation (8) is that it reaches a constant—set by the second term—in the shorttime limit while the exact distribution vanishes as t → 0. This feature can be simply removed by multiplying the second term by a Heaviside step function Θ(t − t_{c}) and reevaluating the normalisation constant. But even in the present form approximation (8) provides a remarkably good insight into the behaviour of the first passage dynamics and can thus be used as an efficient and easytohandle fit formula for data analysis or for explicit analytical derivations of followup processes.
Reaction depth
Lastly, we point out that the contributions of the four different regimes separated by the time scales t_{mp}, t_{c}, and t_{mean} can be further quantified by the corresponding reaction depths defining which fraction of trajectories reacted up to a given time. We thus focus now on the cumulative distribution function of reaction times
with the evident property F_{r}(∞) = 1 in a bounded domain in which H(r, t) is normalised, and thus shows explicitly which fraction of trajectories have reacted up to time t. The reaction depth is illustrated in the Methods section. Table 1 summarises the values of the reaction depths of the three characteristic regions of the RTD: the humplike region around t_{mp}, the plateau region, and the exponential tail. We realise that for r/R = 0.2 the least amount of the reaction events happens within the humplike region: it is of order of just 4% for perfect reactions, and this fraction rapidly diminishes upon a decrease of κ. In turn, a much larger amount of the reaction events is collected within the final exponential region. It is typically of order of almost 37%, independently of the value of κ, meaning that for such a value of the ratio r/R roughly one third of all realisations remain unreacted at time t = t_{mean}. However, most of realisations of the reaction events occur within the plateaulike regime—it amounts to roughly 59% for perfect reactions, and becomes even bigger for smaller values of κ. The situation becomes different for a smaller release radius: r/R = 0.02. Here, for perfect reactions the majority of trajectories (49% such that t_{c} is close to the median time) react within the humplike region, while the plateau region and the final exponential tail contribute only 20 and 30%, respectively. Upon lowering κ, the humplike region is no longer representative, and more reaction events take place during the exponential tail (~37%) and the plateaulike regions (~63%), respectively. In conclusion, the plateau region appears to be the most important part of the RTD, which contributes most to the overall number of reaction events, except for the case r/R ≪ 1 and \(\kappa R/D \gg 1\), for which the humplike region becomes the dominant one. Concurrently, this plateau is the region of the strongest defocusing effect, in particular for increased reactioncontrol.
Discussion
Many molecular signalling processes in living biological cells run off at minute concentrations. Similarly in vitro experiments tracking the motion of colloidal particles employ only few particles. Individual first passage events in such situations are defocused, that is, possible reaction times are spread over a vast range comprising orders of magnitude. In particular, this implies that any pair of reaction events will be characterised by highly disparate reaction times. The quantitative description of the reaction time to a target in this scenario therefore cannot simply be based on the mean reaction time. As we showed, the resulting broad distribution of reaction times is due to a conspiracy between geometrycontrol and reactioncontrol effects which cannot be disentangled.
We analysed this phenomenon in detail for a generic spherical geometry, concentrating on several main features. (i) The reaction time density consists of four regions with distinct asymptotic behaviour. (ii) These time regions are separated by three characteristic time scales, which means that there is no unique time scale characterising the kinetic behaviour exhaustively well and the reaction times are defocused. In consequence, the textbook notions of diffusion versus kinetic control, which are appropriate for reactions operating at abundant concentrations, are not applicable in our case. We explicitly determined these times scales and also the associated reaction depths. (iii) A finite reactivity broadens an intermediate regime characterised by an extended plateau region. We showed that the plateau emerges due to a time scale separation of the lowest and the next eigenvalues of the diffusioncontrolling Laplace operator. The fundamental parameter that we found to quantify this intermediate regime is the reactioncontrol represented by the dimensionless reactivity κR/D. A majority of the reaction events occur within this region, except for the case \(r/R \ll 1\) and \(\kappa R/D \gg 1\). In turn, for perfect reactions with a reactant starting very close to the target site the most important part of the RTD is the humplike region which contributes with almost 50% of the reaction events. (iv) The geometry control of the initial particletotarget distance strongly affects the position and the amplitude of the maximum of the reaction time distribution and thus the most probable reaction time. (v) We came up with a simple and thus practical approximative formula for the full reaction time distribution. In particular, we demonstrated that this approximation captures both the most probable and mean reaction times. While the derivation relied on the rotation symmetry of the considered geometric domain, this approximation is expected to be valid in more complex confinements, as long as the target site is far enough from the surrounding outer boundary. Our main conclusion is that reactioncontrol with finite reactivity leads to even stronger reaction time defocusing, stressing the necessity to know the full RTD. This conclusion will serve as a benchmark for the behaviour in geometrically more complex situations^{29} when, e.g., the target site is on the wall or bound to some geometrical structure within the domain, and a fully analytical solution is impossible.
Methods
Exact distribution of reaction times
We consider a diffusion process in a threedimensional domain \({\mathrm{\Omega }} = \left\{ {{\mathbf{x}} \in {\Bbb R}^3:\rho < \left\ {\mathbf{x}} \right\ < R} \right\}\) between two concentric spheres—a small target and a bounding surface of radii ρ and R, respectively. Although the solution of the underlying diffusion problem is well known^{27,38}, we rederive it here for completeness and to highlight several practical points discussed in the main text. In fact, the Laplace transformed probability density function \(\tilde H({\mathbf{x}},p)\) satisfies the modified Helmholtz equation
subject to the boundary conditions
Here, Δ is the Laplace operator, D is the diffusion coefficient, κ is the intrinsic reactivity, and ∂_{n} is the normal derivative directed outward from the domain Ω.
The rotational symmetry of the domain reduces the partial differential Eq. (11) to an ordinary differential equation with respect to the radial coordinate r,
where primes denote derivatives with respect to r. The solution of this equation is
where
with \(\xi = (R  r)\sqrt {p/D}\). It follows that
The mean reaction time is obtained from the Laplacetransformed density as
from which Eq. (5) follows.
In the limit R → ∞, Eqs (14–16) yield
Due to the transient character of threedimensional diffusion, the related distribution is not normalised to unity, but \(\tilde H_\infty (r,p = 0) = (\rho /r)/(1 + D/(\kappa \rho )) < 1\) is the probability of reacting with the target before escaping to infinity. The inverse Laplace transform of Eq. (18) yields Eq. (6). Using the relation \(\tilde S_\infty (r,p) = (1  \tilde H_\infty (r,p))/p\) and Eq. (18), one can also compute the survival probability \(S_\infty (r,t)\) in the time domain
Now we come back to a bounded case with R < ∞. The Laplace inversion of Eq. (14) can be performed by identifying the poles of the function \(\tilde H(r,p)\) in the complex plane \(p \in {\Bbb C}\), that is, by finding the zeros of the function
For convenience, we introduce dimensionless Laplace variable s = (R − ρ)^{2}p/D, so that
where we defined the dimensionless “dilatoriness” parameter μ as
The perfectly reactive target with κ = ∞ corresponds to μ = 0. In other words, for high reactivity κ the value of the dilatoriness μ is small and reactions occur more likely on first encounter, and vice versa. Note that a fully reflecting target with κ = 0 is excluded from our analysis because the reaction time would be infinite. In other words, we always consider 0 ≤ μ < ∞.
The solutions of the equation F(p) = 0 lie on the negative real axis. Setting s = −α^{2}, one gets the trigonometric equation
This equation has infinitely many positive solutions that we denote as α_{n}, with n = 0,1,2,… Since the function on the righthand side has the slope \(\frac{{\rho R + \mu (R  \rho )^2}}{{\rho (R  \rho ) + \mu (R  \rho )^2}} > 1\) near α = 0, the smallest solution α_{0} lies in the interval (0, π/2). More generally, the nth solution lies in the interval (πn, π(n + 1/2)) and tends, for any fixed κ, to the left boundary of the interval as n → ∞. Note that α = 0 (or p = 0) is not a pole of the function \(\tilde H(r,p)\).
Once the poles are identified, we determine the residues by taking the derivative of F(p) at the poles. Applying the theorem of residues to compute the inverse Laplace transform, we finally deduce the exact expression for the probability density H(r, t) of the reaction time for a particle starting at a distance r−ρ from the target,
with
where the expansion coefficients c_{n} are given explicitly by the residues as
Longtime behaviour of the RTD
When either the target radius ρ is small or the dilatoriness parameter μ is large, the slope of the righthand side of Eq. (23) is close to unity and thus the smallest eigenvalue α_{0} is close to zero. Expanding both sides of Eq. (23) into Taylor series one finds the firstorder approximation
In particular, for small target radius, ρ → 0, at fixed dilatoriness μ, we see that \(\alpha _0 \simeq \sqrt 3 (\rho /R)\mu ^{  1/2}\). In turn, when μ → ∞ with fixed ρ,
In both cases, α_{0} is proportional to ρ and inversely proportional to \(\sqrt \mu\). As a consequence, the term with the slowest decay time behaves as
where in the intermediate approximation we ignored terms of order ρ/R and higher, and we introduced the surface area S_{ρ} = 4πρ^{2} of the target and the volume of the domain V ≈ 4πR^{3}/3.
We also note that the approximation c_{0} ≈ 3(ρ/R)^{2}/(μ + 3ρ/(2R)) holds for \(\rho \ll R\), and thus \(c_0/\alpha _0^2 \simeq 1/(1 + 3\rho /(2\mu R))\), i.e., it is close to unity as long as the dilatoriness μ is not too small. Therefore the survival probability can be accurately approximated as \(S(r,t) \simeq {\mathrm{exp}}(  Dt\alpha _0^2/R^2)\) for intermediate and large times. In this case the median reaction time becomes
from which the relation t_{median} ≈ t_{meanx}ln2 follows. This median value is close to the mean reaction time, in which the limit \(\rho \ll R\) has the dominant behaviour as R^{3}/(3κρ^{2}) according to Eq. (5). In turn, the most probable reaction time, which is determined by the higherorder eigenmodes, is orders of magnitude smaller. This behaviour is, however, only present for weakly reactive targets. In contrast, the median time for perfect reactions is usually close to the crossover time t_{c}, while t_{mean} is orders of magnitude larger.
Most probable reaction time
One may deduce from Fig. 1 that the region around the most probable reaction time is well described by the function in (6), which corresponds to the solution in the limit R → ∞. Hence, the most probable reaction time t_{mp} can be obtained with a good accuracy by merely differentiating this function with respect to t and setting the result equal to zero:
where z is defined implicitly as the solution of the following, rather complicated transcendental equation
where erfcx(x) is the scaled complementary error function, and
We denote the solution of this equation as z_{β}. When β tends to 0, a Taylor expansion of the lefthand side of (33) yields z^{2} − 9 + O(β), from which \(z_0 = \sqrt 3\). In the opposite limit β → ∞, one uses the asymptotic behaviour of the function erfcx(x) to get
With some technical efforts, one can prove that z_{β} is a monotonously decreasing function of β (see Fig. 5). We conclude that z_{β} is bounded between \(\sqrt 3\) and 1 so that the most probable time t_{mp} lies between (r − ρ)^{2}/(6D) (for \(\kappa \rho \gg 1\)) and (r − ρ)^{2}/(2D) (for \(\kappa \rho \ll 1\)). In other words, the most probable reaction time shows remarkably weak dependence on the reactivity κ, as illustrated by Fig. 5.
Moments of the reaction time
As we have already remarked in the main text, the positive moments of the RTD of an arbitrary order are dominated by the integration over the plateaulike region such that their values appear close to the onset of the crossover to the final region—the exponential decay of the RTD. The exact values of the positive moments of the random reaction time τ can be accessed directly by a mere differentiation of \(\tilde H(r,p)\) with respect to the Laplace parameter p and subsequently taking the limit p = 0:
For instance, a lengthy but straightforward calculation yields the exact formula for the variance of the reaction time:
from which one also gets the coefficient of variation, \(\gamma = \sqrt {\left\langle {\tau ^2} \right\rangle  \left\langle \tau \right\rangle ^2} /\left\langle \tau \right\rangle\), which characterises fluctuations of the random reaction time τ around its mean, i.e., the effective broadness of the reaction time density. As compared with ref. ^{30}, the expressions Eqs. (5) and (37) permit us to quantify the effect of both ratecontrolling factors.
For a perfectly reactive target, the coefficient of variation diverges as the starting point r approaches ρ, in particular, one gets
when the target is small or the confining domain is large \(\left( {\rho \ll R} \right)\). In turn, for a partially reactive target, the squared coefficient of variation is finite in the limit r → ρ and for a small target reads
The coefficient of variation γ in Eqs (38) and (39) exceeds 1, allowing one to classify this distribution as broad, according to the standard terminology in statistics^{24,25,26}. In both cases, the asymptotic behaviour of γ does not depend on the size of the confining domain, R.
We turn next to the negative order moments of the RTD which are clearly dominated by the region close to the origin and hence, probe the left tail of the distribution. The computation of negative moments (with ν>0) involves integration:
Although this integral is expressed in terms of the explicitly known Laplace transform \(\tilde H(r,p)\) from Eq. (14), its analytical evaluation does not seem to be feasible.
In turn, the integral takes a more tractable form in the limit R → ∞ corresponding to diffusion in the exterior of a partially reactive target of radius ρ. Due to the transient character of diffusion in three dimensions, the probability density H_{∞}(r, t) is not normalised to 1 as the molecule can escape to infinity. The integral of the density H(r, t) yields thus the probability of reacting at the target:
The negative order moments of the renormalised density H_{∞}(r, t)/q are
where β was defined in (34). In the limit κ → ∞, one finds
While the mean reaction time diverges for the exterior problem, the negative order moments are welldefined and can thus characterise the reaction process. In particular, the harmonic mean reaction time, defined as
is deduced from (42) for ν = 1:
where \({\mathrm{Ei}}_1(z) = {\int}_1^\infty dx{\kern 1pt} e^{  zx}/x\) is the exponential integral. The dependence of the harmonic mean on the reactivity κ is fully captured via β. In the limit κ → ∞, this mean approaches
and is thus of the order of the most probable time, representing the relevant time scale of the problem. In the opposite limit κ → 0, β approaches a constant, and the harmonic mean reaction time also reaches a constant. One can check that t_{harm} monotonously decreases as β (or κ) grows.
Figure 6 illustrates by dashed lines the behaviour of the function in (45), in particular, its approach to the limiting expression (46) as κ increases. One can appreciate a very weak dependence of the harmonic mean reaction time for the exterior problem on the reactivity κ. We also show the harmonic mean reaction time in the concentric domain, obtained by a numerical integration in Eq. (40) with ν = 1. This mean significantly depends on κ and behaves as 1/κ for small κ. Given that the probability density H(r, t) for the concentric domain can be accurately approximated by H_{∞}(r, t) at small times (see Eq. (8)), the harmonic mean reaction time for the concentric domain can be approximated by the expression in (40), multiplied by the reaction probability q. This approximation, shown by solid lines, turns out to be remarkably accurate when the target radius ρ is small as compared to the radius R of the confining domain. We can also conclude that the significant variations of t_{harm} with κ for the concentric domain come from those of q with κ.
Finally, we consider the time scale
(where t_{0} is an arbitrary time scale), based on the mean logarithm of the reaction time—an important characteristic of the reaction process, which emphasises the typical values of t, i.e., values observed in most of experiments. Indeed, the logarithm is a slowly varying function and its average is supported by the most frequently encountered values of t with the rare anomalously long or shortreaction times being effectively filtered out. The estimates based on t_{typ} are widely used in the analysis of stochastic reactiondiffusion or transport process in random environments (see, e.g., refs. ^{40,41} and references therein). Such an averaged value can be formally computed as
where γ ≈ 0.5772… is the Euler constant, from which
where u_{n}(r) are given by (26).
To get a more explicit dependence on the initial radius r, one can again consider the exterior problem (R = ∞). Rewriting Eq. (42) as
in order to get a Taylor expansion as ν → 0, one finds
where the expectation is computed with respect to the renormalised density H_{∞}(r, t)/q. We obtain thus the logarithmic mean time
In the limit κ → ∞, e^{β}Ei_{1}(β) vanishes as 1/β, so that for a perfectly reactive target one gets
which signifies that in the limit κ = ∞ the logarithmic mean time is comparable to the most probable reaction time t_{mp}.
Figure 7 shows the logarithmic mean reaction time, t_{typ}, as a function of the dimensionless reactivity κR/D. As for the harmonic mean in Fig. 6, the results for a bounded concentric domain (R = 1) and for the exterior problem (R = ∞) differ significantly. The particular definition of the logarithmic time does not allow one to easily renormalise t_{typ} for the exterior domain to get an approximation for the bounded domain.
Finally, Fig. 8 compares several mean reaction times for the concentric domain. One can see that the behaviour of the median, the harmonic and the logarithmic means resembles that of the conventional (arithmetic) mean FPT. In particular, all these means behave as 1/κ at small κ, indicating that the reaction is limited by the kinetics. Only the most probable FPT exhibits a very different behaviour and shows almost no depedence on the reactivity κ, as discussed above.
Reaction depth
The reaction depth in Eq. (10) is shown in Fig. 9. Note first that the reaction depths corresponding to the shortest characteristic time t_{mp} are evidently the shortest, amounting to only about 4% for perfect reactions and r close to ρ. For finite κ or for starting points further away from the target, the reaction depth F_{r}(t_{mp}) diminishes. In turn, in all cases the reaction depth connected to the intermediate plateau is dominant, increasingly so due to the reactioncontrol at lower reactivities.
Data and code availability
All figures have been prepared by means of Matlab software. The plotted quantities have been computed by explicit formulas provided in the paper by using custom routines for Matlab software. While the explicit form makes these numerical computations straightforward, custom routines are available from the corresponding author upon request.
References
 1.
Calef, D. F. & Deutch, J. M. Diffusioncontrolled reactions. Ann. Rev. Phys. Chem. 34, 493–524 (1983).
 2.
Weiss, G. H. Overview of theoretical models for reaction rates. J. Stat. Phys. 42, 3–36 (1986).
 3.
Shoup, D. & Szabo, A. Role of diffusion in ligand binding to macromolecules and cellbound receptors. Biophys. J. 40, 33–39 (1982).
 4.
Lindenberg, K., Oshanin, G. & Tachiya, M. Reaction kinetics beyond the textbook: fluctuations, many particle effects and anomalous dynamics. J. Phys.: Condens. Matter 19, 1 (2007).
 5.
Hänggi, P., Talkner, P. & Borkovec, M. Reactionrate theory: fifty years after Kramers. Rev. Mod. Phys. 62, 251–341 (1990).
 6.
Krapivsky, P. L., Redner, S., & BenNaim, E. A Kinetic View of Statistical Physics (Cambridge University Press, Cambridge, 2010).
 7.
Oshanin, G. & Burlatsky, S. F. Fluctuationinduced kinetics of reversible coagulation. J. Phys. A: Math. Gen. 22, L973–L976 (1989).
 8.
Oshanin, G., Ovchinnikov, A. A. & Burlatsky, S. F. Fluctuationinduced kinetics of reversible reactions. J. Phys. A: Math. Gen. 22, L977–L982 (1989).
 9.
Yuste, S. B., Oshanin, G., Lindenberg, K., Bénichou, O. & Klafter, J. Survival probability of a particle in a sea of mobile traps: A tale of tails. Phys. Rev. E 78, 021105 (2008).
 10.
Grebenkov, D. S. Searching for partially reactive sites: analytical results for spherical targets. J. Chem. Phys. 132, 034104 (2010).
 11.
Grebenkov, D. S. Subdiffusion in a bounded domain with a partially absorbingreflecting boundary. Phys. Rev. E 81, 021128 (2010).
 12.
Smoluchowski, M. Versuch einer matematischen theorie der koagulationskinetik kolloider lösungen. Z. Phys. Chem. 92, 129–168 (1917).
 13.
Collins, F. C. & Kimball, G. E. Diffusioncontrolled reaction rates. J. Colloid Sci. 4, 425–437 (1949).
 14.
Sapoval, B. General formulation of Laplacian transfer across irregular surfaces. Phys. Rev. Lett. 73, 3314–3316 (1994).
 15.
Grebenkov, D. S. Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena, in “Focus on Probability Theory” (ed. Velle L. R.) 135–169 (Nova Science Publishers, New York, 2006).
 16.
Holcman, D. & Schuss, Z. Control of flux by narrow passages and hidden targets in cellular biology. Phys. Progr. Rep. 76, 074601 (2013).
 17.
GudowskaNowak, E., Lindenberg, K. & Metzler, R. Marian Smoluchowski’s 1916 paper – a century of inspiration. J. Phys. A: Math. Theor. 50, 380301 (2017).
 18.
Kolesov, G., Wunderlich, Z., Laikova, O. N., Gelfand, M. S. & Mirny, L. A. How gene order is influenced by the biophysics of transcription regulation. Proc. Natl. Acad. Sci. USA 104, 13948–13953 (2007).
 19.
Pulkkinen, O. & Metzler, R. Distance matters: the impact of gene proximity in bacterial gene regulation. Phys. Rev. Lett. 110, 198101 (2013).
 20.
Alberts, B. et al. Molecular Biology of the Cell. (Garland Science, New York, NY, 2014).
 21.
Snustad, D. P. & Simmons, M. J. Principles of Genetics. (Wiley, New York, 2000).
 22.
Godec, A. & Metzler, R. Universal proximity effect in target search kinetics in the fewencounter limit. Phys. Rev. X 6, 041037 (2016).
 23.
Godec, A. & Metzler, R. First passage time distribution in heterogeneity controlled kinetics: going beyond the mean first passage time. Sci. Rep. 6, 20349 (2016).
 24.
MejíaMonasterio, C., Oshanin, G., & Schehr, G. First passages for a search by a swarm of independent random searchers, J. Stat. Mech. P06022 (2011). https://doi.org/10.1088/17425468/2011/06/P06022
 25.
Mattos, T., MejíaMonasterio, C., Metzler, R. & Oshanin, G. First passages in bounded domains: When is the mean first passage time meaningful? Phys. Rev. E 86, 031143 (2012).
 26.
Mattos, T., MejíaMonasterio, C., Metzler, R., Oshanin, G., & Schehr, G. Trajectorytotrajectory fluctuations in firstpassage phenomena in bounded domains, in “FirstPassage Phenomena and Their Applications”, (eds. Metzler R., Oshanin G., & Redner S.) 203–225. (World Scientific, Singapore, 2014).
 27.
Redner, S. A Guide to First Passage Processes. (Cambridge University press, Cambridge, 2001).
 28.
Bénichou, O., Chevalier, C., Klafter, J., Meyer, B. & Voituriez, R. Geometrycontrolled kinetics. Nat. Chem. 2, 472–477 (2010).
 29.
Bénichou, O. & Voituriez, R. From firstpassage times of random walks in confinement to geometrycontrolled kinetics. Phys. Rep. 539, 225–284 (2014).
 30.
Moffitt, J. R., Chemla, Y. R. & Bustamante, C. Mechanistic constraints from the substrate concentration dependence of enzymatic fluctuations. Proc. Natl Acad. Sci. USA 107, 15739–15744 (2010).
 31.
Shoup, D., Lipari, G. & Szabo, A. Diffusioncontrolled bimolecular reaction rates. Biophys. J. 36, 697–714 (1981).
 32.
von Hippel, P. H. & Berg, O. Facilitated target location in biological systems. J. Biol. Chem. 264, 675–678 (1989).
 33.
van Zon, J. S. & ten Wolde, P. R. Simulating biochemical networks at the particle level and in time and space: green’s function reaction dynamics. Phys. Rev. Lett. 94, 128103 (2005).
 34.
Grebenkov, D. S. & Oshanin, G. Diffusive escape through a narrow opening: new insights into a classic problem. Phys. Chem. Chem. Phys. 19, 2723–2739 (2017).
 35.
Grebenkov, D. S., Metzler, R. & Oshanin, G. Effects of the target aspect ratio and intrinsic reactivity onto diffusive search in bounded domains. New J. Phys. 19, 103025 (2017).
 36.
Grebenkov, D. S., Metzler, R. & Oshanin, G. Towards a full quantitative description of singlemolecule reaction kinetics in biological cells. Phys. Chem. Chem. Phys. 20, 16393–16401 (2018).
 37.
Gardiner, C. W. Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences (Berlin, Springer, 1985).
 38.
Carslaw, H. S. & Jaeger, J. C. Conduction of Heat in Solids (Oxford University Press, Oxford, UK, 1959).
 39.
Oshanin, G., Popescu, M. & Dietrich, S. Active colloids in the context of chemical kinetics. J. Phys. A: Math. Theor. 50, 134001 (2017).
 40.
Evans, M. R. & Majumdar, S. N. Diffusion with stochastic resetting. Phys. Rev. Lett. 106, 160601 (2011).
 41.
Dean, D. S., Gupta, S., Oshanin, G., Schehr, G. & Rosso, A. Diffusion in periodic, correlated random forcing landscapes. J. Phys. A: Math. Theor. 47, 372001 (2014).
 42.
Meyer, B., Chevalier, C., Voituriez, R. & Bénichou, O. Universality classes of firstpassagetime distribution in confined media. Phys. Rev. E 83, 051116 (2011).
 43.
Isaacson, S. A. & Newby, J. Uniform asymptotic approximation of diffusion to a small target. Phys. Rev. E 88, 012820 (2013).
Acknowledgements
D.S.G. acknowledges the support under Grant No. ANR13JSV5000601 of the French National Research Agency. RM acknowledges funding from Deutsche Forschungsgemeinschaft (project ME 1535/61) and from the Foundation for Polish Science within a Humboldt Polish Honorary Research Scholarship. The authors acknowledge the support of Deutsche Forschungsgemeinschaft (German Research Foundation) and Open Access Publication Fund of Potsdam University.
Author information
Affiliations
Contributions
D.S.G., R.M. and G.O. formulated the problem. D.S.G. performed the mathematical derivations and prepared the figures. D.S.G., R.M. and G.O. analysed the obtained results and wrote the paper.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Grebenkov, D.S., Metzler, R. & Oshanin, G. Strong defocusing of molecular reaction times results from an interplay of geometry and reaction control. Commun Chem 1, 96 (2018). https://doi.org/10.1038/s420040180096x
Received:
Accepted:
Published:
Further reading

A firstpassage approach to diffusioninfluenced reversible binding and its insights into nanoscale signaling at the presynapse
Scientific Reports (2021)

A Probabilistic Approach to Extreme Statistics of Brownian Escape Times in Dimensions 1, 2, and 3
Journal of Nonlinear Science (2020)

Distribution of extreme first passage times of diffusion
Journal of Mathematical Biology (2020)

Theory of transcription bursting: stochasticity in the transcription rates
Journal of Mathematical Chemistry (2020)

Burst statistics in an early biofilm quorum sensing model: the role of spatial colonygrowth heterogeneity
Scientific Reports (2019)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.