Strong defocusing of molecular reaction times results from an interplay of geometry and reaction control

Text-book concepts of diffusion- versus kinetic-control are well-defined for reaction-kinetics involving macroscopic concentrations of diffusive reactants that are adequately described by rate-constants -- the inverse of the mean-first-passage-time to the reaction-event. In contradistinction, an open important question is whether the mean-first-passage-time 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 conspiratory effect of diffusion and chemical reaction-limitations on the full reaction-time distribution. We show that it has a complex structure with four distinct regimes delimited by three characteristic time scales spanning a window of several decades. Consequently, the reaction-times are defocused: no unique time-scale characterises the reaction-process, diffusion- and kinetic-control can no longer be disentangled, and it is imperative to know the full reaction-time distribution. We introduce the concepts of geometry- and reaction-control, and also quantify each regime by calculating the corresponding reaction depth.

R eactions 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-3 -apart from some singular and well-known reaction schemes which exhibit anomalous, fluctuation-induced 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 rate-controlling 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 well-studied 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 higher-than-average 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 few-body 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 HðtÞ % qδðtÞ þ ð1 À qÞexpðÀt=t mean Þ=t mean ; 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 rate-controlling 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 so-called self-consistent 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 rate-controlling 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 zero-current boundary condition on the outer boundary of the domain, and the "radiation", or partially-reflecting boundary condition [1][2][3]39 D ∂Sðr; tÞ ∂r 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 so-called 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 cut-off 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 power-law decrease; third, an extended plateau region beyond t c which stretches up to the mean reaction time t mean ; and fourth, an ultimate long-time exponential cut-off. The shape of the RTD for varying reactivities highlighting the geometry-controlled 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  Fig. 1 Reaction control. Reaction time density H(r, t) for a reaction on an inner target of radius ρ/R = 0.01, with starting point a r/R = 0.2 and b r/R = 0.02 for four progressively decreasing (from top to bottom) values of the dimensionless reactivity κ′ = κR/D indicated in the plot. Note that κ′ includes R and D such that smaller values of κ′ can also be achieved at a fixed κ upon lowering R or by increasing the values of D. The coloured vertical arrows indicate the mean reaction times for these cases. The vertical black dashed line indicates the crossover time t c = 2(R − ρ) 2 /(Dπ 2 ) above which the contribution of higher order Laplacian eigenmodes become negligible. This characteristic time marks the end of the hump-like region (Lévy-Smirnov region specific to an unbounded system, see below and the Methods section for more details) and indicates the crossover to a plateau region with equiprobable realisations of the reaction times. This plateau region spans a considerable window of reaction times, especially for lower reactivity values. Thin coloured lines show the reaction time density H ∞ (r, t) from Eq. (6) for the unbounded case (R → ∞). Length and time units are fixed by setting R = 1 and R 2 /D = 1. Note the extremely broad range of relevant reaction times (the horizontal axis) spanning over 12 orders of magnitude for the panel b. Coloured bar-codes c, d indicate the cumulative depths corresponding to four considered values of κ′ in decreasing order from top to bottom. Each bar-code is split into ten regions of alternating brightness, representing ten 10%-quantiles of the distribution (e.g., the first dark blue region of the top bar-code in panel c indicates that 10% of reaction events occur till Dt=R 2 ' 1) calculated explicitly (see the Methods section) and has the approximate form t mp % ðr À ρÞ 2 =ð6DÞ: 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 geometry-controlled trajectories 22 to systems with reaction control. Note that expression (3) is different from the diffusion-controlled 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 hump-like Lévy-Smirnov region specific to an Reaction versus geometry control. Impact of the finite reactivity (reaction control, a) and of the distance to the target (geometry control, b) onto the reaction time density shown as a "heatmap", in which the value of the reaction time density (in arbitrary units) is determined by the colour code. Blue and white lines respectively show the mean and the most probable reaction times that differ by orders of magnitude. The grey vertical line indicates the crossover time t c that does depend neither on reactivity, nor on distance. a When the reactivity κ decreases (with r/R = 0.2 being fixed), the distribution becomes much broader and extends towards longer reaction times. b When the distance to the target decreases (with κR/D = 1 fixed), the most probable reaction time shifts to the left, whereas the mean reaction time remains constant 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 hump-like 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 hump-like 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 clear-cut 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 reaction-control. 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 reaction-control 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 geometry-control and reaction-control. 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 cut-off at short reaction times and the position of the maximum of the distribution is geometry-controlled, 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 geometry-controlled 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 reaction-control 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 close-by 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 long-time behaviour. We now turn to the discussion of the short-and long-time tails of H(r, t). The long-time 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 rotation-invariant 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λ 1 Þ ( t ( 1=ðDλ 0 Þ, for which the contribution of all higher-order eigenmodes vanishes, that is, e ÀDtλ n ( 1, whereas the contribution of the lowest eigenmode is almost constant in time, e ÀDtλ 0 % 1. This is precisely the reason why the intermediate, plateau-like 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 cut-off, the short-time behaviour of the reaction time density H(r, t) is determined by other eigenmodes. Since the limit of a small target ρ ( R ð Þ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 erfcxðxÞ ¼ e x 2 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-Smirnov-type density is given precisely by Eq. (3), as intuitively expected.
Approximate form of the full distribution. Combining the short-and long-time contributions, we arrive at the following approximate formula for the reaction time density Hðr; tÞ % H 1 ðr; tÞ þ ð1 À qÞ e Àt=t mean t mean ; ð8Þ 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 partially-reactive case, i.e., for arbitrary finite values of κ, but also emphasises and provides an explicit form for the contribution from the hump-like 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 cut-off 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 short-time 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 re-evaluating 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 easyto-handle fit formula for data analysis or for explicit analytical derivations of follow-up 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  Table 1 summarises the values of the reaction depths of the three characteristic regions of the RTD: the hump-like 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 hump-like 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 plateau-like 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 hump-like region, while the plateau region and the final exponential tail contribute only 20 and 30%, respectively. Upon lowering κ, the hump-like region is no longer representative, and more reaction events take place during the exponential tail (~37%) and the plateau-like 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 κR=D ) 1, for which the hump-like region becomes the dominant one. Concurrently, this plateau is the region of the strongest defocusing effect, in particular for increased reaction-control.

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 geometry-control and reaction-control 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 diffusion-controlling Laplace operator. The fundamental parameter that we found to quantify this intermediate regime is the reaction-control represented by the dimensionless reactivity κR/D. A majority of the reaction events occur within this region, except for the case r=R ( 1 and κR=D ) 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 hump-like region which contributes with almost 50% of the reaction events. (iv) The geometry control of the initial particle-to-target 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 reaction-control 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 Ω ¼ x 2 R 3 : ρ< x k k<R È É between two concentric spheresa 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 functionHð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, H′′ðr; pÞ þ 2 rH ′ðr; pÞ À p DH ðr; pÞ ¼ 0; H′ðR; pÞ ¼ 0;  Relative weights (in percents) of three characteristic regions of the reaction time density for ρ/R = 0.01: the hump-like region around the most probable reaction time tmp, extending from 0 till tc = 2(R−ρ) 2 /(π 2 D) (and thus merging two subregions discussed in the text: the exponential tail left to tmp and the power-law decay right to tmp); the plateau-like region stretching from tc to the mean reaction time tmean; and the exponential tail which persists beyond t = tmean. Two starting points r/R and four values of dimensionless reactivity κ′ = κR/D are used, corresponding to Fig. 1 where with ξ ¼ ðR À rÞ ffiffiffiffiffiffiffiffi ffi p=D p . It follows that The mean reaction time is obtained from the Laplace-transformed density as from which Eq. (5) follows.
In the limit R → ∞, Eqs (14-16) yield Due to the transient character of three-dimensional diffusion, the related distribution is not normalised to unity, butH 1 ðr; p ¼ 0Þ ¼ ðρ=rÞ=ð1 þ D=ðκρÞÞ<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 relationS 1 ðr; pÞ ¼ ð1 ÀH 1 ðr; pÞÞ=p and Eq. (18), one can also compute the survival probability S 1 ð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 functionHðr; pÞ in the complex plane p 2 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 right-hand side has the slope ρRþμðRÀρÞ 2 ρðRÀρÞþμðRÀρÞ 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 functionHð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 λ n ¼ α 2 n =ðR À ρÞ 2 ; ð25Þ u n ðrÞ ¼ c n D ðR À ρÞ 2 Rα n cos α n RÀr RÀρ À ðR À ρÞsin α n RÀr RÀρ where the expansion coefficients c n are given explicitly by the residues as c n ¼ 2ρ 2 α 2 n ðρR þ μðR 2 þ ρ 2 ÞÞα n sin α n þ ρðμRα 2 n À ρÞ cos α n : Long-time behaviour of the RTD. When either the target radius ρ is small or the dilatoriness parameter μ is large, the slope of the right-hand 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 first-order approximation In particular, for small target radius, ρ → 0, at fixed dilatoriness μ, we see that α 0 ' ffiffi ffi 3 p ðρ=RÞμ À1=2 . In turn, when μ → ∞ with fixed ρ, In both cases, α 0 is proportional to ρ and inversely proportional to ffiffi ffi μ p . 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 ρ ( R, and thus c 0 =α 2 0 ' 1=ð1 þ 3ρ=ð2μ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Þ ' expðÀDtα 2 0 =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 ρ ( 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 higher-order 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 left-hand side of (33) yields z 2 − 9 + O(β), from which z 0 ¼ ffiffi ffi 3 p . 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 ffiffi ffi 3 p and 1 so that the most probable time t mp lies between (r − ρ) 2 /(6D) (for κρ ) 1) and (r − ρ) 2 /(2D) (for κρ ( 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 plateau-like 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 ofHð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, γ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi τ 2 h i À τ h i 2 q = τ h i, 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 ρ ( R ð Þ . 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 transformHð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 well-defined and can thus characterise the reaction process. In particular, the harmonic mean reaction time, defined as is deduced from (42) for ν = 1: where Ei 1 ðzÞ ¼ R 1 1 dx 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 reaction-diffusion or transport process in random environments (see, e.g., refs. 40 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 reaction-control at lower reactivities. Dt /R 2 Dt /R 2 R /D = 1 R /D = 0.1 Fig. 9 Cumulative distribution function of reaction times, F r (t). It is evaluated for the reaction on an inner target of radius ρ/R = 0.01, with the starting point a r/R = 0.2 and b r/R = 0.02 and varying reactivity κ. Symbols indicate the relevant characteristic times: most probable time t mp (asterisks), harmonic mean t harm (squares), logarithmic mean t log (diamonds), median t median (triangles), and mean time (circles). Note that some most probable times are not seen at this scale