Distribution of first-reaction times with target regions on boundaries of shell-like domains

We study the probability density function (PDF) of the first-reaction times between a diffusive ligand and a membrane-bound, immobile imperfect target region in a restricted"onion-shell"geometry bounded by two nested membranes of arbitrary shapes. For such a setting, encountered in diverse molecular signal transduction pathways or in the narrow escape problem with additional steric constraints, we derive an exact spectral form of the PDF, as well as present its approximate form calculated by help of the so-called self-consistent approximation. For a particular case when the nested domains are concentric spheres, we get a fully explicit form of the approximated PDF, assess the accuracy of this approximation, and discuss various facets of the obtained distributions. Our results can be straightforwardly applied to describe the PDF of the terminal reaction event in multi-stage signal transduction processes.


Introduction
A completed reaction event between a diffusive particle and a specific target site often represents an intermediate yet crucial step in diverse biochemical and biophysical processes. In many realistic situations a particle diffuses in a shell-like region delimited by impermeable outer and inner boundaries and reacts with an immobile target region (e.g., a catalytic site; in the remainder, we simply refer to "the target") placed on either of the boundaries. In some applications, this target is located on the inner boundary ( figure 1(a)). This is a common situation in chemoreception processes [1,2] as well as, more generally, in cellular signal transduction pathways [3][4][5]. Here the shell-like domain can be an extracellular medium and the inner boundary represents the (outer) plasma membrane of a cell. In such a setting the particle is commonly referred to as the "ligand" or, in the literature on signal transduction, as the first "messenger". The target is a receptor that undergoes a conformational change when the first messenger binds to it, stimulating then a synthesis of the second messenger which moves inside the cell itself, i.e., within the inner domain. Similarly the particle may cross the cell wall through Figure 1. Sketch of our geometrical setup. A domain Ω 2 with an impermeable boundary ∂Ω 2 is nesting a smaller domain Ω 1 enclosed by an impermeable boundary ∂Ω 1 . A particle, whose starting position is indicated by the blue point, diffuses within the shell-like domain Ω = Ω 2 \Ω 1 delimited by the two boundaries and seeks an immobile target (drawn here as an interval in red) placed on either the inner boundary (a) or the outer boundary (b). membrane pores. In a different scenario, the shell-like domain can be the intracellular medium (cytoplasm). Then the outer boundary is the cellular membrane and the particle can be, e.g., the second messenger, which searches diffusively for a specific target on the inner boundary, e.g., the nuclear membrane, then launching a cascade of processes upon binding to this site. In other situations, the target can be located on the outer boundary ( figure 1(b)). It can be a tiny aperture-an escape window, in which case the shell-like domain delimited by two boundaries can be regarded as the cortex region, while the inner domain may represent, e.g., the centrosome, as studied in [6,7]. Such a geometrical setup differs from (and is more complicated than) usually studied geometrical settings of the by-now classical Narrow Escape Problem (NEP) [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] due to the presence of a centrosome. From a different perspective, this case can be viewed as an initial step in cell-to-cell communication processes [2][3][4][5] (see also recent results in [24,25]). In this important situation, a particle is the first signalling molecule emitted at a specified location within the intracellular medium, which then has to engage with the target on the cellular membrane. After a reaction with this region, the cell secretes a ligand that moves diffusively in the extracellular domain until it binds to another cell. The binding event is followed by the so-called internalisation process: the signal propagates within the second cell in a cascade of successive reactions, as described above.
Understanding the kinetic behaviour of such multi-stage processes comprising a specific reaction event as one of its controlling factors is impossible without the knowledge of how long such a single reaction stage lasts, starting with the launch of a diffusive particle and terminating at the instant of a successful reaction event. The duration of this stage is a random variable, which in what follows we call the first reaction time (FRT). Its distribution can be rather broad even in bounded systems with a simple geometry [26,27]. From the mathematical point of view, the derivation of the probability density function (PDF) of FRTs is well-established [28][29][30]: it consists of solving the diffusion equation within the domain under consideration, subject to the appropriate boundary conditions. This solution determines the so-called survival probability, i.e., the probability that the particle did not react up to the time instant t, and the desired FPT PDF is deduced as the derivative of the survival probability with respect to t, taken with the minus sign. However, the domain may have a complicated shape, and even more crucially, the appropriate boundary conditions for chemical reactions are the so-called mixed boundary conditions: a zero-current boundary condition on the impermeable surfaces and a reactive boundary condition on the target. Thus, in general, this problem has no explicit exact solutions, except for asymptotical results obtained for the NEP in the limit of a vanishingly small size of the escape window [8][9][10] as well as several spectral exact solutions derived for simple geometries [31][32][33]. In the general case, one therefore either has to resort to stochastic simulations (ranging from lattice random walks or basic Monte Carlo schemes to more advanced techniques such as enhanced Green's Function Reaction Dynamics (eGFRD) [34]), or to a numerical analysis of the boundary value problem via standard discretisation schemes (such as, e.g., finite difference or finite element methods). Here, we follow another direction which consists in developing approximate analytical methods. The latter are, of course, more advantageous because they show how the pertinent properties depend on the system parameters, a dependence that can be verified by comparison with simulations results or experimental data. The predictions of approximate methods can be also tested against available exact solutions in simple geometries.
One such approximate method is the self-consistent approximation (SCA) originally developed in [35] to calculate the Smoluchowski-type constant for reactions with a small centre situated on the otherwise impenetrable surface of a spherical domain. In essence, one replaces the mixed Robin-Neumann boundary condition by an inhomogeneous Neumann boundary condition, and then establishes an appropriate closure relation for the current through the surface of the target. This approach has been subsequently invoked to study several reaction-diffusion problems. Specifically, it was used to calculate (i) the mean velocity of a directional motion of a colloid decorated with a catalytic patch, which prompts reactions in the embedding medium and creates a self-propulsive force [36]; (ii) the mean FRT in the NEP with long-ranged interactions with the boundary and an entropic barrier at the aperture [23]; and (iii) the mean FRT for particle binding to a specific site on a stretched DNA [37]. For the latter case the SCA was also compared against the predictions of another approximate approach, the so-called boundary homogenisation method, and was shown to be more accurate than the latter [38]. Extending the SCA beyond the mean rates, the full PDF of the FRT was determined for the NEP in circular and spherical domains [39] and for the binding kinetics to a specific site on an impermeable cylinder, which mimics an elongated DNA, in a bigger cylindrical domain [40]. Moreover, upon a comparison with the numerical solution of the mixed boundary-value problem it was demonstrated that the SCA is a very reliable approximation, whose predictions agree very well with the numerical results obtained through the finite elements method.
Here we study the statistics of the FRTs in the restricted "onion-shell" geometry depicted in figure 1 which consists of two nested bounded domains: an inner domain Ω 1 placed inside a larger domain Ω 2 . A diffusive point-like particle is launched at time t = 0 from an arbitrary fixed position within the shell-like domain Ω = Ω 2 \Ω 1 delimited by the impermeable boundaries ∂Ω 1 and ∂Ω 2 . The particle then searches for the immobile target T located at either of the delimiting boundaries. Such settings correspond to many realistic situations encountered in molecular signal transduction or in the NEP with additional steric constraints [6,7]. Considering the reaction between the particle and the target we pursue the general case in which a reaction (or binding event) is not perfect and takes place only with a finite probability. This defines the intrinsic chemical reactivity κ with 0 < κ < ∞ (see, e.g., [23,29,35,36,[41][42][43][44][45][46]). When κ = ∞ we are in the case of perfect reactions, occurring immediately upon first contact. The FRT in this case is then exactly the first-passage time to the target, similar to Smoluchowski's original assumption [47]. For finite κ, the reaction is not instantaneous and may not complete upon the first encounter of the diffusive particle with the target, thus necessitating repeated diffusive loops and reaction attempts. Clearly, the FRT is always longer than the first-passage time and strongly depends on the value of κ. Our aim here is to calculate the PDF of the FRT.
We proceed as follows. We start with the general situation in which the nested domains have arbitrary shape with (sufficiently smooth) impermeable boundaries. Capitalising on the recent analysis [31], we present a formally exact spectral solution of the problem and then develop an SCA for domains of arbitrary shapes. Note that the SCA has only been worked out previously for some particular geometries. Here, we establish a general theoretical framework which includes previous geometrical settings as particular cases. Moreover, all steps involved in this general approach are clearly identified and will thus be useful and instructive for the analysis of the FRT statistics in other systems. We then apply the developed framework to the case when the domains Ω 1 and Ω 2 are concentric balls, such that Ω has the form of a spherical shell. For this particular case, we present explicit forms of the FRT PDF, discuss its detailed behaviour and also compare it against the formally exact spectral solution, in which the entering matrices are inverted numerically (see below). We note that the obtained FRT PDF can be considered as a "building block" in more complex signal transduction pathways taking place in nested bounded domains (see, e.g. [48] for more details).
The paper is organised as follows. In section 2, we present the mathematical formulation of the FRT problem and its formal spectral solution as obtained in [31]. Mainly we derive a general form of the SCA, which is valid for domains of arbitrary shape and connectivity, even including unbounded domains with a bounded boundary. ‡ In section 3 the developed framework is applied to spherical shell domains, for which we evaluate the novel exact spectral solution and also provide an explicit prediction for the generating function of the FRT. The result is also discussed for some particular limiting ‡ The mathematical term "bounded boundary" is used in the sense of a compact (loosely, "finite") boundary.
cases. Finally, we analyse the corresponding FRT PDF via numerical inversion of the Laplace transform. In section 4 we present a brief discussion of the general form of the obtained FRT PDF and its asymptotic behaviour. In addition, we show how the shape of the PDF depends on the local curvature of the boundary in the vicinity of the target and also on the radius of the inner domain. We conclude in section 5 with a brief summary of our results and outline some perspectives for future research. Details of intermediate calculations and some of the results are relegated to Appendices. In Appendix A we discuss the accuracy of the SCA in spherical-shell domains as compared to the exact spectral solution, in Appendix B we determine the mean FRT, while in Appendix C we sketch an extension of our results to a planar circular annulus domain. Appendix D yields complementary insights via the analysis of the distribution of reaction event locations.

Spectral solution and the SCA in Laplace domain
Consider a point-like particle with diffusivity D which starts at time zero at position x and diffuses within a bounded, d-dimensional Euclidean domain Ω ⊂ R d . The boundary ∂Ω of Ω is assumed to be sufficiently smooth and is reflecting everywhere, except for the target denoted by T . Upon hitting the target, a particle reacts with it with a finite probability which defines the intrinsic chemical reactivity κ [23,35,36]. The FRT τ is a random variable, that is distributed according to the PDF Here S(t|x) = P x {τ > t} is the survival probability, that satisfies the diffusion equation where ∆ is the Laplace operator, subject to the initial condition and the mixed Robin-Neumann boundary conditions Here ∂/∂n is the normal derivative oriented outwards from the domain. The first relation in (4) states that the diffusive flux toward the target T is equal to the reaction flux on that target. In turn, the second relation indicates zero diffusive flux on the remaining part of the impermeable boundary. These boundary conditions can be compactly written as where I T (x) is the indicator function of the subset T that is I T (x) = 1 if x ∈ T and 0 otherwise. The parameter q = κ/D (in units of inverse length) ranges from 0 to infinity and quantifies the interplay between bulk diffusive transport (characterised by D) and surface reaction (characterised by κ). After multiplication by an appropriate length scale R of the confining domain, one can distinguish reaction-limited (small qR) and diffusion-limited (large qR) processes [49][50][51]. The inverse of q can also be interpreted as the size of a typical region around the first arrival point on the boundary, in which the reaction occurs (see [52][53][54] and Appendix D for details).

Spectral solution in Laplace domain
We focus on the Laplace-transform of the first-reaction time PDF, which is related to the Laplace-transformed survival probability viaH(p|x) = 1 −pS(p|x) and thus satisfies the modified Helmholtz equation subject to the mixed boundary conditions As discussed in [31] a formal spectral solution of this problem can be obtained by use of the Dirichlet-to-Neumann operator M p , which associates to a given function f (x) on the boundary ∂Ω another function g(x) on that boundary as follows: In other words, for a given function f (x), one solves the modified Helmholtz equation with Dirichlet boundary condition u(x) = f (x) and then evaluates the normal derivative of u(x). It is known that M p is a pseudo-differential self-adjoint operator which represents the action of the normal derivative onto a solution of the modified Helmholtz equation [55][56][57][58][59]. As a consequence the mixed boundary condition (8) can be written as Keeping the same notation I T for the operator of multiplication by the function I T (x), one can formally invert this operator equation to get from which the solutionH(p|x) can be extended to the whole domain Ω by using the Dirichlet Green's function (see below). Since the boundary ∂Ω was assumed to be bounded, the spectrum of the self-adjoint operator M p with p ≥ 0 is discrete, i.e., there is an infinite sequence of eigenvalues, µ n (x)} forming an orthonormal basis of the space L 2 (∂Ω) of square integrable functions on ∂Ω. Using these eigenfunctions, one can represent the above solution asH where M (p) is the diagonal matrix of eigenvalues µ n (δ n,n ′ denoting the Kronecker symbol), and is the matrix representation of the operator I T with respect to the eigenbasis {v (p) n }. In the above the asterisk denotes the complex conjugate. This construction is an exact solution of the problem defined by (7) and (8), which requires, however, the inversion of the infinite-dimensional matrix M (p) /q + K (p) . We stress that it is valid for any Euclidean domain Ω with a sufficiently smooth bounded boundary and target T of any shape, not necessarily small (it may, in fact, cover the entire boundary ∂Ω), nor necessarily simply-connected; in fact, the solution holds for multiple target, as well.
While the solution (12) was derived for x ∈ ∂Ω, its extension to any x ∈ Ω can be obtained in a standard way by solving the Dirichlet boundary value problem for the modified Helmholtz equation (7), whereG ∞ (x ′ , p|x) is the Dirichlet Green's function for equation (7). Here, the expression in parentheses is the Laplace transform of the probability flux density j ∞ (s, t|x), that yields a probabilistic interpretation to this extension: a diffusing particle first hits the boundary ∂Ω at some point s (the first-passage problem described by j ∞ (s, t|x)), from which it continues searching the target (the first-reaction time problem described by H(t|s)). Their convolution in time domain takes the form of a product in (15) in Laplace domain. Substituting the spectral decomposition (12), one finally gets, with In this way, one extends the eigenfunctions v (p) n (s) of the Dirichlet-to-Neumann operator, defined on the boundary ∂Ω, into the whole domain Ω. Such extensions can also be understood as the eigenfunctions of the associated Steklov problem.
At a first glance, the spectral representation (12) may look useless, as it expresses an unknown but intuitively appealing quantityH(p|x) in terms of several unknown and less clear quantities (eigenfunctions v (p) n , matrices M (p) and K (p) ). In section 3 we will discuss that in some geometric settings the latter quantities can be found explicitly thus rendering the above formal solution suitable for both numerical computations and analytical studies.

Self-consistent approximation
We now turn to the SCA as developed and applied in [23,[35][36][37][38][39][40]. This approximate method consists in replacing the mixed Robin-Neumann boundary condition (5) by an effective Neumann boundary condition. The latter preserves a zero flux boundary condition at the reflecting part ∂Ω\T of the boundary and stipulates that the current through the target T is a constant that does not depend on the spatial coordinates within the target. In other words, one aims at solving the modified problem where the subscript "app" highlights that the solution of this problem,S app (p|x), is meant to approximateS(p|x). The adjustable parameter J has to be determined from the self-consistent closure condition, which requires that the original Robin boundary condition on the target T is satisfied on average: Using the modified boundary condition (19), one gets where |T | is the area of the target T . In turn, setting the problem defined by equations (18) to (19) reads while the closure relation (21) becomes The modified boundary value problem, equations (22) to (23), is generally much simpler than the original one. In particular, it can be solved in an explicit exact form for some geometric settings, see the examples in [23,[35][36][37][38][39][40] and in section 3. In general, one can express its solution by using the eigenbasis of the Dirichlet-to-Neumann operator. In fact, as M p represents the normal derivative in (23), its inversion immediately yields where the C n are given by (13). Substituting this solution into (24), we get We emphasise that the solution (25) with (26) does not require any matrix inversion and is thus explicit, once the eigenbasis of the Dirichlet-to-Neumann operator is known (e.g., in the cases of an interval and of rotation-invariant domains, see the examples in [60] and in section 3). Note also that the expressions for the SCA derived in [23,[35][36][37][38][39][40] could be directly found from this general solution for each considered geometric setting. An extension of the expression (25) to any x ∈ Ω is simplỹ with V (p) n (x) given by (17). The comparison of the spectral solution (12) with (25) indicates that the SCA can formally be understood as a sort of approximate inversion of the matrix M (p) /q + K (p) . More precisely, the SCA is retrieved if (a) (b) Figure 2. Sketch of the geometrical setup in section 3: a spherical shell domain Ω is delimited by the impermeable boundaries of two concentric balls, Ω 1 with radius R 1 and Ω 2 with radius R 2 (R 1 < R 2 ). The small red ball indicates the starting position of the diffusing particle, with spherical coordinates (r, θ). Note that the azimuthal angle φ does not matter in such a geometry due to the axial symmetry. This specific form of approximate inversion is quite unexpected. The first factor is the result of the inversion of the diagonal matrix M (p) /q, if the matrix K (p) was neglected.
In turn, the second factor, which does not depend on n, is a common correction to all diagonal elements due to the matrix K (p) . It is worthwhile noting that this approximation differs from the so-called diagonal approximation developed in [19][20][21][22] in the context of surface-mediated intermittent diffusion, in which the "correction" matrix (like K (p) here) was approximated by keeping only its diagonal elements, allowing to obtain explicit approximate solutions. We also stress that the limit q → ∞ of a perfectly reactive target is rather challenging for the exact spectral solution in (12) because the factor 1/q stands in front of the matrix M (p) representing the dominant contribution of an unbounded operator M p , as compared to the bounded operator I T represented by the matrix K. Concurrently, this limit is trivial within the SCA.
In the next section we apply this general approach to spherical shell domains, whose rotational symmetry greatly simplifies the above solutions and permits to obtain explicit expressions.

Spherical shell domains
We consider two variants of the FRT problem in a spherical shell domain Ω = {x ∈ R 3 : R 1 < |x| < R 2 }, between two concentric spheres of radii R 1 < R 2 , respectively (see figure 2). The boundary ∂Ω of the domain is fully reflecting except for a circular cap region at the North pole, the target, defined in spherical coordinates (r, θ, φ) by the inequality θ < ε. This target can be located either on the boundary of the inner sphere (Problem I) or on the boundary of the outer sphere (Problem II). In our further analysis we concentrate on Problem I with the target on the boundary of the inner sphere. The analysis for the Problem II is very similar and we will merely present the final results without detailed derivation.
Since the boundary ∂Ω of the spherical shell domain consists of the boundaries of two disjoint spheres, ∂Ω = ∂Ω 1 ∪ ∂Ω 2 , one of which is fully reflecting it is convenient to modify the definition (9) of the Dirichlet-to-Neumann operator by restricting its action only on the sphere which contains the target [60,61]. For instance, for Problem I, one defines the operator M p that associates to a given function f (x) on ∂Ω 1 another function g(x) on the same boundary, . (29) 3.1. Problem I. Target on the inner sphere We start with the exact spectral solution discussed in the previous section. Since the Dirichlet-to-Neumann operator does not depend on the target location, we take advantage of the rotational invariance of the shell domain to determine the eigenbasis of M p along with the matrices M (p) and K (p) [31,46,60]. Moreover, as the circular target preserves the axial symmetry of the problem, the solution does not depend on the angle φ. One can therefore keep only the eigenfunctions that are independent of φ, see [60] for details, where P n (z) are Legendre polynomials, and with Here α = p/D, i n (z) = π/(2z)I n+1/2 (z) and k n (z) = 2/(πz)K n+1/2 (z) are the modified spherical Bessel functions of the first and second kind. The prime here and henceforth denotes the derivative with respect to the argument. The radial functions g n (r) satisfy the second-order differential equation with g n (R 1 ) = 1 and g ′ n (R 2 ) = 0. Note that here the eigenfunctions v n (θ) are independent of p. Using the explicit form of the Dirichlet Green's function from [60], one also gets V (p) n (x) = g n (r)v n (θ).
One can also compute explicitly the matrix elements K n,n ′ (see equation (D12) of [31]), which are independent of the Laplace parameter p, where Note that a more general case of multiple non-overlapping circular targets was also considered in [31]. Finally, the integral in (13) can be easily performed to yield Substituting these expressions into (16) we get with the standard convention P −1 (z) ≡ 1. Together with the explicit formulas for the matrices M (p) and K, this is the exact solution of the original problem in the Laplace domain. However, we are unable to invert the infinite-dimensional matrix M (p) /q + K analytically and have to resort to a numerical inversion. In doing so we truncate the matrices M (p) and K at some order and exercise care afterwards that finite-size effects do not matter, such that our numerical solution is valid with any prescribed accuracy.
While our main focus is on the case of a fixed starting point x, it is instructive to consider two other common situations, in which the starting point is uniformly distributed, either in the bulk, or on a sphere of some radius r. The former case can be relevant when the particle is produced inside the domain in a random location. In turn, the second case accounts for situations when a cell has many membrane channels for a given molecular species, and thus the release point to the shell of interest on the inner membrane surface can be viewed as randomly located. In both cases the above solution in (37) is simplified. The surface-averaged solution reads where the orthogonality of the Legendre polynomials leads to the cancellation of all terms with n > 0. The volume-averaged solution involves an additional integral over r, First-reaction times with target regions on boundaries of shell-like domains 13 Multiplying (33) by r 2 and integrating over r from R 1 to R 2 , one gets from which We next turn to the SCA developed in the previous section. In the geometry considered here our equation (27) becomes where the parameter J is determined from (26) as and we took into account the fact that the area |T | of the spherical cap is given explicitly by |T | = 2πR 2 1 (1 − cos ε). The surface-averaged and volume-averaged approximations are particularly simple,H andH In Appendix A, we illustrate the remarkable agreement between the prediction (41) of the SCA and the exact solution (37), even in the case when the target covers half of the inner sphere. Below we consider two particular cases in which the system reduces to previously studied models.
The entire boundary of the inner sphere is a target. In the particular case when the target extends over the whole boundary of the inner sphere, i.e., when ε = π, the matrix K is the identity matrix due to the orthonormality of the eigenfunctions {v (p) n }. Moreover in this case all terms in (37) with n > 0 vanish, such that the exact spectral solution becomesH where g 0 (r) and µ (p) 0 in (31) and (32) are explicitly given by We thus retrieve the textbook solution for a spherical target surrounded by a larger reflecting sphere (see, e.g., [62] for a more detailed discussion). Interestingly, the SCA result (41) also yields the exact result (37) when ε → π. As a consequence, the SCA appears to be accurate not only in the limit of a vanishingly small target, ε → 0, but also exact in the opposite limit ε → π. This is a direct consequence of the self-consistent closure condition and one of the reasons why the SCA yields accurate results even for the intermediate case of large targets.
Limit R 2 → ∞. In the limit R 2 → ∞ the outer reflecting boundary extends to infinity and one retrieves a common setting of a single circular target on a sphere, explored by a particle diffusing in the unbounded space Ω = {x ∈ R 3 : |x| > R 1 }. This is precisely the classical problem of binding of a ligand that diffuses in an extracellular medium to a finite-sized receptor on an impermeable boundary (see, e.g., [1,35]). Recall, however, that former studies of this problem concentrated exclusively on the mean FRTs (or mean reaction rates). Our analysis below shows that the full PDF of the FRTs can be evaluated for such a geometrical setting. Even though the domain itself is unbounded, its boundary ∂Ω is bounded, and the above solution is still applicable. From the asymptotic behaviour of the modified spherical Bessel functions i n (z) and k n (z) one finds that the radial functions g n (r) converge to g n (r) → − k n (αr) k n (αR 1 ) .
These determine the eigenvalues µ (p) n via equation (31). The other quantities remain unchanged. To our knowledge, this is the first time when exact and approximate expressions for the Laplace-transformed PDFH(p|x) are presented in this geometric setting.

Problem II. Target on the boundary of the outer sphere
As we already remarked, the case of the target on the outer sphere is fairly similar to the previously considered Problem I. Therefore, here we just list the modifications by using the known form of the Dirichlet-to-Neumann eigenbasis [60], and with satisfying g n (R 2 ) = 1 and g ′ n (R 1 ) = 0.S In turn, the matrix K remains unchanged, as well as the formulae (37), (41), and (42) forH(p|x),H app (p|x), and J.
In the limit R 1 → 0, the reflecting inner sphere shrinks to a point, and hence, does not hinder anymore the particle dynamics. One thus retrieves the NEP for an escape from a sphere through a circular hole on its boundary. In this limit, radial functions g n (r) converge to and we therefore recover the SCA result derived in [39]  The remarkable accuracy of the SCA prediction is illustrated in Appendix A.

The PDF in time domain
The solution H(t|x) in time domain can be obtained via inverse Laplace transformation by using the residue theorem. As the spherical shell domain is bounded, the Laplace operator with mixed boundary conditions governing the diffusion-reaction dynamics in equations (2) to (4) has a discrete spectrum, i.e.,H(p|x) has infinitely many poles {p k } lying on the negative axis in the complex plane p ∈ C. Formally, these poles are determined as zeros of the determinant of the matrix M (p) /q + K (p) when the matrix is not invertible, yielding a singular, pole-like behaviour ofH(p|x). If {γ where the c n (x) are determined by these residues. In practice it is more convenient to search the poles {p ′ k } ofH app (p|x), that can be used as approximations of {p k }. The approximate poles are determined by an explicit equation on p, at which the parameter J from (42) diverges, The computation of the residue at each pole p ′ k and their properties are discussed for the case of a sphere in [39]. Instead of this analysis we will invert the Laplace transform numerically by using the Talbot algorithm.
S Note that there is a misprint in equation (B2) of [60], in which g ′ n (R) should be replaced by g ′ n (L).   Figure 3 illustrates the behaviour of the PDF H(t|x) and its self-consistent approximation H app (t|x) for a small target located on the inner sphere. Expectedly, the PDF is broader when the target is smaller and less reactive. In particular, for the case q = 1 and ε = 0.1, the PDF spans over 8 orders of magnitude in time. Figure 4 presents the same quantities for the target located at the outer sphere. While we kept the same angular sizes of the target, their linear sizes are R 2 /R 1 = 10 times larger here, which explains why the distributions are narrower in figure 4 in comparison to figure 3. We discuss these PDFs in the next section.

Discussion
In this section we discuss several features of the calculated PDF H(t|x) such as the different regimes it exhibits as well as the asymptotic behaviour. We also evoke the question of the "screening" effects of the boundary on, e.g., a perfectly-reactive target (κ = ∞) and finally analyse the "concealment" effect of the inner sphere.

The structure of the PDF in time domain
One may notice that the PDF H(t|x) evaluated in the previous section as well as the exact spectral solution exhibits four different temporal regimes delimited by three  characteristic time scales. Namely, an initial regime that extends from time zero up to the most probable (typical) FRT t mp ; a second regime in which the PDF descends from the maximal value; then a plateau-like regime and, ultimately, an exponential decay at long times. We note that this is quite a generic behaviour of the PDFs of the FRTs in bounded domains, which has been amply discussed in previous papers, see, e.g., [39,40,62]. Here, for the sake of completeness, we present a succinct account of the different temporal regimes and the corresponding asymptotic behaviour. At short times H(t|x) is characterised by the singular Lévy-Smirnov-type form [28] where A t is a computable amplitude and ρ is the shortest distance between the starting point and the target. This regime is thus dominated by the "direct" trajectories [62][63][64], which go straight from the initial point to the target. As a consequence, in this regime the dimensionality of the embedding space and of the boundaries does not come into play, unless a particle starts close to the South pole while the target is located at the North pole, such that the optimal path is "interacting" with either of the boundaries. Note that A t in (55) is t-independent for perfect target (κ = ∞), and is an algebraic function of time, A t ∼ κt, for imperfect reactions. In the latter case the particle may be reflected from the target and react with it only upon some subsequent arrival to its location. Therefore, in this latter case H(t|x) is smaller for short times than in the case of perfect reactions. The second regime corresponds to a power-law descent from the maximum. This regime is also universal, in the sense that it is independent of the actual dimensionality of space, and it terminates at the characteristic time t c when a particle first engages with any of the boundaries; thus realising that it moves in a bounded domain. Note that this intermediate power-law decay is responsible for an effective broadness of the PDF. As demonstrated earlier in [26,27], in situations when this regime lasts sufficiently long, the FRTs observed for two realisations of the process may become disproportionally different.
Next a plateau-like regime with a very slow variation of H(t|x) emerges due to the gap between the first and second eigenvalues λ 1 and λ 2 (which are related to the poles, i.e., the solutions of equation (54) above). In this regime all values of the FRT are nearly equally probable.
Finally, at times t ∼ 1/(Dλ 1 ) the PDF crosses over to an exponential decay of the form In summary, there is an appreciable amount of the trajectories that lead to much earlier reaction times, see the discussions in [39,40,62,63,65]. In the case of very small target, the mean FRT is close to 1/(Dλ 1 ) (see Appendix B), and thus the mean FRT controls the long-time behaviour of the PDF [66]. However, the mean FRT, which is typically orders of magnitude longer than the typical reaction time t mp , is therefore insufficient to describe the rich structure of the PDF, especially in the limit of low concentrations typical for many biochemical situations [67]. Strong defocusing of the FRTs is a generic feature of diffusion-controlled reactions, which was commonly ignored in former studies.

Screening effects of the boundary
As we have already remarked, the main mathematical difficulty of the problem at hand is represented by the mixed boundary conditions: a zero-flux boundary condition imposed on the reflecting parts of the boundary and a reactive condition at the target. From a physical point of view, this also implies that the reflecting part of the surface has its effect on the reaction and may, to some unknown extent, "screen" the target. To get some (at least partial) understanding of an emerging effective screening, we compare here the PDFs calculated for two cases: the particle starts at the North Pole of the outer sphere and searches for the target region at the North Pole of the inner sphere (Problem I), and the particle starts at the North pole of the inner sphere and searches for the target region at the North Pole of the outer sphere (Problem II). The initial distance here is the same by definition, the domain in which a particle moves is the same, but in one case the surface in the vicinity of the target is concave, while in the other it is convex. There is, however, a small subtlety concerning the size of the target: to render the two cases identical, should one take the same area of the target or the same angular size? We consider both situations. In figure 5 we compare the FRT PDFs H(t|x) for situations (i) and (ii) in the case of a high intrinsic reactivity, q = 100, which permits us to disentangle the effects of a finite reaction probability from the screening effects of the geometry. Clearly, for such a high reactivity, a particle most likely reacts with the target upon first encounter such that H(t|x) should be very close to the first-passage time PDF. We observe that, indeed, a curvature in the vicinity of the target does matter, regardless of whether we fix the area of the target or its angular size. When the area is fixed ( figure 5(a)), the PDF is slightly broader for the outer-to-inner case (Problem I) than for the inner-to-outer case (Problem II), and also the height of the maximum is lower. As a consequence, transmitting signals from the membrane to the nucleus, under similar biophysical parameters, should usually take longer than in the opposite direction. This seems highly relevant to cell signalling in general, and possibly provides another rationale for the presence of concentric cytoskeletal tracks in the cells. In the case when ε is identical ( figure 5(b)), the effect is more pronounced and the difference is more appealing: the distribution is much broader and the maximum is an order of maximum lower in the outer-to-inner case than in the inner-to-outer case. Therefore, the FRTs in the case when the target is concave are more focused around the most probable value. We note parenthetically that the effect of a local curvature on the first-passage times was already studied within the context of the NEP. In particular, the mean first-passage time was discussed in [68][69][70] for the case when the escape window is located at a corner, or at a cusp in the boundary and on a Riemannian manifold. It was shown there that the very functional dependence of the mean first-passage time on the angular size ε of the aperture can be different, as compared to the case when the boundary is smooth. Our results in figure 5 demonstrate that even in the case of a smooth boundary the effect of the local curvature is present.

Concealment effect of the inner sphere
We turn here to the situation in which the target is located at the North pole of the outer sphere, whereas the starting point is on the South pole of the outer sphere. Evidently, this model can be viewed as the NEP in which the particle dynamics is concealed by the impermeable inner sphere (an obstacle). However, the concealment effect is not evident a priori : On the one hand, one may expect that the presence of such an obstacle should slow down the search process, because the particle now cannot move along "direct trajectories" from the South to the North pole of the outer sphere and has to bypass the obstacle, which makes its trajectories somewhat longer. On the other hand, when the radius of the inner sphere becomes comparable to that of the outer one (i.e., R 1 is close to R 2 ), the particle diffuses in a thin spherical shell and thus undergoes effectively two-dimensional diffusion that may speed up the search process. Which argument is valid here? For a different geometric setting [71], it was argued that obstacles cannot speed up the search process. We thus pose the question of the effect that the radius R 1 of the inner sphere has on the shape of the PDF. In figure 6 we depict H(t|x) for the particular case when the starting point of the particle is located at the South pole of the outer sphere, the radius of the outer sphere is fixed, and the radius of the inner sphere is varied. We find that, in general, the PDF gets more focused around the most probable FRT upon an increase of the radius of the inner sphere, the peak value of the distribution at this point increases, and the distribution becomes more narrow, which implies that the variability of the FRT reduces in the presence of an obstacle.

Conclusions and perspectives
Experimental progress in monitoring biochemical reactions has been massive in the last two decades. Thus it is now possible to monitor gene expression events such as realtime protein production "one protein molecule at a time" [72], providing unprecedented insight into, e.g., bursty protein production [72], specific binding of transcription factor proteins to their binding site on the DNA [73], or elucidating the molecular origins of bacterial individuality [74]. Similarly, single-molecule biochemical reactions in the membranes of living cells can be monitored and rationalised [75]. The control of gene expression is often running off at extremely small concentrations of the involved molecules. Thus, the well-studied Lac and phage lambda repressor proteins occur at copy numbers of only few to few tens in a single, relatively small bacteria cell, corresponding to nanomolar concentrations. It was found that even in such bacteria cells the distance between a gene encoding for a protein A and another gene, to which A is supposed to bind to controls its expression, is relevant [76,77]. Therefore it was concluded that a single reaction rate, as used since Smoluchowski's seminal work [47], cannot account for the complexity of the involved diffusive search [78] and biochemical reactions themselves [79]. Concurrently the role of heterogeneity in cells is being recognised as a key element in understanding intracellular transport [63,[80][81][82][83].
Here we studied the full statistics of first-reaction times to an imperfect immobile target in a typical setting for intracellular reactions [82]. In our model a particle starts from a fixed location and diffuses in a bounded region between two nested impermeable domains, and the target is assumed to be placed on either boundary, a typical situation for many molecular regulation processes naturally running off in biological cells. As in previous work [39,40,48,62,63,84] we obtained clear evidence that the associated reaction times have a broad distribution and thus the mean reaction time (or its inverse, the reaction rate) are not representative for individual reaction events, especially in the low concentration limit.
For the case in which the nested domains have an arbitrary shape with sufficiently smooth impermeable boundaries we presented a formally exact spectral solution of the problem and on this basis developed a general, self-consistent approximation, which has previously been used for several particular geometries. We thus established a general theoretical framework which includes previous geometrical settings as particular cases and in which all steps involved are clearly identified. This framework will be useful and instructive for the analysis of FRT statistics in other systems.
We then demonstrated how to apply the self-consistent approximation to the case when the domains are concentric balls, such that the inter-domain region, in which the particle diffuses, has the form of a spherical shell. For such a geometry we presented explicit forms of the Laplace-transformed FRT PDF evaluated within the self-consistent approximation and showed that it agrees exceptionally well with the numerically-evaluated exact spectral solution.
As we mentioned in the Introduction, our model is quite realistic and appears indeed in many systems of biophysical and biochemical interest. In particular, it corresponds to intermediate steps in either intracellular reaction pathways or can be seen as the initial stage in cell-to-cell communication processes [5], see also the recent work [48]. In view of the relevance to these important fields, several generalisations of the analysis presented here may be important: (i) In many applications, a target may not be unique, as we supposed here, but there rather exists an array of such targets, and the particle of interest may react with either of them in order to trigger the same desired response. For instance, receptors on the cellular membrane may be present in sufficiently big amounts (see, e.g., [5]) and binding to either of them will result in the launch of the second messenger. The general form of the self-consistent approximation derived here allows for an immediate extension to this setting. Actually, the only change consists in replacing the matrix K (which is determined by the shape of the target), and the elements of this matrix were calculated in [33] for the case of multiple non-overlapping targets of circular shape.
(ii) In some instances, the target may themselves perform a slow diffusive motion along the boundary of the domain. Recently, some analytical approaches have been developed to study the survival of a diffusive particle in the presence of diffusive mobile sites in unbounded domains [85][86][87][88][89][90], which ultimately gives an access to the first-reaction time PDF. An extension of these approaches to the case in which a particle diffuses within the inter-boundary region while the target perform diffusive motion along the boundaries represents an interesting problem in its own right.
(iii) At some stages of the intercellular signalling, the signal can be "amplified" [5]. In other words, instead of a single particle, initially some amount N of identical particles are launched. This is an important aspect when an extreme event in which the arrival of the first out of N particles matters. The impact of such an extreme event on the FRTs and their PDF has been rather extensively discussed within the recent years, and several analytical analyses have been proposed (see, e.g., [84,[91][92][93][94] and references therein). An extension of such analyses for the geometrical setting studied here should be of interest.
(iv) Finally we mention that inhomogeneities of the environment should ultimately be included in the description. Such conditions may either be represented by positiondependent diffusivities, by local patches with their diffusion coefficient, or lead to anomalous diffusion statistics. We also mention scenarios of two-channel diffusion of a particle, that switches between inert and reactive "channels" with two different probabilities [65], or more general switching diffusion models that may account for reversible binding to other molecules or organelles [67,95,96]. Yet another extension consists in the incorporation of radial active transport on concentric cytoskeletal tracks (e.g. microtubules), in which the signal molecules are transported by a mixed diffusiveconvective motion (with convective motion in radial directions). All of these points will deserve attention in future work.

Acknowledgments
DG acknowledges the Alexander von Humboldt Foundation for support within a Bessel Prize award. RM acknowledges the German Science Foundation (DFG, grant no. ME 1535/12-1) for support. RM also acknowledges the Foundation for Polish Science (Fundacja na rzecz Nauki Polskiej, FNP) for support within an Alexander von Humboldt Honorary Polish Research Scholarship.

Appendix A. Validation of the self-consistent approximation
In this Appendix, we discuss the accuracy of the self-consistent approximation in Laplace domain.
We start by considering the case when the target is located on the inner sphere. Figure A1 illustrates the remarkable agreement between the prediction (41) of the SCA and the exact solution (37). Note that both expressions were truncated at n max = 25, i.e., the size of the (infinite-dimensional) matrices M (p) and K was set to (n max + 1) × (n max + 1), and also in (41) we truncated the series at n = 25. We verified that this truncation does not affect the accuracy of the exact solution (e.g., an increase of n max from 25 to 50 does not change the solution visibly). The two left panels (a) and (c) showH(p|x) for a small target (ε = 0.1) with moderate (q = 1) and high (q = 100) reactivity, respectively. Expectedly the SCA appears to be very accurate for small targets, regardless of the actual value of q. Most surprisingly, however, the SCA turns out to be accurate even for a large target (ε = π/2, i.e., half of the inner sphere), see the right panels (b) and (d). We note that the starting point is relatively far from the target which is the most "favourable" configuration for the applicability of the SCA. Yet, this is also a typical situation in many applications. For instance, for the signal transduction processes a particle starts on one of the boundaries and has to reach a target on the other boundary, i.e., it appears initially quite far from the target. The relative error of the self-consistent approximation is shown in figure A2 for the same set of parameters. Expectedly, the error is smaller when the target is smaller and less reactive. But even in the worst case of a large, highly reactive target (ε = π/2, q = 100), the relative error remains below 5% for a very broad range of p-values, in spite of the fact thatH(p|x) varies over more than 10 orders of magnitude.
The accuracy of the SCA remains excellent for Problem II when the target is located on the outer boundary. Figure A3 compares the SCA prediction against the exact spectral solution, whereas figure A4 illustrates its relative error. The agreement is again remarkably good, with the relative error raising only up to 8% for large values of p, which corresponds to the short-t behaviour of the PDF.

Appendix D. Spread harmonic measure
In this Appendix, we provide some insights into the role of the parameter q characterising the reactivity of the target. For the sake of brevity, we focus on the case when the target covers the whole inner sphere. After the first arrival onto the inner sphere, the particle does not necessarily react immediately but may be reflected to continue its bulk diffusion, return to the target and be reflected again, and so on, until an eventual reaction after a number of failed attempts. As a consequence, the location of the reaction event is in general different from the location of the first arrival. The distribution of reaction locations is called the spread harmonic measure [53]. Some properties of this measure and the related interpretation of the parameter q were discussed in [52][53][54] for a restricted class of domains. In general, the spread harmonic measure density ω q (s|s 0 ), characterising the probability of the reaction event at a boundary point s after the arrival onto the boundary point s 0 , can be expressed in terms of the eigenfunctions of the Dirichlet-to-Neumann operator [46] ω q (s|s 0 ) = where sin θ comes from using spherical coordinates. The integral of ω q (θ ′ ) from 0 to θ is the probability that the reaction location is on the spherical cap of angle θ, P q (θ) = θ 0 dθ ′ ω q (θ ′ ) = 1 2 ∞ n=0 P n−1 (cos θ) − P n+1 (cos θ) Setting P q (Θ t ) = 1/2 determines implicitly the angular size Θ t of the spherical cap, on which half of the reaction events occur. Figure D1(a) illustrates the behaviour of this probability for three values of q.
(D. 8) In the limit of high reactivity, q → ∞, the spread harmonic measure density converges to a Dirac distribution around the first arrival point. Since c k ∝ k −3 , one gets the asymptotic behaviour Θ m ∝ ln(qR 1 )/(qR 1 ). In other words, for large q, the inverse of q is proportional to the mean geodesic distance R 1 Θ m , up to a logarithmic correction. In the opposite limit of low reactivity, q → 0, the spread harmonic measure density is getting uniform, and one gets Θ m → π/2. This reflects the fact that the inner sphere, as well as all geodesic distance on it, are bounded. In other words, the mean geodesic distance cannot grow as 1/q forever, as it was the case for the plane [52,54]. The dependence of Θ m on q is illustrated on figure D1(b).