A molecular relay race: sequential first-passage events to the terminal reaction centre in a cascade of diffusion controlled processes

We consider a sequential cascade of molecular first-reaction events towards a terminal reaction centre in which each reaction step is controlled by diffusive motion of the particles. The model studied here represents a typical reaction setting encountered in diverse molecular biology systems, in which, e.g., a signal transduction proceeds via a series of consecutive"messengers": the first messenger has to find its respective immobile target site triggering a launch of the second messenger, the second messenger seeks its own target site and provokes a launch of the third messenger and so on, resembling a relay race in human competitions. For such a molecular relay race taking place in infinite one-, two- and three-dimensional systems, we find exact expressions for the probability density function of the time instant of the terminal reaction event, conditioned on preceding successful reaction events on an ordered array of target sites. The obtained expressions pertain to the most general conditions: number of intermediate stages and the corresponding diffusion coefficients, the sizes of the target sites, the distances between them, as well as their reactivities are arbitrary.


Introduction
In 1916 Smoluchowski calculated the typical time it takes a diffusing molecule to locate a reaction centre in a three-dimensional setting [1]. His rate picture of the diffusion limit to a molecular reaction still remains a cornerstone of physical chemistry [2][3][4][5]. Adam and Delbrück [6] as well as Eigen [7] were the first to realise that the diffusive search of a transcription factor protein for its specific binding site on the DNA molecule is more complex than the Smoluchowski picture, involving intermittent threedimensional excursions in the bulk volume and one-dimensional motion along the DNA, thus considerably optimising the search time. This facilitated diffusion picture was its opposite extremity.
We here investigate the probability density function (PDF) of the time instant at which the terminal reaction event in such a sequential reaction cascade takes place. In a basic setting, this random time can be considered as the sum of independent firstpassage times (FPTs) associated with diffusion-reaction steps in a relay race fashion. As a consequence, its PDF is the N-fold convolution of PDFs of each FPT, where N is the number of steps. Using Laplace transform techniques, one can access the PDF through the inverse Laplace transform of the product of the generating functions of individual FPTs. Even though the problem can be considered as formally solved from the mathematical point of view, most properties of such terminal FPTs necessitating a defined sequence of preceding events are yet unknown, to large extent. Particularly, we are interested in the shape of such FPT densities (FPTDs) in systems of different dimensions. What are their asymptotic short-time and long-time behaviours? What are the relevant time scales of the dynamics, and how do they depend on the system parameters? And, importantly, how does the spatial arrangement of targets influence the FPT PDF? Despite of the great biological relevance of such cascade reactions, the answers to these questions remain elusive.
We develop a mathematical framework for the general case with an arbitrary number of intermediate reactions taking place on sequentially labelled target sites located at arbitrary fixed positions in space (figure 1). Moreover we consider both the case of perfect reactions of the messengers with their respective target sites that happen with unit probability upon first mutual encounter, and the case of imperfect reactions. For the latter, an elementary reaction act takes place with a finite probability on each encounter, and thus in general necessitates repeated reaction attempts interspersed with excursions away from the target. Concurrently we consider here a somewhat simplified geometrical setup assuming that the cascade of reactions takes place in an unbounded d-dimensional systems. Such a simplification renders the derivations rather straightforward and the resulting expressions are compact, permitting us to highlight the impact of the intermediate messengers on the terminal reaction event in a transparent way, and also to set an instructive framework for further analysis. Ultimately we obtain explicit and rather simple formulae for the FPTD to the terminal point conditioned by previous passages to an ordered set of positions in space at ordered time instants. In some sense this result can be viewed as a kind of an analogue, formulated in terms of extreme events, of the Wiener measure that defines the probability that the trajectory of a Brownian motion visits some point at time instant t, given that earlier it visited a fixed set of positions at fixed preceding times. At the same time, it is worth stressing that the PDFs of the first-reaction events in unbounded systems do not posses integer moments of any order, while they do in bounded domains, and the latter define important characteristic time-scales. However, in bounded domains the PDFs of the intermediate first-reaction times admit general spectral expansions whose forms depend on the relative distances from the target sites to the domain boundaries (see, e.g., [47][48][49][50][51][52]). This fact makes the analysis rather cumbersome and thus less transparent. A discussion of the cascade reactions for bounded geometries will be presented elsewhere.
This paper is outlined as follows. In section 2 we formulate our model, introduce basic notations and present general results. In section 3 we concentrate on onedimensional systems with perfect or imperfect diffusion-controlled reactions, while in section 4 we present analogous results for two-dimensional systems, followed by a discussion of the behaviour in three-dimensional systems in section 5 for which, in particular, we discuss the dependence of the probability of incomplete reactions as a function of the number of intermediate stages. Lastly in section 6 we conclude with a brief summary of our results and outline some perspectives for further research.

Model and general results
Consider a d-dimensional space, that is infinite in all directions, with an arbitrary number N of immobile target sites (see figure 1 for an example with three target sites). The targets are labelled sequentially by the index j = 1, 2, 3, . . . , N. We refer to them as the first, the second, the third target site, and so on. The target sites are assumed to be spheres with the respective radii r j , and R j are the vectors connecting the origin with their centres. The distance ρ k between the centres of the kth and the (k − 1)th target sites is thus defined by ρ k = |R k − R k−1 |. We use the convention that R 0 = 0 (corresponding to the origin) and hence, ρ 1 = |R 1 | is the distance from the origin to the centre of the first target site.
Suppose next that point-like particle (e.g., the first messenger) launches from the origin at time instant t 0 = 0 and moves diffusively with diffusion coefficient D 1 until it hits (the surface of) the first target site. Then, with probability p 1 , the first messenger binds to this target site and triggers a release of the second point-like particle (e.g., the second messenger) from the centre of the first target site. At this moment, the first messenger is no longer relevant for the further signal relay. § If this reaction event does not complete, with probability 1 − p 1 , the first messenger keeps on diffusing and hitting the surface of the first target site again and again, until a successful binding event followed by the launch of the second messenger. We denote this random time instant as τ 1 . The second messenger diffuses with diffusion coefficient D 2 and seeks the second target site. Once it arrives to this site, a binding event takes place with probability p 2 and the third messenger is released. Similar to before, the fate of the second particle is no longer of interest. In the case of an unsuccessful reaction attempt, with probability 1 − p 2 the diffusive search process is repeated until a successful reaction event takes place. The random duration of the second step is denoted by τ 2 . Such a cascade of diffusion-reaction processes proceeds until the random moment in time when the Nth particle successfully binds to (with probability p N ) the surface of the Nth target site. We are interested here in the form of the probability density function H(t) of the total § Note that some signalling molecules may be inactivated after a certain lifetime [53], while others may have multiple functions, such as "global" regulatory proteins [54].
In chemistry terms successful reactions require the passage of an activation energetic barrier.  Figure 1. Illustration of the molecular relay race on the plane with three circular perfect immobile targets. The first particle starts from the origin (black cross) and moves diffusively with the diffusion coefficient D 1 (blue curve represents an individual simulated trajectory) until it hits the first target (blue); then, the second particle is launched and diffuses (green curve represents an individual simulated trajectory) with the diffusion coefficient D 2 until it finds the second (green) target provoking a release of the third particle. The relay race is terminated when the third particle, diffusing (see the red curve) with the diffusion coefficient D 3 , arrives at the last target (red). Note that for imperfect reactions each particle may need to revisit the corresponding target many times until a reaction event takes place.
time T = τ 1 + τ 2 + . . . + τ N , from the initial release of the first messenger to the terminal reaction at the Nth target site: P{T ∈ (t, t + dt)} = H(t)dt. For simplicity, we will ignore the excluded volume of the targets, i.e., at each stage j, the particle can freely diffuse through all targets i = j. This simplifying assumption allows one to reduce the N-target problem to a sequence of single target problems. Let Ψ j (t), which is a function of the parameters D j , p j , r j and ρ j , denote the PDF of the event that it took the jth particle exactly the time t to bind to the jth target site and, hence, to trigger the release of the (j + 1)th particle, i.e., P{τ j ∈ (t, t+dt)} = Ψ j (t)dt. Correspondingly, letΨ j (p) = L {Ψ j (t)} = ∞ 0 dt e −pt Ψ j (t) denote its Laplace transform, i.e., the generating function of the FPT τ j . As the τ j are considered as independent random variables, then, quite formally, one can write down the desired PDF H(t) as a convolution which implies that the generating function of T is simply the product of the generating functions of τ j ,H Several remarks at this point are in order: (i) Expressions (1) and (2) are fairly general and valid not only for standard diffusive motion but also for anomalous diffusion.
(ii) A natural reset of memory that happens with reactions occurring at intermediate steps within the cascade makes this model applicable to non-Markovian processes with a memory, such as, e.g., a fractional Brownian motion.
(iii) In some settings, a release of the (j + 1)th particle upon arrival of the jth messenger to the corresponding target site does not occur instantaneously but may require some random time δ j (with PDF ψ j (t)) to produce the corresponding effector and to eventually synthesise the (j + 1)th particle. In this case expression (1) is to be appropriately generalised to include such intermediate stages but will still maintain the form of a convolution. As a result expression (2) is to be multiplied by the product N −1 j=2ψ j (p), in whichψ j (p) denotes the generating function of δ j . We note, however, that these purely chemical processes usually do not involve a diffusive transport stage and hence happen at much shorter time scales, which can thus be safely discarded. (iv) Expression (2) implies a simple relation between the nth order moment of the total time T and the moments of the first-reaction times τ j of the intermediate stages within a cascade, if all of them exist, as it happens, e.g., for search processes in bounded domains. This follows straightforwardly from the independence of the first-reaction times of the intermediate stages and can be derived, e.g., by simply differentiating n times both sides of this expression with respect to p and ultimately setting p = 0. Quite trivially, one has where the overbar denotes averaging. Differentiating further, one finds that also the variance of T is just the sum of variances of the first-reaction times τ j and moreover, a cumulant of an arbitrary order of the PDF of T is the sum of cumulants of τ j of the same order. Note, however, that even in bounded domains the moments of first-passage or reaction times are often not representative of actual behaviour [55,56] and stem from anomalously long searching trajectories, corresponding to long-t tails of the PDF. As a consequence, the moments or cumulants alone do not tell us much, if anything about a typical behaviour. We will return to this question at the end of this paper.
(v) There are some subtleties of the behaviour in high-dimensional infinite systems (d ≥ 3), in which the fractal dimension of trajectories (equal to 2 in the case of standard Brownian motion) is smaller than the embedding spatial dimension. In such systems, at each stage there is a finite probability that the reaction does not take place at all because a finite fraction of trajectories travels to infinity and never reaches the corresponding target site. We will address this question in section 5, while sections 3 and 4 focus on the form of the PDF for one-and two-dimensional systems, in which the reaction events happen with probability 1.
In the following, we consider the general case of imperfect reactions such that all (or some of) reaction probabilities p j are less than unity, and thus the particle may have repeated collisions with the target until a successful reaction occurs. The calculation of the corresponding probability density function of the first reaction time T thus amounts to solving the diffusion equation with radiation boundary condition (also called Robin or mixed boundary condition) on the surface of the target site, see, e.g., [57][58][59]. Here the diffusive current through the surface is balanced by the probability density right at the surface, the proportionality factor being the reactivity κ j . This reactivity is related to the reaction probability p j via κ j = p j v/(1 − p j ), where v has a dimension of velocity and equals the thickness of the reaction zone around the target site multiplied by the frequency of reaction attempts upon a contact with the surface, see [60][61][62][63][64] for more details. The reactivity is infinite when the corresponding reaction probability is unity, p j = 1. In turn, the reactivity equals zero when p j = 0, i.e., the reaction is completely inhibited. As the cascade reaction never occurs if at least one of the p j is zero, we assume that all p j (and thus all κ j ) are strictly positive: κ j > 0. We also assume that all diffusion coefficients are finite and strictly positive: 0 < D j < ∞. While targets can be point-like in a one-dimensional setting, one has to consider extended targets in higher dimensions so that r j > 0. We also assume that ρ j > r j .

One-dimensional systems
We start with one-dimensional systems, in which two situations have to be distinguished: (i) in the so-called one-sided case, a particle arriving at the target site is reflected when a reaction attempt is rejected (with probability 1 − p j ) in the direction from which it arrived to this site -in other words, a particle cannot pass through the target site; (ii) in the two-sided case, when the reaction at the target site does not occur, a particle passes through this site and may approach the same target from the opposite direction. In the one-sided case the generating function and the PDF corresponding to the jth stage of the cascade reaction process are given by (see, e.g., [50]) where erfcx(z) = e z 2 erfc(z) is the scaled complementary error function, while erfc(z) is the complementary error function [65].
In order to determine explicitly the PDF of the time of the terminal reaction event in a cascade of imperfect reactions in one-dimensional systems, we first assume that all ratios D j /κ j are different from one another (the case when all of them are equal will be considered below). In this case it is convenient to expand the generating function into partial fractions, with and the naturally emerging characteristic time scale This characteristic time tends to infinity when any of the diffusion coefficients vanishes (or any of the ρ j tends to infinity); in turn, the contribution of the jth stage vanishes when either ρ j = 0 (the target sites overlap) or D j = ∞, as it should be. The Laplace transform can now be readily inverted to give This is an exact expression which holds for arbitrary t and arbitrary sets {D j } and {κ j } (or equivalently {p j }), conditioned by the constraint that all ratios D j /κ j are unequal to each other. Albeit equation (9) defines the PDF explicitly, it contains special functions, and it might be useful to make the behaviour more explicit in the limits t → ∞ and t → 0. In the limit t → ∞, the PDF behaves to leading order as where we have used the identities N j=1 π j = 1 and N j=1 D j π j /κ j = N j=1 D j /κ j . Note that we use here and henceforth the symbol ∼ for asymptotic equivalence. We emphasise two consequences: (i) It is well known (see, e.g., [60,61]) that for diffusion-controlled reactions in low-dimensional systems the asymptotic behaviour of the apparent reaction constants, calculated within a suitably generalised Smoluchowski approach, is independent of the intrinsic chemical rate κ. This is a direct consequence of the fact that in such systems Brownian motion visits each point in the system many times, i.e., oversamples the space. Hence, as time evolves it becomes progressively more difficult to transport a particle diffusively on larger and larger spatial scales.P As a result the diffusive delivery of a particle to the target site poses an increasingly more dominant resistance to the reaction process than the resistance due to the reaction activation barriers embodied in the κ j . In the first-passage time formulation, however, this appears not to be the case -even the long-time behaviour of the PDF of the first reaction event time contains an explicit dependence on all κ j . (ii) For fixed ρ j and κ j , the prefactor in equation (10) is a non-monotonic function of D j , which diverges when either D j → ∞ or D j → 0. The prefactor attains a minimal value when D j = D * j = κ j ρ j . (iii) Finally, the H(t) ≃ t −3/2 scaling is the same as for the Lévy-Smirnov density for first-passage in a semi-infinite domain [60], independent of N.
In the short-time limit, one can expand the exact expression (9) at small t and take into account a rather interesting property of π j : namely, the sums N j=1 π j ( D j /κ j ) n = 0 for any positive integer n < N, and N j=1 π j ( D j /κ j ) N = (−1) N +1 N j=1 D j /κ j . The PDF (9) to leading order behaves as where is another time scale, which is the geometric mean of times D j /κ 2 j . While T characterised exclusively the diffusive steps, T κ accounts for the interplay between diffusion and reaction on each target. Here we see that the exponential cutoff exp(−T /[4t]), similar to the Lévy-Smirnov form, is modulated by the N-dependent power t N −3/2 .
Lastly, we consider the particular case when all D j = D and κ j = κ -the case in which the ratio in the first line in equation (6) cannot be decomposed into elementary fractions. In this special case expression (6) attains the form where R = j ρ j . Inverting the latter expression we find that here the PDF is given by where D −α (z) is the parabolic cylinder function [65]. Note that in both limits t → 0 and t → ∞ the argument of the parabolic cylinder function becomes infinitely large. The asymptotic behaviour of D −α (z) in the large-z limit to leading order follows [65] such that H(t) attains the asymptotic form which is valid for both t → ∞ and t → 0. Keeping only the leading terms we get and Thus, despite the fact that the distributions (9) and (14) pertain to different physical situations and thus have different functional forms, the one in equation (9) defines exactly the same asymptotic behaviour in the limits t → 0 and t → ∞ as expression (14). This can be readily seen by setting D j = D and κ j = κ in equations (10) and (11). We now turn to the two-sided case with imperfect reactions, in which a particle can pass through a target. The expression analogous to expression (6) then reads Inspecting the latter formula and comparing it against expression (6), we conclude that all the above results obtained for the one-sided case are valid upon the simple replacement κ j → κ j /2 , i.e., each intrinsic reactivity just gets reduced by a factor 1/2. At first glance, such a reduction of the reactivity may seem to be somewhat counterintuitive. Indeed, permitting for the reaction to happen from both sides of the target site, one effectively increases the size of the target which, in principle, should ease reaction events. A plausible explanation is that permitting the particle to pass through the target also gives it the possibility to escape from it to both plus and minus infinity, while in the one-sided case it can go to infinity only in one direction. Finally, the above expressions are considerably simplified in the ideal case of perfect reactions when a successful reaction event, resulting in the release of the subsequent messenger, takes place instantaneously upon the first encounter with the surface, p j = 1. Hence all τ j are merely the first-passage times from R j−1 to R j . In this case we obtaiñ with T given by equation (8). These expressions can be either deduced from expression (4) in the limit κ j → ∞ or found directly from the Lévy-Smirnov PDF and generating function of the first-passage time τ j to the perfect target, One can see that the PDF H(t) of the terminal reaction event preserves the form of the individual components, corresponding to the intermediate stages in reaction cascade. In fact expression (20) has exactly the Lévy-Smirnov form in which the characteristic time T includes the parameters representing the intermediate stages.

Two-dimensional systems
The generating function of the first-passage time to a partially reactive circular target in a plane is (see, e.g., [50]) where K ν (z) is the modified Bessel function of the second kind [65]. + Even though the inverse Laplace transform of expression (22) formally yields the PDF, this inversion cannot be performed analytically and cannot even be written in terms of common special functions. One thus has to resort to numerical analysis or concentrate on the asymptotic behaviour of the PDF in the limits t → 0 and t → ∞ (see, e.g., [50,68]). For these purposes it is convenient to represent the inverse Laplace transform as the Bromwich integral and to exploit the regularity property of the K ν (z) in order to change the integration variable and to deform the integration contour in the complex plane. These operations allow one to express Ψ j (t) as the integral over the positive real halfaxis (see [69] for details), Here Im denotes the imaginary part of the expression and where J ν (z) and Y ν (z) are the Bessel functions of the first and second kind, respectively [65]. Following [69] we next generalise expression (24) for the reaction cascade case. + One immediately sees from K 0 (0) = ∞ that the first-passage to a point-like target in two dimensions is impossible.

According to expression (2) the generating function of the first-passage time to the terminal reaction event in a cascade reads
Since the product of ratios of modified Bessel functions of the second kind has no poles one can apply the same technique as above to get the following result: Equation (28) is exact for any t and any values of the parameters, and generalises the Bromwhich-type representation of the PDF of first-reaction times for the problem with a single imperfect target site (see, e.g., [50]) to the case of a cascade of diffusion controlled reactions with an array of targets. The short-time behaviour of the PDF (28) can be obtained by analysing the largep limit of the generating functionH(p) in expression (27). In this limit, K ν (z) ∼ e −z π/(2z) [65] and hencẽ where is a characteristic time scale which only slightly differs from our result (8) above, N ′ is the number of finite κ j , and ′ j denotes the product extending over such values of j for which the corresponding κ j is finite. The inverse Laplace transform yields the short-time asymptotic form When all elementary reaction steps are imperfect, that is, κ j < ∞ for all j, one has N ′ = N and the corresponding short-time behaviour follows from equation (31) by replacing N ′ by N. In the opposite situation when all reaction steps are perfect, i.e., κ j = ∞ for all j, N ′ = 0 and, hence, the product ′ j κ j / D j = 1 which yields Therefore, the short-time tail of H(t) in two-dimensional systems with perfect reactions has the form of the Lévy-Smirnov density, analogous to one-dimensional systems, and differs from that case only by the dependence of the prefactor and of the characteristic time on system parameters. The analysis of the long-time behaviour of H(t) is generally much more subtle (see, e.g., [69]) and here the general integral representation (28) appears to be particularly useful. Following [69] we notice that the long-time limit corresponds to z → 0, for which where γ is the Euler-Mascheroni constant. Substituting these expressions into equation (28) and taking advantage of some standard arguments (see [69]) we find the following long-t asymptotic where We note that the asymptotic form (34) turns out to be remarkably accurate for sufficiently large t, as evidenced by comparison with the numerical inversion of the Laplace transform of expression (27), see figure 2. Expression (34) contains, however, an imaginary part of some rational function, which somewhat obscures its actual time dependence, and it can thus be instructive to calculate this imaginary part explicitly.
To this end, we rewrite B j (t) formally as which implies that B j (t) tends to zero as t → ∞. Expanding then the product in equation (34) in inverse powers of this small parameter and keeping only the leadingorder terms we find which defines the long-time limit explicitly. This asymptotic form is valid for arbitrary values of {ρ j }, {r j }, {κ j }, and {D j }. For perfect reactions, when all κ j = ∞ (and D j > 0), the first term in the numerator vanishes such that the summands becomes proportional to the logarithm of the ratio ρ j /r j and logarithmically weakly depend on the diffusion coefficients.

Three-dimensional systems
In three dimensions the generating function of τ j reads and the PDF has the closed form Rewriting expression (38) formally as we observe that result (40) indeed has exactly the same functional form asΨ j (p) in the one-dimensional case (see equation (4)) if one replaces ρ j by ρ ′ j = ρ j − r j , κ j by κ ′ j = κ j + D j /r j , and multiplies expression (4) by (r j /ρ j )(1 + D j /(κ j r j )) −1 . As a consequence, we can use the results from section 3 for the PDF H(t) and its asymptotic behaviour in the three-dimensional case. For instance, once all D j /κ ′ j are distinct, equation (9) in the three-dimensional case has the form in which and the time scale T is given by equation (30). Moreover π j is determined by equation (7) with κ j replaced by κ ′ j . Finally, one retrieves the same short-time and long-time asymptotic behaviour of H(t) as in the one-dimensional case, where the prefactor is corrected by P react , and where T κ is given by (12) with κ j replaced by κ ′ j . Despite of these similarities the three-dimensional behaviour is crucially different due to the fact that the diffusing messenger can escape to infinity at each diffusionreaction stage of a cascade. Indeed, setting p → 0 in expression (38) one gets the reaction probabilityΨ j (0) = (r j /ρ j )/(1 + D j /(κ j r j )), which is strictly smaller than 1 as soon as r j < ρ j . In other words, the PDF Ψ j (t) is not normalised to unity. As a consequence, the PDF H(t) is also not normalised, Therefore, P react in equation (42) is an important characteristic property which defines, for an arbitrary set of system parameters, the fraction of trajectories in the cascade that eventually arrive to the terminal target site, while 1 − P react is the fraction of trajectories which do not lead to the final reaction event. Since all multipliers in result (42) are less than unity, P react exponentially decreases with the number of stages in the reaction cascade. We also dwell some more on the asymptotic behaviour of expression (42), which harbours some subtle features. One observes that P react vanishes as soon as either of the r j becomes equal to zero meaning that the terminal reaction event does not take place at all. In this limit the corresponding target site becomes point-like, such that it cannot be found by a search process in a three-dimensional space -this stage in a cascade cannot be accomplished successfully. When either of the ρ j becomes infinitely large, P react also vanishes, as it should. For perfect reactions (all κ j → ∞) and D j > 0, the probability of the eventual reaction depends only on the geometrical parameters, and is independent of the diffusion coefficients. In this case, one gets a very simple formula, which is almost identical to (20), For finite κ j , however, P react depends explicitly on all diffusion coefficients. Naturally, P react also vanishes when either of the κ j equals zero, i.e., the corresponding intermediate stage within the cascade is inhibited. The above asymptotic behaviour is physically plausible and intuitively clear. Conversely, in the limit when either of the diffusion coefficients D j vanishes, or, in contrast, tends to infinity, expression (42) produces a somewhat counter-intuitive behaviour: for D j → 0 the corresponding intermediate stage within the cascade cannot be accomplished such that the actual P react has to vanish, while expression (42) does not. In turn, when D j → ∞, meaning that the transport within the corresponding stage takes place instantaneously, the contribution of this stage has to be equal to 1, while expression (42) predicts that P react = 0. The point is that in the derivation of equation (42) we relied on the diffusion equation with the Robin boundary condition, which, in turn, is based on the assumption that the characteristic times of diffusion and reaction attempts are the same. Hence, the diffusion coefficient D j is linked to κ j on the level of the microscopic parameters; in this setting, D j is thus not an independent parameter and cannot be tuned independently of κ j (while the latter, in contrast, can be tuned independently of D j by varying the reaction probability). This observation resolves the apparently controversial behaviour of P react . Extensions of the Robin boundary condition were recently proposed in [64].
In figure 3 we depict our exact result in (41) together with the asymptotic forms (43) and (44), and compare them against the numerical inversion of the Laplace transform (ILT) of the full expression (27). This comparison shows excellent agreement between (41) and the numerical inversion of the Laplace transform, and also indicates the limits in which the asymptotic forms (43) and (44) describe accurately the tails of the PDF.

Conclusion
We established a framework for molecular reaction cascades. Similar to a relay race in human competitions, an incoming signal is here relayed by stepwise reactions to a terminal target site. At each intermediate target, activated through reaction with an incoming diffusive messenger molecule, a new messenger molecule is produced, that in turn needs to locate and react with the next reaction target, until the final messenger molecule reacts with the last target. We obtained the first-reaction time PDF H(t) for the successful terminal reaction event in a cascade of N reaction steps in infinite one, two, and three spatial dimensions along with the asymptotic behaviours in the short and long time limits. We identified the time scale T that is characterised by the target-target distances and the diffusivities of the respective messenger molecules. This specific form of T defines the most typical terminal reaction time and thus mirrors the "geometrycontrol" of direct trajectories from messenger release to the associated reaction target, as observed earlier in one-step reaction settings [48,50] (see also [70] for a "geometrical optics" interpretation). Similar to the Lévy-Smirnov form for the first-passage density we obtain an exponential cutoff at short reaction times, modified by some prefactors, and the long-time behaviour shows the expected power-law scaling of the Lévy-Smirnov density. The dependence of the first-reaction time PDF H(t) on the system parameters depends on the embedding dimension.
We established the reaction cascade framework here in an infinite-space setting without boundaries. This scenario offers relatively simple solutions, allowing us to spotlight the essential dynamic features of reaction cascades. An intriguing question is what will be different if the cascade (or parts thereof) will be placed in finite reaction volumes. Infinite volumes may represent specific cell-to-cell signals, e.g., from bacteria to predator cells such as amoeba or to other bacteria in a growing biofilm. From a broader perspective, one may envisage their utility in robotic search algorithms in which an encounter-propagation may also be of relevance, as well as directed rumour spreading in societal models. Finite volumes, in turn, are relevant in scenarios such as internal signalling in cells. Indeed, in finite volumes, additional time scales spanning considerable time ranges will enter the specific form of H(t) as known for single-step reactions [48,50]. In finite domains, the PDF of the first-reaction time with a single immobile imperfect target generically consists of three temporal domains delimited by characteristic time scales. These domains are: a hump-like region for (relatively) short times, peaked at the most probable reaction time and terminated at the instant of time when a particle first engages with a boundary realising that it moves in a bounded space. This short-time behaviour is followed by a plateau regime, in which all the first-reaction times are equally probable. Mathematically, such a regime emerges due to the gap between the first and the second eigenvalues in the eigenvalue problem which can be distinctly different and show a different dependence on the parameters, in particular, on the chemical reactivity. Finally, the plateau crosses over to an exponential decay, associated with the mean first-reaction time. Now, turning to the typical cellular domains (with their geometrical peculiarities) and signal transduction processes, which necessitate intermediate messengers operating between the subdomains, we may expect that the PDF of the first-reaction time at the terminal target site will have essentially the same functional form. Clearly, in the long-t limit one will find an exponential form with the characteristic time scale which will be close to the mean terminal first-reaction time, equal to the sum of the mean first-reaction times at the intermediate stages within the cascade of reactions, (3). Further on, at relatively short times, at which the "direct" trajectories dominate, one will find either the Lévy-Smirnov t-dependence for perfect reactions or the t-dependence in our equations (11) and (31) for imperfect ones, which we thus expect to be generic. Conversely, a careful analysis for particular, geometrically-relevant settings is required to calculate the amplitudes in both asymptotic laws and their dependence on the system parameters and also the functional form of the eigenvalues, which define the duration of the plateau region. Such an analysis, which is indispensable for understanding the relative importance of different regimes, as indicated, e.g., by the fraction of reaction events happening at each stage, and eventually, the functioning of the signal transduction processes, will be presented elsewhere.
Apart from the above mentioned geometrical particularities of many realistic cascade processes, as exemplified, e.g., by a signal transduction in cellular environments, our analytical description of such molecular relay races can be extended further in several other directions. In particular, upon an interaction of a given messenger with a corresponding target site, often not a single but multiple messengers are released. This is an important aspect and its impact on the first-reaction times and their distribution for a search for a single target site has been rather extensively discussed within the recent years (see, e.g., [71][72][73][74] and [49] and references therein). Second, a target site for a given messenger may not be unique, as we have supposed here, but rather there will exist an array of such sites, each capable to launch a subsequent messenger. For instance, receptors on the cellular membrane may be present in sufficiently big amounts (see, e.g., [14,15]). Next, in biophysical applications, environments in which the particles move are typically extremely complex such that the latter interact with immobile obstacles of different kinds as well as with other mobile particle, which form a dense dynamical background but are not participating in reactions themselves. Such frenetic molecular crowding environments are known to entail a non-Gaussian dynamics, at least at short time scales (see, e.g., [75][76][77]), which will certainly affect the PDF of the terminal reaction time. Moreover, the diffusing molecules may (de-)polymerise on their way [78][79][80] or significantly shift their shape [81], both effects leading to changes of their diffusivities as function of time. Lastly, as interesting perspective consists in studying the effect of the excluded volume of the target sites that was ignored in the present work. This question is apparently not very relevant for the signal transduction processes, for which different messengers operate in different subdomains such that the target sites involved in previous and subsequent intermediate stages have no effect on the dynamics of a given messenger, but may become relevant in infinite systems. In fact, as the particle searches for the jth target, the other targets would typically present inert impermeable obstacles to its motion. When the targets are small and uniformly spread in the medium, their hard-core exclusion effect is expected to be weak, at least at long times. In turn, as the particle starts the search at the jth stage on the surface of the (j − 1)th target, the hindering effect might be notable at short times. While exact solutions are not available in such geometric configurations, asymptotic methods [27,63,82] and semi-analytical approaches [83][84][85][86][87] can be applied.