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

Denis S Grebenkov1, Ralf Metzler2 & Gleb Oshanin3 1Laboratoire de Physique de
la Matière Condensée (UMR 7643), CNRS – Ecole Polytechnique, IP Paris, 91120
Palaiseau, France  2Institute of Physics and Astronomy, University of Potsdam,
14476 Potsdam-Golm, Germany  3Sorbonne Université, CNRS, Laboratoire de
Physique Théorique de la Matière Condensée (UMR CNRS 7600), 4 Place Jussieu,
75252 Paris Cedex 5, France

###### Abstract

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.

##  1 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 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.

Figure 1: Sketch of our geometrical setup. A domain $\Omega_{2}$ with an
impermeable boundary ${\partial\Omega}_{2}$ is nesting a smaller domain
$\Omega_{1}$ enclosed by an impermeable boundary ${\partial\Omega}_{1}$. A
particle, whose starting position is indicated by the blue point, diffuses
within the shell-like domain $\Omega=\Omega_{2}\backslash\Omega_{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).

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 [9, 10, 8] 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 $\Omega_{1}$ placed inside a larger domain $\Omega_{2}$. A
diffusive point-like particle is launched at time $t=0$ from an _arbitrary_
fixed position within the shell-like domain
$\Omega=\Omega_{2}\backslash\Omega_{1}$ delimited by the impermeable
boundaries ${\partial\Omega}_{1}$ and ${\partial\Omega}_{2}$. The particle
then searches for the immobile target ${\cal 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 $\kappa$ with $0<\kappa<\infty$ (see, e.g.,
[42, 35, 36, 23, 29, 41, 43, 44, 45, 46]). When $\kappa=\infty$ 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 $\kappa$, 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 $\kappa$. 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
$\Omega_{1}$ and $\Omega_{2}$ are concentric balls, such that $\Omega$ 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.111The mathematical term ”bounded boundary” is used in the
sense of a compact (loosely, ”finite”) 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 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 A we discuss the accuracy of the SCA in spherical-shell domains as compared
to the exact spectral solution, in B we determine the mean FRT, while in C we
sketch an extension of our results to a planar circular annulus domain. D
yields complementary insights via the analysis of the distribution of reaction
event locations.

##  2 Spectral solution and the SCA in Laplace domain

Consider a point-like particle with diffusivity $D$ which starts at time zero
at position $\mathbf{x}$ and diffuses within a bounded, $d$-dimensional
Euclidean domain $\Omega\subset\mathbb{R}^{d}$. The boundary
${\partial\Omega}$ of $\Omega$ is assumed to be sufficiently smooth and is
reflecting everywhere, except for the target denoted by $\mathcal{T}$. Upon
hitting the target, a particle reacts with it with a finite probability which
defines the intrinsic chemical reactivity $\kappa$ [35, 36, 23]. The FRT
$\tau$ is a random variable, that is distributed according to the PDF

$H(t|\mathbf{x})=-\frac{\partial S(t|\mathbf{x})}{\partial t}.$ (1)

Here $S(t|\mathbf{x})=\mathbb{P}_{\mathbf{x}}\\{\tau>t\\}$ is the survival
probability, that satisfies the diffusion equation [28]

$\frac{\partial S(t|\mathbf{x})}{\partial t}=D\Delta
S(t|\mathbf{x})\quad(\mathbf{x}\in\Omega),$ (2)

where $\Delta$ is the Laplace operator, subject to the initial condition

$S(0,\mathbf{x})=1$ (3)

and the mixed Robin-Neumann boundary conditions

$\left\\{\begin{array}[]{ll}-D\partial S(t,\mathbf{x})/\partial n=\kappa
S(t,\mathbf{x})&(\mathbf{x}\in\mathcal{T})\\\ -D\partial
S(t,\mathbf{x})/\partial
n=0&(\mathbf{x}\in{\partial\Omega}\backslash\mathcal{T})\end{array}\right..$
(4)

Here $\partial/\partial n$ is the normal derivative oriented outwards from the
domain. The first relation in (4) states that the diffusive flux toward the
target $\mathcal{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

$-\frac{\partial S(t|\mathbf{x})}{\partial
n}=q\mathbb{I}_{\mathcal{T}}(\mathbf{x})S(t|\mathbf{x})\quad(\mathbf{x}\in{\partial\Omega}),$
(5)

where $\mathbb{I}_{\mathcal{T}}(\mathbf{x})$ is the indicator function of the
subset $\mathcal{T}$ that is $\mathbb{I}_{\mathcal{T}}(\mathbf{x})=1$ if
$\mathbf{x}\in\mathcal{T}$ and $0$ otherwise.

The parameter $q=\kappa/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 $\kappa$). 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 D for details).

###  2.1 Spectral solution in Laplace domain

We focus on the Laplace-transform

$\tilde{H}(p|\mathbf{x})=\int\limits_{0}^{\infty}dte^{-pt}H(t|\mathbf{x})$ (6)

of the first-reaction time PDF, which is related to the Laplace-transformed
survival probability via $\tilde{H}(p|\mathbf{x})=1-p\tilde{S}(p|\mathbf{x})$
and thus satisfies the modified Helmholtz equation

$(p-D\Delta)\tilde{H}(p|\mathbf{x})=0\quad(\mathbf{x}\in\Omega)$ (7)

subject to the mixed boundary conditions

$\frac{\partial\tilde{H}(p|\mathbf{x})}{\partial
n}+q\mathbb{I}_{\mathcal{T}}(\mathbf{x})\tilde{H}(p|\mathbf{x})=q\mathbb{I}_{\mathcal{T}}(\mathbf{x})\quad(\mathbf{x}\in{\partial\Omega}).$
(8)

As discussed in [31] a formal spectral solution of this problem can be
obtained by use of the Dirichlet-to-Neumann operator $\mathcal{M}_{p}$, which
associates to a given function $f(\mathbf{x})$ on the boundary
${\partial\Omega}$ another function $g(\mathbf{x})$ on that boundary as
follows:

$\mathcal{M}_{p}~{}:~{}f(\mathbf{x})\to
g(\mathbf{x})=\left.\left(\frac{\partial u(\mathbf{x})}{\partial
n}\right)\right|_{{\partial\Omega}},\mbox{ where
}\left\\{\begin{array}[]{ll}(p-D\Delta)u(\mathbf{x})=0&(\mathbf{x}\in\Omega)\\\
u(\mathbf{x})=f(\mathbf{x})&(\mathbf{x}\in{\partial\Omega})\end{array}\right..$
(9)

In other words, for a given function $f(\mathbf{x})$, one solves the modified
Helmholtz equation with Dirichlet boundary condition
$u(\mathbf{x})=f(\mathbf{x})$ and then evaluates the normal derivative of
$u(\mathbf{x})$. It is known that $\mathcal{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

$\mathcal{M}_{p}\tilde{H}(p|\mathbf{x})+q\mathbb{I}_{\mathcal{T}}(\mathbf{x})\tilde{H}(p|\mathbf{x})=q\mathbb{I}_{\mathcal{T}}(\mathbf{x}).$
(10)

Keeping the same notation $\mathbb{I}_{\mathcal{T}}$ for the operator of
multiplication by the function $\mathbb{I}_{\mathcal{T}}(\mathbf{x})$, one can
formally invert this operator equation to get

$\tilde{H}(p|\mathbf{x})=\bigl{(}\mathcal{M}_{p}/q+\mathbb{I}_{\mathcal{T}}\bigr{)}^{-1}\mathbb{I}_{\mathcal{T}}(\mathbf{x})\qquad(\mathbf{x}\in{\partial\Omega}),$
(11)

from which the solution $\tilde{H}(p|\mathbf{x})$ can be extended to the whole
domain $\Omega$ by using the Dirichlet Green’s function (see below). Since the
boundary ${\partial\Omega}$ was assumed to be bounded, the spectrum of the
self-adjoint operator $\mathcal{M}_{p}$ with $p\geq 0$ is discrete, i.e.,
there is an infinite sequence of eigenvalues,
$\mu_{0}^{(p)}\leq\mu_{1}^{(p)}\leq\ldots\nearrow+\infty$, associated to
eigenfunctions $\\{v_{n}^{(p)}(\mathbf{x})\\}$ forming an orthonormal basis of
the space $L_{2}({\partial\Omega})$ of square integrable functions on
${\partial\Omega}$. Using these eigenfunctions, one can represent the above
solution as

$\tilde{H}(p|\mathbf{x})=\sum\limits_{n,n^{\prime}=0}^{\infty}v_{n}^{(p)}(\mathbf{x})\biggl{[}\bigl{(}\mathbf{M}^{(p)}/q+\mathbf{K}^{(p)}\bigr{)}^{-1}\biggr{]}_{n,n^{\prime}}C_{n^{\prime}}^{*}\qquad(\mathbf{x}\in{\partial\Omega}),$
(12)

where $\mathbf{M}^{(p)}$ is the diagonal matrix of eigenvalues
$\mu_{n}^{(p)}$,
$\mathbf{M}^{(p)}_{n,n^{\prime}}=\delta_{n,n^{\prime}}\mu_{n}^{(p)}$
($\delta_{n,n^{\prime}}$ denoting the Kronecker symbol),

$C_{n}=\int\limits_{{\partial\Omega}}d\mathbf{x}\mathbb{I}_{\mathcal{T}}(\mathbf{x})v_{n}^{(p)}(\mathbf{x})=\int\limits_{\mathcal{T}}d\mathbf{x}v_{n}^{(p)}(\mathbf{x}),$
(13)

and

$\bigl{[}\mathbf{K}^{(p)}\bigr{]}_{n,n^{\prime}}=\int\limits_{{\partial\Omega}}d\mathbf{x}v_{n}^{(p)}(\mathbf{x})\mathbb{I}_{\mathcal{T}}(\mathbf{x})[v_{n^{\prime}}^{(p)}(\mathbf{x})]^{*}=\int\limits_{\mathcal{T}}d\mathbf{x}v_{n}^{(p)}(\mathbf{x})[v_{n^{\prime}}^{(p)}(\mathbf{x})]^{*}$
(14)

is the matrix representation of the operator $\mathbb{I}_{\mathcal{T}}$ with
respect to the eigenbasis $\\{v_{n}^{(p)}\\}$. 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 $\mathbf{M}^{(p)}/q+\mathbf{K}^{(p)}$. We stress
that it is valid for any Euclidean domain $\Omega$ with a sufficiently smooth
bounded boundary and target $\mathcal{T}$ of any shape, not necessarily small
(it may, in fact, cover the entire boundary ${\partial\Omega}$), nor
necessarily simply-connected; in fact, the solution holds for multiple target,
as well.

While the solution (12) was derived for $\mathbf{x}\in{\partial\Omega}$, its
extension to any $\mathbf{x}\in\Omega$ can be obtained in a standard way by
solving the Dirichlet boundary value problem for the modified Helmholtz
equation (7),

$\tilde{H}(p|\mathbf{x})=\int\limits_{{\partial\Omega}}d\mathbf{s}\tilde{H}(p|\mathbf{s})\underbrace{\left(-D\frac{\partial}{\partial
n}\tilde{G}_{\infty}(\mathbf{x}^{\prime},p|\mathbf{x})\right)_{\mathbf{x}^{\prime}=\mathbf{s}}}_{=\tilde{j}_{\infty}(\mathbf{s},p|\mathbf{x})},$
(15)

where $\tilde{G}_{\infty}(\mathbf{x}^{\prime},p|\mathbf{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_{\infty}(\mathbf{s},t|\mathbf{x})$, that yields a probabilistic
interpretation to this extension: a diffusing particle first hits the boundary
${\partial\Omega}$ at some point $\mathbf{s}$ (the first-passage problem
described by $j_{\infty}(\mathbf{s},t|\mathbf{x})$), from which it continues
searching the target (the first-reaction time problem described by
$H(t|\mathbf{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,

$\tilde{H}(p|\mathbf{x})=\sum\limits_{n,n^{\prime}=0}^{\infty}V_{n}^{(p)}(\mathbf{x})\biggl{[}\bigl{(}\mathbf{M}^{(p)}/q+\mathbf{K}^{(p)}\bigr{)}^{-1}\biggr{]}_{n,n^{\prime}}C_{n^{\prime}}^{*}\qquad(\mathbf{x}\in\Omega),$
(16)

with

$V_{n}^{(p)}(\mathbf{x})=\int\limits_{{\partial\Omega}}d\mathbf{s}v_{n}^{(p)}(\mathbf{s})\tilde{j}_{\infty}(\mathbf{s},p|\mathbf{x}).$
(17)

In this way, one extends the eigenfunctions $v_{n}^{(p)}(\mathbf{s})$ of the
Dirichlet-to-Neumann operator, defined on the boundary ${\partial\Omega}$,
into the whole domain $\Omega$. 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 quantity
$\tilde{H}(p|\mathbf{x})$ in terms of several unknown and less clear
quantities (eigenfunctions $v_{n}^{(p)}$, matrices $\mathbf{M}^{(p)}$ and
$\mathbf{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.

###  2.2 Self-consistent approximation

We now turn to the SCA as developed and applied in [35, 36, 23, 37, 38, 40,
39]. 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
${\partial\Omega}\backslash\mathcal{T}$ of the boundary and stipulates that
the current through the target $\mathcal{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

$\displaystyle(p-D\Delta)\tilde{S}_{\rm app}(p|\mathbf{x})$
$\displaystyle=1\quad(\mathbf{x}\in\Omega),$ (18)
$\displaystyle-\frac{\partial\tilde{S}_{\rm app}(p|\mathbf{x})}{\partial n}$
$\displaystyle=\frac{J}{p}\mathbb{I}_{\mathcal{T}}(\mathbf{x})\quad(\mathbf{x}\in{\partial\Omega}),$
(19)

where the subscript "app" highlights that the solution of this problem,
$\tilde{S}_{\rm app}(p|\mathbf{x})$, is meant to approximate
$\tilde{S}(p|\mathbf{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 $\mathcal{T}$ is satisfied on average:

$\int\limits_{\mathcal{T}}d\mathbf{x}\left(-\frac{\partial\tilde{S}_{\rm
app}(p|\mathbf{x})}{\partial
n}\right)=\int\limits_{\mathcal{T}}d\mathbf{x}\bigl{(}q\tilde{S}_{\rm
app}(p|\mathbf{x})\bigr{)}.$ (20)

Using the modified boundary condition (19), one gets

$J=\frac{qp}{|\mathcal{T}|}\int\limits_{\mathcal{T}}d\mathbf{x}\tilde{S}_{\rm
app}(p|\mathbf{x}),$ (21)

where $|\mathcal{T}|$ is the area of the target $\mathcal{T}$. In turn,
setting

$\tilde{H}_{\rm app}(p|\mathbf{x})=1-p\tilde{S}_{\rm
app}(p|\mathbf{x})=J\tilde{h}_{\rm app}(p|\mathbf{x}),$

the problem defined by equations (18) to (19) reads

$\displaystyle(p-D\Delta)\tilde{h}_{\rm app}(p|\mathbf{x})$
$\displaystyle=0\quad(\mathbf{x}\in\Omega),$ (22)
$\displaystyle\frac{\partial\tilde{h}_{\rm app}(p|\mathbf{x})}{\partial n}$
$\displaystyle=\mathbb{I}_{\mathcal{T}}(\mathbf{x})\quad(\mathbf{x}\in{\partial\Omega}),$
(23)

while the closure relation (21) becomes

$J=q\left(1-\frac{J}{|\mathcal{T}|}\int\limits_{\mathcal{T}}d\mathbf{x}\tilde{h}_{\rm
app}(p|\mathbf{x})\right),$

from which

$J=\left(\frac{1}{q}+\frac{1}{|\mathcal{T}|}\int\limits_{\mathcal{T}}d\mathbf{x}\tilde{h}_{\rm
app}(p|\mathbf{x})\right)^{-1}.$ (24)

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 [35, 36, 23, 37,
38, 40, 39] and in section 3. In general, one can express its solution by
using the eigenbasis of the Dirichlet-to-Neumann operator. In fact, as
$\mathcal{M}_{p}$ represents the normal derivative in (23), its inversion
immediately yields

$\tilde{H}_{\rm
app}(p|\mathbf{x})=J\mathcal{M}_{p}^{-1}\mathbb{I}_{\mathcal{T}}(\mathbf{x})=\sum\limits_{n=0}^{\infty}\frac{v_{n}^{(p)}(\mathbf{x})}{\mu_{n}^{(p)}/J}C_{n}^{*}\qquad(\mathbf{x}\in{\partial\Omega}),$
(25)

where the $C_{n}$ are given by (13). Substituting this solution into (24), we
get

$J=\left(\frac{1}{q}+\frac{1}{|\mathcal{T}|}\sum\limits_{n=0}^{\infty}\frac{|C_{n}|^{2}}{\mu_{n}^{(p)}}\right)^{-1}.$
(26)

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 [35, 36, 23, 37, 38, 40, 39] could be
directly found from this general solution for each considered geometric
setting. An extension of the expression (25) to any $\mathbf{x}\in\Omega$ is
simply

$\tilde{H}_{\rm
app}(p|\mathbf{x})=\sum\limits_{n=0}^{\infty}\frac{V_{n}^{(p)}(\mathbf{x})}{\mu_{n}^{(p)}/J}C_{n}^{*},$
(27)

with $V_{n}^{(p)}(\mathbf{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
$\mathbf{M}^{(p)}/q+\mathbf{K}^{(p)}$. More precisely, the SCA is retrieved if

$\biggl{[}\bigl{(}\mathbf{M}^{(p)}/q+\mathbf{K}^{(p)}\bigr{)}^{-1}\biggr{]}_{n,n^{\prime}}\approx\frac{\delta_{n,n^{\prime}}}{\mu_{n}^{(p)}/q}\left(1+\frac{q}{|\mathcal{T}|}\sum\limits_{n=0}^{\infty}\frac{|C_{n}|^{2}}{\mu_{n}^{(p)}}\right)^{-1}.$
(28)

This specific form of approximate inversion is quite unexpected. The first
factor is the result of the inversion of the diagonal matrix
$\mathbf{M}^{(p)}/q$, if the matrix $\mathbf{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 $\mathbf{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
$\mathbf{K}^{(p)}$ here) was approximated by keeping only its diagonal
elements, allowing to obtain explicit approximate solutions. We also stress
that the limit $q\to\infty$ 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 $\mathbf{M}^{(p)}$ representing the dominant
contribution of an unbounded operator $\mathcal{M}_{p}$, as compared to the
bounded operator $\mathbb{I}_{\mathcal{T}}$ represented by the matrix
$\mathbf{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.

##  3 Spherical shell domains

Figure 2: Sketch of the geometrical setup in section 3: a spherical shell
domain $\Omega$ is delimited by the impermeable boundaries of two concentric
balls, $\Omega_{1}$ with radius $R_{1}$ and $\Omega_{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,\theta)$. Note that the
azimuthal angle $\phi$ does not matter in such a geometry due to the axial
symmetry. (a) The target (dark blue region) is located on the inner surface.
(b) The target (light red region) is located on the outer surface.

We consider two variants of the FRT problem in a spherical shell domain
$\Omega=\\{\mathbf{x}\in\mathbb{R}^{3}~{}:~{}R_{1}<|\mathbf{x}|<R_{2}\\}$,
between two concentric spheres of radii $R_{1}<R_{2}$, respectively (see
figure 2). The boundary ${\partial\Omega}$ of the domain is fully reflecting
except for a circular cap region at the North pole, the target, defined in
spherical coordinates $(r,\theta,\phi)$ by the inequality
$\theta<\varepsilon$. 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 ${\partial\Omega}$ of the spherical shell domain consists
of the boundaries of two disjoint spheres,
${\partial\Omega}={\partial\Omega}_{1}\cup{\partial\Omega}_{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 $\mathcal{M}_{p}$ that associates to a given function
$f(\mathbf{x})$ on ${\partial\Omega}_{1}$ another function $g(\mathbf{x})$ on
the same boundary,

$\mathcal{M}_{p}:f(\mathbf{x})\to g(\mathbf{x})=\left.\left(\frac{\partial
u(\mathbf{x})}{\partial n}\right)\right|_{{\partial\Omega}_{1}},\mbox{where
}\left\\{\begin{array}[]{ll}(p-D\Delta)u(\mathbf{x})=0&(\mathbf{x}\in\Omega)\\\
u(\mathbf{x})=f(\mathbf{x})&(\mathbf{x}\in{\partial\Omega}_{1})\\\ \partial
u(\mathbf{x})/\partial
n=0&(\mathbf{x}\in{\partial\Omega}_{2})\end{array}\right..$ (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 $\mathcal{M}_{p}$ along with the matrices
$\mathbf{M}^{(p)}$ and $\mathbf{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 $\phi$. One can therefore keep only the eigenfunctions
that are independent of $\phi$, see [60] for details,

$v_{n}(\theta)=\frac{1}{R_{1}}\sqrt{\frac{2n+1}{4\pi}}P_{n}(\cos\theta),$ (30)

where $P_{n}(z)$ are Legendre polynomials, and

$\mu_{n}^{(p)}=-g^{\prime}_{n}(R_{1}),$ (31)

with

$g_{n}(r)=\frac{k^{\prime}_{n}(\alpha R_{2})i_{n}(\alpha
r)-i^{\prime}_{n}(\alpha R_{2})k_{n}(\alpha r)}{k^{\prime}_{n}(\alpha
R_{2})i_{n}(\alpha R_{1})-i^{\prime}_{n}(\alpha R_{2})k_{n}(\alpha R_{1})}.$
(32)

Here $\alpha=\sqrt{p/D}$, $i_{n}(z)=\sqrt{\pi/(2z)}I_{n+1/2}(z)$ and
$k_{n}(z)=\sqrt{2/(\pi 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

$g^{\prime\prime}_{n}(r)+\frac{2}{r}g^{\prime}_{n}(r)-\frac{n(n+1)}{r^{2}}g_{n}(r)-\alpha^{2}g_{n}(r)=0,$
(33)

with $g_{n}(R_{1})=1$ and $g^{\prime}_{n}(R_{2})=0$. Note that here the
eigenfunctions $v_{n}(\theta)$ are independent of $p$. Using the explicit form
of the Dirichlet Green’s function from [60], one also gets
$V_{n}^{(p)}(\mathbf{x})=g_{n}(r)v_{n}(\theta)$.

One can also compute explicitly the matrix elements
$\mathbf{K}_{n,n^{\prime}}$ (see equation (D12) of [31]), which are
independent of the Laplace parameter $p$,

$\mathbf{K}_{n,n^{\prime}}=\sqrt{(n+1/2)(n^{\prime}+1/2)}\sum\limits_{k=0}^{\min\\{n,n^{\prime}\\}}B_{nn^{\prime}}^{k}\frac{P_{n+n^{\prime}-2k-1}(\cos\varepsilon)-P_{n+n^{\prime}-2k+1}(\cos\varepsilon)}{2(n+n^{\prime}-2k)+1},$
(34)

where

$B_{nn^{\prime}}^{k}=\frac{A_{k}A_{n-k}A_{n^{\prime}-k}}{A_{n+n^{\prime}-k}}\frac{2n+2n^{\prime}-4k+1}{2n+2n^{\prime}-2k+1},\qquad
A_{k}=\frac{\Gamma(k+1/2)}{\sqrt{\pi}\Gamma(k+1)}.$ (35)

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

$C_{n}=R_{1}\frac{\sqrt{\pi}}{\sqrt{2n+1}}\bigl{(}P_{n-1}(\cos\varepsilon)-P_{n+1}(\cos\varepsilon)\bigr{)}.$
(36)

Substituting these expressions into (16) we get

$\tilde{H}(p|\mathbf{x})=\sum\limits_{n,n^{\prime}=0}^{\infty}\sqrt{2n+1}g_{n}(r)P_{n}(\cos\theta)\biggl{[}\bigl{(}\mathbf{M}^{(p)}/q+\mathbf{K}\bigr{)}^{-1}\biggr{]}_{n,n^{\prime}}\frac{P_{n^{\prime}-1}(\cos\varepsilon)-P_{n^{\prime}+1}(\cos\varepsilon)}{2\sqrt{2n^{\prime}+1}},$
(37)

with the standard convention $P_{-1}(z)\equiv 1$. Together with the explicit
formulas for the matrices $\mathbf{M}^{(p)}$ and $\mathbf{K}$, this is the
exact solution of the original problem in the Laplace domain. However, we are
unable to invert the infinite-dimensional matrix
$\mathbf{M}^{(p)}/q+\mathbf{K}$ analytically and have to resort to a numerical
inversion. In doing so we truncate the matrices $\mathbf{M}^{(p)}$ and
$\mathbf{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 $\mathbf{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

$\displaystyle\overline{\tilde{H}(p|r)}$ $\displaystyle=\frac{1}{4\pi
r^{2}}\int\limits_{|\mathbf{x}|=r}d\mathbf{x}\tilde{H}(p|\mathbf{x})$
$\displaystyle=\sum\limits_{n^{\prime}=0}^{\infty}g_{0}(r)\biggl{[}\bigl{(}\mathbf{M}^{(p)}/q+\mathbf{K}\bigr{)}^{-1}\biggr{]}_{0,n^{\prime}}\frac{P_{n^{\prime}-1}(\cos\varepsilon)-P_{n^{\prime}+1}(\cos\varepsilon)}{2\sqrt{2n^{\prime}+1}},$
(38)

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$,

$\displaystyle\overline{\tilde{H}(p)}$
$\displaystyle=\frac{1}{4\pi(R_{2}^{3}-R_{1}^{3})/3}\int\limits_{\Omega}d\mathbf{x}\tilde{H}(p|\mathbf{x})=\frac{3}{R_{2}^{3}-R_{1}^{3}}\int\limits_{R_{1}}^{R_{2}}drr^{2}\overline{\tilde{H}(p|r)}.$

Multiplying (33) by $r^{2}$ and integrating over $r$ from $R_{1}$ to $R_{2}$,
one gets

$\int\limits_{R_{1}}^{R_{2}}drr^{2}g_{0}(r)=\frac{R_{2}^{2}g^{\prime}_{0}(R_{2})-R_{1}^{2}g^{\prime}_{0}(R_{1})}{\alpha^{2}}=\frac{R_{1}^{2}\mu_{0}^{(p)}}{p/D},$
(39)

from which

$\overline{\tilde{H}(p)}=\frac{3DR_{1}^{2}\mu_{0}^{(p)}}{p(R_{2}^{3}-R_{1}^{3})}\sum\limits_{n^{\prime}=0}^{\infty}\biggl{[}\bigl{(}\mathbf{M}^{(p)}/q+\mathbf{K}\bigr{)}^{-1}\biggr{]}_{0,n^{\prime}}\frac{P_{n^{\prime}-1}(\cos\varepsilon)-P_{n^{\prime}+1}(\cos\varepsilon)}{2\sqrt{2n^{\prime}+1}}.$
(40)

We next turn to the SCA developed in the previous section. In the geometry
considered here our equation (27) becomes

$\tilde{H}_{\rm
app}(p|\mathbf{x})=J\sum\limits_{n=0}^{\infty}g_{n}(r)P_{n}(\cos\theta)\frac{1}{\mu_{n}^{(p)}}\frac{P_{n-1}(\cos\varepsilon)-P_{n+1}(\cos\varepsilon)}{2},$
(41)

where the parameter $J$ is determined from (26) as

$J=\left(\frac{1}{q}+\frac{1}{2(1-\cos\varepsilon)}\sum\limits_{n=0}^{\infty}\frac{\bigl{(}P_{n-1}(\cos\varepsilon)-P_{n+1}(\cos\varepsilon)\bigr{)}^{2}}{(2n+1)\mu_{n}^{(p)}}\right)^{-1},$
(42)

and we took into account the fact that the area $|\mathcal{T}|$ of the
spherical cap is given explicitly by $|\mathcal{T}|=2\pi
R_{1}^{2}(1-\cos\varepsilon)$. The surface-averaged and volume-averaged
approximations are particularly simple,

$\overline{\tilde{H}_{\rm
app}(p|r)}=J\frac{g_{0}(r)}{\mu_{0}^{(p)}}\frac{1-\cos\varepsilon}{2}$ (43)

and

$\overline{\tilde{H}_{\rm
app}(p)}=J\frac{3R_{1}^{2}D}{p(R_{2}^{3}-R_{1}^{3})}\frac{1-\cos\varepsilon}{2}=J\frac{D|\mathcal{T}|}{p|\Omega|}.$
(44)

In 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 $\varepsilon=\pi$, the matrix $\mathbf{K}$ is the
identity matrix due to the orthonormality of the eigenfunctions
$\\{v_{n}^{(p)}\\}$. Moreover in this case all terms in (37) with $n>0$
vanish, such that the exact spectral solution becomes

$\tilde{H}_{\varepsilon=\pi}(p|\mathbf{x})=\frac{g_{0}(r)}{1+\mu_{0}^{(p)}/q},$
(45)

where $g_{0}(r)$ and $\mu_{0}^{(p)}$ in (31) and (32) are explicitly given by

$\displaystyle g_{0}(r)$
$\displaystyle=\frac{R_{1}}{r}\frac{e^{\alpha(R_{1}-r)}(e^{2\alpha r}(1+\alpha
R_{2})-e^{2\alpha R_{2}}(1-\alpha R_{2}))}{e^{2\alpha R_{1}}(1+\alpha
R_{2})-e^{2\alpha R_{2}}(1-\alpha R_{2})},$ (46) $\displaystyle\mu_{0}^{(p)}$
$\displaystyle=\frac{e^{2\alpha R_{1}}(1-\alpha R_{1})(1+\alpha
R_{2})-e^{2\alpha R_{2}}(1+\alpha R_{1})(1-\alpha
R_{2})}{R_{1}\bigl{(}e^{2\alpha R_{1}}(1+\alpha R_{2})-e^{2\alpha
R_{2}}(1-\alpha R_{2})\bigr{)}}.$ (47)

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
$\varepsilon\to\pi$. As a consequence, the SCA appears to be accurate not only
in the limit of a vanishingly small target, $\varepsilon\to 0$, but also exact
in the opposite limit $\varepsilon\to\pi$. 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}\to\infty$.

In the limit $R_{2}\to\infty$ 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
$\Omega=\\{\mathbf{x}\in\mathbb{R}^{3}~{}:~{}|\mathbf{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 ${\partial\Omega}$ 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)\to-\frac{k_{n}(\alpha r)}{k_{n}(\alpha R_{1})}.$ (48)

These determine the eigenvalues $\mu_{n}^{(p)}$ 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 PDF
$\tilde{H}(p|\mathbf{x})$ are presented in this geometric setting.

###  3.2 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],

$v_{n}(\theta)=\frac{1}{R_{2}}\sqrt{\frac{2n+1}{4\pi}}P_{n}(\cos\theta)$ (49)

and

$\mu_{n}^{(p)}=g^{\prime}_{n}(R_{2}),$ (50)

with

$g_{n}(r)=\frac{k^{\prime}_{n}(\alpha R_{1})i_{n}(\alpha
r)-i^{\prime}_{n}(\alpha R_{1})k_{n}(\alpha r)}{k^{\prime}_{n}(\alpha
R_{1})i_{n}(\alpha R_{2})-i^{\prime}_{n}(\alpha R_{1})k_{n}(\alpha R_{2})},$
(51)

satisfying $g_{n}(R_{2})=1$ and $g^{\prime}_{n}(R_{1})=0$.222Note that there
is a misprint in equation (B2) of [60], in which $g_{n}^{\prime}(R)$ should be
replaced by $g_{n}^{\prime}(L)$. In turn, the matrix $\mathbf{K}$ remains
unchanged, as well as the formulae (37), (41), and (42) for
$\tilde{H}(p|\mathbf{x})$, $\tilde{H}_{\rm app}(p|\mathbf{x})$, and $J$.

In the limit $R_{1}\to 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

$g_{n}(r)\to\frac{i_{n}(\alpha r)}{i_{n}(\alpha R_{2})},$ (52)

and we therefore recover the SCA result derived in [39] [indeed, equation
(A.21) from [39], together with equations (A.13, A.17, A.18), is identical
with equations (41), (42)]. The remarkable accuracy of the SCA prediction is
illustrated in A.

###  3.3 The PDF in time domain

The solution $H(t|\mathbf{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., $\tilde{H}(p|\mathbf{x})$ has infinitely many poles
$\\{p_{k}\\}$ lying on the negative axis in the complex plane
$p\in\mathbb{C}$. Formally, these poles are determined as zeros of the
determinant of the matrix $\mathbf{M}^{(p)}/q+\mathbf{K}^{(p)}$ when the
matrix is not invertible, yielding a singular, pole-like behaviour of
$\tilde{H}(p|\mathbf{x})$. If $\\{\gamma_{n}^{(p)}\\}$ are the eigenvalues of
this matrix, then the poles correspond to such values of $p$ at which at least
one $\gamma_{n}^{(p)}$ is zero. The poles determine the eigenvalues of the
Laplace operator: $\lambda_{n}=-p_{n}/D$. Evaluating the residues of
$\tilde{H}(p|\mathbf{x})$, one formally gets the spectral expansion of the
solution in time domain:

$H(t|\mathbf{x})=\sum\limits_{n=1}^{\infty}c_{n}(\mathbf{x})e^{-Dt\lambda_{n}},$
(53)

where the $c_{n}(\mathbf{x})$ are determined by these residues.

In practice it is more convenient to search the poles $\\{p^{\prime}_{k}\\}$
of $\tilde{H}_{\rm app}(p|\mathbf{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,

$\frac{1}{2(1-\cos\varepsilon)}\sum\limits_{n=0}^{\infty}\frac{\bigl{(}P_{n-1}(\cos\varepsilon)-P_{n+1}(\cos\varepsilon)\bigr{)}^{2}}{(2n+1)\mu_{n}^{(p)}}=-\frac{1}{q}.$
(54)

The computation of the residue at each pole $p^{\prime}_{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.

Figure 3: Comparison between the exact solution $H(t|\mathbf{x})$ (filled
circles) and the SCA result $H_{\rm app}(t|\mathbf{x})$ (solid line) both of
which are obtained numerically via inverse Laplace transformation of
$\tilde{H}(p|\mathbf{x})$ and $\tilde{H}_{\rm app}(p|\mathbf{x})$ by using the
Talbot algorithm. The latter quantities were calculated with truncation at
$n_{\rm max}=25$. The target is located at the inner sphere, with parameters:
$R_{1}=0.1$, $R_{2}=1$, $r=0.45$, $\theta=0$, and $D=1$. (a) $q=1$,
$\varepsilon=0.1$; (b) $q=1$, $\varepsilon=\pi/2$; (c) $q=100$,
$\varepsilon=0.1$; (d) $q=100$, $\varepsilon=\pi/2$.

Figure 3 illustrates the behaviour of the PDF $H(t|\mathbf{x})$ and its self-
consistent approximation $H_{\rm app}(t|\mathbf{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
$\varepsilon=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.

Figure 4:  Comparison between the exact spectral solution for
$H(t|\mathbf{x})$ (filled circles) and the SCA result $H_{\rm
app}(t|\mathbf{x})$ (solid line) both of which are obtained numerically via
the inverse Laplace transform of $\tilde{H}(p|\mathbf{x})$ and $\tilde{H}_{\rm
app}(p|\mathbf{x})$ by help of the Talbot algorithm. The latter quantities
were calculated with truncation at $n_{\rm max}=25$. The target is located at
the outer sphere, with parameters: $R_{1}=0.1$, $R_{2}=1$, $r=0.45$,
$\theta=0$, and $D=1$. (a) $q=1$, $\varepsilon=0.1$; (b) $q=1$,
$\varepsilon=\pi/2$; (c) $q=100$, $\varepsilon=0.1$; (d) $q=100$,
$\varepsilon=\pi/2$.

##  4 Discussion

In this section we discuss several features of the calculated PDF
$H(t|\mathbf{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 ($\kappa=\infty$) and
finally analyse the "concealment" effect of the inner sphere.

###  4.1 The structure of the PDF in time domain

One may notice that the PDF $H(t|\mathbf{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_{\mathrm{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|\mathbf{x})$ is characterised by the singular Lévy-
Smirnov-type form [28]

$\displaystyle H(t|\mathbf{x})\sim
A_{t}\frac{\exp\left(-\rho^{2}/(4Dt)\right)}{t^{3/2}}\quad(t\to 0),$ (55)

where $A_{t}$ is a computable amplitude and $\rho$ is the shortest distance
between the starting point and the target. This regime is thus dominated by
the "direct" trajectories [63, 64, 62], 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
($\kappa=\infty$), and is an algebraic function of time, $A_{t}\sim\kappa 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|\mathbf{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|\mathbf{x})$
emerges due to the gap between the first and second eigenvalues $\lambda_{1}$
and $\lambda_{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\sim 1/(D\lambda_{1})$ the PDF crosses over to an
exponential decay of the form

$\displaystyle H(t|\mathbf{x})\sim\exp\left(-Dt\lambda_{1}\right).$ (56)

In summary, there is an appreciable amount of the trajectories that lead to
much earlier reaction times, see the discussions in [63, 65, 62, 40, 39]. In
the case of very small target, the mean FRT is close to $1/(D\lambda_{1})$
(see 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_{\rm 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.

###  4.2 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.

Figure 5: Comparison of the FRT PDFs $H(t|\mathbf{x})$ for outer-to-inner case
(Problem I, when the particle starts from the North Pole of the outer sphere
($r=R_{2},~{}\theta=0$) and the target is at the North pole of the inner
sphere; solid blue curve) and inner-to-outer case (Problem II, when the
particle starts from the North pole of the inner sphere
($r=R_{1},~{}\theta=0$) and the target is at the North Pole of the outer
sphere; red dashed curve). Parameters: $q=100$, $D=1$, $R_{1}=0.1$, and
$R_{2}=1$. (a): Areas of the target are equal for problems I and II. The
angular size of the target is $\varepsilon_{1}=0.1$ for I and
$\varepsilon_{2}=0.01$ for II. (b): Angular size of the target site
$\varepsilon=0.1$ is the same for problems I and II.

In figure 5 we compare the FRT PDFs $H(t|\mathbf{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|\mathbf{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.333This effect is related to the
”centrifugal drift” or ”geometric spurious drift” mentioned in [28, 64]. In
the case when $\varepsilon$ 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 $\varepsilon$ 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.

###  4.3 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.

Figure 6: FRT PDF to a target on the outer sphere, with $R_{2}=1$, $D=1$,
$\varepsilon=0.1$, $q=100$ (a), and $q=1$ (b), the starting point at the South
pole of the outer sphere ($r=R_{2}$, $\theta=\pi$), and four values of
$R_{1}$. Note that some curves are not shown for too short times due to
numerical artifacts of the Laplace transform inversion.

In figure 6 we depict $H(t|\mathbf{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.

##  5 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 real-time 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
[80, 81, 82, 63, 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 [63,
62, 39, 40, 48, 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 $\mathbf{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., [91, 92, 93, 94, 84] 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 position-dependent 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 [95, 96, 67]. 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 diffusive-convective motion (with convective motion in
radial directions). All of these points will deserve attention in future work.

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 7 illustrates the remarkable agreement between the prediction
(41) of the SCA and the exact solution (37). Note that both expressions were
truncated at $n_{\rm max}=25$, i.e., the size of the (infinite-dimensional)
matrices $\mathbf{M}^{(p)}$ and $\mathbf{K}$ was set to $(n_{\rm
max}+1)\times(n_{\rm 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_{\rm max}$ from $25$ to $50$ does not
change the solution visibly). The two left panels (a) and (c) show
$\tilde{H}(p|\mathbf{x})$ for a small target ($\varepsilon=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 ($\varepsilon=\pi/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 8 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 ($\varepsilon=\pi/2$, $q=100$), the relative error remains below $5\%$
for a very broad range of $p$-values, in spite of the fact that
$\tilde{H}(p|\mathbf{x})$ varies over more than $10$ orders of magnitude.

Figure 7: Comparison between the SCA result (41) truncated at $n_{\rm max}=25$
(solid line) and the exact solution (37) shown by filled circles, in which the
matrices $\mathbf{M}^{(p)}$ and $\mathbf{K}$ are truncated at the size
$(n_{\rm max}+1)\times(n_{\rm max}+1)$. The target is located on the inner
sphere. The parameters are: $R_{1}=0.1$, $R_{2}=1$, $r=0.45$, $\theta=0$, and
$D=1$. (a) $q=1$, $\varepsilon=0.1$; (b) $q=1$, $\varepsilon=\pi/2$; (c)
$q=100$, $\varepsilon=0.1$; (d) $q=100$, $\varepsilon=\pi/2$. Figure 8:
Relative error of the SCA result (41) as compared to the exact spectral
solution. Parameters are the same as shown in figure 7.

The accuracy of the SCA remains excellent for Problem II when the target is
located on the outer boundary. Figure 9 compares the SCA prediction against
the exact spectral solution, whereas figure 10 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.

Figure 9: Comparison between the SCA result (41) truncated at $n_{\rm max}=25$
(solid line) and the exact solution (37) shown by filled circles, in which the
matrices $\mathbf{M}^{(p)}$ and $\mathbf{K}$ are truncated at the size
$(n_{\rm max}+1)\times(n_{\rm max}+1)$. The target is located at the outer
sphere. The parameters are: $R_{1}=0.1$, $R_{2}=1$, $r=0.45$, $\theta=0$, and
$D=1$. (a) $q=1$, $\varepsilon=0.1$; (b) $q=1$, $\varepsilon=\pi/2$; (c)
$q=100$, $\varepsilon=0.1$; (d) $q=100$, $\varepsilon=\pi/2$. Figure 10:
Relative error of the SCA result (41) as compared to the exact spectral
solution. Parameters are the same as shown in figure 9.

##  Appendix B Mean first-reaction time

The moments of the first-reaction time can be accessed directly from the
generating function

$\mathbb{E}_{\mathbf{x}}\\{\tau^{k}\\}=(-1)^{k}\lim\limits_{p\to
0}\frac{\partial^{k}}{\partial p^{k}}\tilde{H}(p|\mathbf{x}).$ (57)

Once again, it is more convenient to use the SCA result $\tilde{H}_{\rm
app}(p|\mathbf{x})$ for this computation, i.e.,

$\mathbb{E}_{\mathbf{x}}\\{\tau^{k}\\}\approx
T_{k}(\mathbf{x})=(-1)^{k}\lim\limits_{p\to 0}\frac{\partial^{k}}{\partial
p^{k}}\tilde{H}_{\rm app}(p|\mathbf{x}).$ (58)

Since the domain is bounded, the smallest eigenvalue $\mu_{0}^{(p)}$ vanishes
as $p\to 0$ [97], whereas other eigenvalues $\mu_{n}^{(p)}$ converge to
strictly positive limits $\mu_{n}^{(0)}$.

For the target on the inner sphere, one has the following expansion

$\mu_{0}^{(p)}=ap-bp^{2}+O(p^{3}),\quad
a=\frac{R_{2}^{3}-R_{1}^{3}}{3DR_{1}^{2}},\quad
b=\frac{5R_{2}^{6}-9R_{2}^{5}R_{1}+5R_{1}^{3}R_{2}^{3}-R_{1}^{6}}{45D^{2}R_{1}^{3}},$
(59)

while [60]

$\mu_{n}^{(0)}=\frac{n(n+1)}{R_{1}}\frac{1-(R_{1}/R_{2})^{2n+1}}{n+(n+1)(R_{1}/R_{2})^{2n+1}}\quad(n>0).$
(60)

As a consequence, we find

$J=J_{1}p-J_{2}p^{2}+O(p^{3}),\qquad J_{1}=\frac{a}{d},\quad
J_{2}=\frac{a^{2}c+bd}{d^{2}},$ (61)

where $d=(1-\cos\varepsilon)/2$ and

$c=\frac{1}{q}+\frac{1}{2(1-\cos\varepsilon)}\sum\limits_{n=1}^{\infty}\frac{(P_{n-1}(\cos\varepsilon)-P_{n+1}(\cos\varepsilon))^{2}}{(2n+1)\mu_{n}^{(0)}}.$
(62)

Using

$\frac{g_{0}(r)}{\mu_{0}^{(p)}}=\frac{1}{ap}+e+O(p),$ (63)

with

$e=R_{1}^{2}\frac{5r^{3}+10R_{2}^{3}-3r(R_{1}^{2}+6R_{2}^{2})}{10r(R_{2}^{3}-R_{1}^{3})}+\frac{9R_{1}^{4}R_{2}^{2}(R_{2}-R_{1})}{5(R_{2}^{3}-R_{1}^{3})^{2}},$
(64)

we finally get

$T_{1}(\mathbf{x})=\frac{ac}{d}+\frac{b}{a}-ae-\frac{a}{d}\sum\limits_{n>0}\frac{g_{n}^{(p=0)}(r)}{\mu_{n}^{(0)}}P_{n}(\cos\theta)\frac{P_{n-1}(\cos\varepsilon)-P_{n+1}(\cos\varepsilon)}{2},$
(65)

where

$g_{n}^{(p=0)}(r)=(R_{1}/r)^{n+1}\frac{n+(n+1)(r/R_{2})^{2n+1}}{n+(n+1)(R_{1}/R_{2})^{2n+1}}.$
(66)

For the target on the outer sphere, one has

$\mu_{0}^{(p)}=ap-bp^{2}+O(p^{3}),\quad
a=\frac{R_{2}^{3}-R_{1}^{3}}{3DR_{2}^{2}},\quad
b=\frac{R_{2}^{6}-5R_{1}^{3}R_{2}^{3}+9R_{2}R_{1}^{5}-5R_{1}^{6}}{45D^{2}R_{2}^{3}},$
(67)

while [60]

$\mu_{n}^{(0)}=\frac{n(n+1)}{R_{2}}\frac{1-(R_{1}/R_{2})^{2n+1}}{n+1+n(R_{1}/R_{2})^{2n+1}}\quad(n>0).$
(68)

As a consequence, we find that

$J=J_{1}p-J_{2}p^{2}+O(p^{3}),\qquad J_{1}=\frac{a}{d},\quad
J_{2}=\frac{a^{2}c+bd}{d^{2}},$ (69)

with $c$ and $d$ as defined above. Using

$\frac{g_{0}(r)}{\mu_{0}^{(p)}}=\frac{1}{ap}+e+O(p),$ (70)

with

$e=R_{2}^{2}\frac{5r^{3}+10R_{1}^{3}-3r(R_{2}^{2}+6R_{1}^{2})}{10r(R_{2}^{3}-R_{1}^{3})}+\frac{9R_{1}^{2}R_{2}^{4}(R_{2}-R_{1})}{5(R_{2}^{3}-R_{1}^{3})^{2}},$
(71)

we get again equation (65) for $T_{1}(\mathbf{x})$, with

$g_{n}^{(p=0)}(r)=(r/R_{2})^{n}\frac{n+1+n(R_{1}/r)^{2n+1}}{n+1+n(R_{1}/R_{2})^{2n+1}}.$
(72)

##  Appendix C A planar circular annulus domain

The general approach from section 2 can also be applied to a planar circular
annulus domain
$\Omega=\\{\mathbf{x}\in\mathbb{R}^{2}~{}:~{}R_{1}<|\mathbf{x}|<R_{2}\\}$,
with the target represented by an arc $(-\varepsilon,\varepsilon)$ either on
the inner circle, or on the outer circle. This problem is also equivalent to
diffusion between two coaxial cylinders of radii $R_{1}$ and $R_{2}$, (i.e.,
in the domain $\Omega^{\prime}=\Omega\times\mathbb{R}\subset\mathbb{R}^{3}$),
towards a vertical infinitely long partially reactive stripe of angular size
$2\varepsilon$. Note that this shape of the target is different from that
considered in [23, 37, 40].

In this case, the eigenbasis of the Dirichlet-to-Neumann operator is well
known. For instance, for the target on the inner circle, one has [60]:
$v_{n}^{(p)}=e^{in\theta}/\sqrt{2\pi R_{1}}$ and
$\mu_{n}^{(p)}=-g^{\prime}_{n}(R_{1})$, where

$g_{n}(r)=\frac{K^{\prime}_{n}(\alpha R_{2})I_{n}(\alpha
r)-I^{\prime}_{n}(\alpha R_{2})K_{n}(\alpha r)}{K^{\prime}_{n}(\alpha
R_{2})I_{n}(\alpha R_{1})-I^{\prime}_{n}(\alpha R_{2})K_{n}(\alpha R_{1})},$
(73)

with $n\in\mathbb{Z}$. One then easily gets

$\mathbf{K}_{n,n^{\prime}}=\int\limits_{-\varepsilon}^{\varepsilon}d(R_{1}\theta)\frac{e^{i(n-n^{\prime})\theta}}{2\pi
R_{1}}=\frac{\sin(n-n^{\prime})\varepsilon}{\pi(n-n^{\prime})}$ (74)

and

$C_{n}=\int\limits_{-\varepsilon}^{\varepsilon}d(R_{1}\theta)\frac{e^{in\theta}}{\sqrt{2\pi
R_{1}}}=\sqrt{2\pi R_{1}}\frac{\sin n\varepsilon}{\pi n}.$ (75)

The exact solution (12) is then

$\tilde{H}(p|\mathbf{x})=\sum\limits_{n,n^{\prime}=-\infty}^{\infty}e^{in\theta}\biggl{[}\bigl{(}\mathbf{M}^{(p)}/q+\mathbf{K}^{(p)}\bigr{)}^{-1}\biggr{]}_{n,n^{\prime}}\frac{\sin
n^{\prime}\varepsilon}{\pi n^{\prime}}\qquad(\mathbf{x}\in{\partial\Omega}),$
(76)

while the self-consistent approximation reads

$\tilde{H}_{\rm
app}(p|\mathbf{x})=J\sum\limits_{n=-\infty}^{\infty}g_{n}(r)e^{in\theta}\frac{1}{\mu_{n}^{(p)}}\frac{\sin
n\varepsilon}{\pi n},$ (77)

with

$J=\left(\frac{1}{q}+\frac{\pi}{\varepsilon}\sum\limits_{n=-\infty}^{\infty}\frac{1}{\mu_{n}^{(p)}}\biggl{(}\frac{\sin
n\varepsilon}{\pi n}\biggr{)}^{2}\right)^{-1}.$ (78)

Note that the above summations can be reduced to nonnegative values of $n$ by
symmetry.

When the target is on the outer circle one has
$v_{n}^{(p)}=e^{in\theta}/\sqrt{2\pi R_{2}}$ and
$\mu_{n}^{(p)}=g^{\prime}_{n}(R_{2})$, and

$g_{n}(r)=\frac{K^{\prime}_{n}(\alpha R_{1})I_{n}(\alpha
r)-I^{\prime}_{n}(\alpha R_{1})K_{n}(\alpha r)}{K^{\prime}_{n}(\alpha
R_{1})I_{n}(\alpha R_{2})-I^{\prime}_{n}(\alpha R_{1})K_{n}(\alpha R_{2})},$
(79)

while relations (76), (77), and (78) remain unchanged.

##  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
$\omega_{q}(\mathbf{s}|\mathbf{s}_{0})$, characterising the probability of the
reaction event at a boundary point $\mathbf{s}$ after the arrival onto the
boundary point $\mathbf{s}_{0}$, can be expressed in terms of the
eigenfunctions of the Dirichlet-to-Neumann operator [46]

$\omega_{q}(\mathbf{s}|\mathbf{s}_{0})=\sum\limits_{n=0}^{\infty}\frac{v_{n}^{(0)}(\mathbf{s})[v_{n}^{(0)}(\mathbf{s}_{0})]^{*}}{1+\mu_{n}^{(0)}/q}.$
(80)

For the considered shell-like domain in which the target covers the whole
inner sphere, the first arrival point $\mathbf{s}_{0}$ can be set to be on the
North pole, without loss of generality. The explicit form (30) of the
eigenfunctions yields

$\omega_{q}(\mathbf{s}|\mathbf{s}_{0})=\frac{1}{4\pi
R_{1}^{2}}\sum\limits_{n=0}^{\infty}\frac{(2n+1)P_{n}(\cos\theta)}{1+\mu_{n}^{(0)}/q},$
(81)

with $\mu_{n}^{(0)}$ given by (60) for $n>0$, and $\mu_{0}^{(0)}=0$.
Integrating over the azimuthal angle $\phi$ and multiplying by $R_{1}^{2}$, we
get the distribution of the polar angle $\theta$ of the reaction point,

$\omega_{q}(\theta)=\frac{\sin\theta}{2}\sum\limits_{n=0}^{\infty}\frac{(2n+1)P_{n}(\cos\theta)}{1+\mu_{n}^{(0)}/q},$
(82)

where $\sin\theta$ comes from using spherical coordinates. The integral of
$\omega_{q}(\theta^{\prime})$ from $0$ to $\theta$ is the probability that the
reaction location is on the spherical cap of angle $\theta$,

$P_{q}(\theta)=\int\limits_{0}^{\theta}d\theta^{\prime}\omega_{q}(\theta^{\prime})=\frac{1}{2}\sum\limits_{n=0}^{\infty}\frac{P_{n-1}(\cos\theta)-P_{n+1}(\cos\theta)}{1+\mu_{n}^{(0)}/q}.$
(83)

Setting $P_{q}(\Theta_{t})=1/2$ determines implicitly the angular size
$\Theta_{t}$ of the spherical cap, on which half of the reaction events occur.
Figure 11(a) illustrates the behaviour of this probability for three values of
$q$.

In addition, one can calculate the mean angle $\Theta_{m}$, which is related
to the mean geodesic distance between the North pole (the first arrival) and
the reaction location,

$\Theta_{m}=\int\limits_{0}^{\pi}d\theta\theta\omega_{q}(\theta)=\frac{1}{2}\sum\limits_{n=0}^{\infty}\frac{(2n+1)}{1+\mu_{n}^{(0)}/q}I_{n},$
(84)

where

$I_{n}=\int\limits_{-1}^{1}dx\arccos(x)P_{n}(x)=\pi\times\left\\{\begin{array}[]{ll}1&(n=0),\\\
-\left(\frac{(n-2)!!}{2^{(n+1)/2}((n+1)/2)!}\right)^{2}&(n\textrm{ odd}),\\\
0&(n\textrm{ even}).\\\ \end{array}\right.$ (85)

We thus obtain

$\Theta_{m}=\frac{\pi}{2}\left(1-\sum\limits_{k=0}^{\infty}\frac{2(2k+1)+1}{1+\mu_{2k+1}^{(0)}/q}c_{k}\right),$
(86)

with

$c_{k}=\biggl{(}\frac{(2k-1)!!}{2^{k+1}(k+1)!}\biggr{)}^{2}\qquad(k=0,1,2,\ldots).$
(87)

In the limit of high reactivity, $q\to\infty$, the spread harmonic measure
density converges to a Dirac distribution around the first arrival point.
Since $c_{k}\propto k^{-3}$, one gets the asymptotic behaviour
$\Theta_{m}\propto\ln(qR_{1})/(qR_{1})$. In other words, for large $q$, the
inverse of $q$ is proportional to the mean geodesic distance
$R_{1}\Theta_{m}$, up to a logarithmic correction. In the opposite limit of
low reactivity, $q\to 0$, the spread harmonic measure density is getting
uniform, and one gets $\Theta_{m}\to\pi/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 $\Theta_{m}$ on $q$ is
illustrated on figure 11(b).

Figure 11: (a) Probability $P_{q}(\theta)$ for the reaction event to occur on
the inner sphere within the spherical cap of angle $\theta$, for four values
of $q$, with $R_{1}=1$ and $R_{2}=2$. Note that
$P_{0}(\theta)=(1-\cos\theta)/2$. The crossing with the horizontal line at
$1/2$ determines the size $R_{1}\Theta_{t}$ of the typical region on which the
half of reaction events occurs. For large $q$, $R_{1}\Theta_{t}\propto 1/q$.
(b) Mean geodesic distance $R_{1}\Theta_{m}$ as function $q$, for $R_{1}=1$
and $R_{2}=2$. The dashed line shows $q\ln(q)$.

## References

## References

  * [1] Berg HC and Purcell EM 1977 Physics of chemoreception Biophys. J. 20 193
  * [2] Lauffenburger DA and Linderman JJ, eds., 1996 Receptors. Models for binding, trafficking, and signalling, (IRL Press, Oxford)
  * [3] Purves D et al, eds., 2001 Neuroscience, 2nd ed., (Sunderland (MA): Sinauer Associates)
  * [4] Bradshaw RA and Dennis EA, eds., 2010 Handbook of Cell Signaling, 2nd ed., (Amsterdam, Netherlands: Academic Press)
  * [5] Alberts B et al 2015 Molecular Biology of the Cell, $6$th ed., (New York: Garland)
  * [6] Hafner A and Rieger H 2018 Spatial Cytoskeleton Organization Supports Targeted Intracellular Transport Biophys. J. 114 1420
  * [7] Mangeat M and Rieger H The Narrow Escape Problem in two-shell spherical domains; arXiv:2104.13125v1
  * [8] Holcman D and Schuss Z 2014 The narrow escape problem SIAM Rev. 56 213
  * [9] Cheviakov AF, Ward MJ and Straube R 2010 An asymptotic analysis of the mean first passage time for narrow escape problems: part II: the sphere Multiscale Model. Simul. 88 36
  * [10] Cheviakov AF, Reimer AS and Ward MJ 2012 Mathematical modeling and numerical computation of narrow escape problems Phys. Rev. E 85 021131
  * [11] Caginalp C and Chen X 2012 Analytical and Numerical Results for an Escape Problem Arch. Rational. Mech. Anal. 203 329-342
  * [12] Marshall JS 2016 Analytical Solutions for an Escape Problem in a Disc with an Arbitrary Distribution of Exit Holes Along Its Boundary J. Stat. Phys. 165 920-952
  * [13] Grebenkov DS 2016 Universal formula for the mean first passage time in planar domains Phys. Rev. Lett. 117 260201
  * [14] Lindsay AE, Bernoff AJ, and Ward MJ 2017 First Passage Statistics for the Capture of a Brownian Particle by a Structured Spherical Target with Multiple Surface Traps, Multiscale Model. Simul. 15 74-109
  * [15] Bernoff AJ and Lindsay AE 2018 Numerical approximation of diffusive capture rates by planar and spherical surfaces with absorbing pores SIAM J. Appl. Math. 78 266-290
  * [16] Bénichou O and Voituriez R 2008 Narrow-escape time problem: time needed for a particle to exit a confining domain through a small window Phys. Rev. Lett. 100 168105
  * [17] Bénichou O and Voituriez R 2014 From first-passage times of random walks in confinement to geometry-controlled kinetics Phys. Rep. 539 225
  * [18] Oshanin G, Tamm M and Vasilyev O 2010 Narrow-escape times for diffusion in microdomains with a particle-surface affinity: Mean-field results J. Chem. Phys. 132 06B607
  * [19] Bénichou O, Grebenkov D S, Levitz P, Loverdo C and Voituriez R 2010 Optimal reaction time for surface-mediated diffusion Phys. Rev. Lett. 105 150606
  * [20] Bénichou O, Grebenkov DS, Levitz P, Loverdo C and Voituriez R 2011 Mean first-passage time of surface-mediated diffusion in spherical domains J. Stat. Phys. 142 657
  * [21] Rupprecht J-F, Bénichou O, Grebenkov D S and Voituriez R 2012 Kinetics of active surface-mediated diffusion in spherically symmetric domains J. Stat. Phys. 147 891
  * [22] Rupprecht J-F, Bénichou O, Grebenkov DS, and Voituriez R 2015 Exit time distribution in spherically symmetric two-dimensional domains J. Stat. Phys. 158 192-230
  * [23] Grebenkov DS and Oshanin G 2017 Diffusive escape through a narrow opening: new insights into a classic problem Phys. Chem. Chem. Phys. 19 2723
  * [24] Berezhkovskii AM, Coppey M, Sealfon SC and Shvartsman S 2008 Cell-to-cell communication: Time and length scales of ligand internalization in cultures of suspended cells J. Chem. Phys. 128 225102
  * [25] Hughes N, Faulkner C, Morris RJ and Tomkins M 2021 Intercellular Communication as a Series of Narrow Escape Problems IEEE Transactions on Molecular, Biological and Multi-Scale Communications 7 89
  * [26] Mejía-Monasterio C, Oshanin G and Schehr G 2011 First passages for a search by a swarm of independent random searchers J. Stat. Mech. P06022
  * [27] Mattos TG, Mejía-Monasterio C, Metzler R and Oshanin G 2012 First passages in bounded domains: When is the mean first passage time meaningful? Phys. Rev. E 86 031143
  * [28] Redner S 2001 A Guide to First Passage Processes (Cambridge: Cambridge University Press), p. 107.
  * [29] Metzler R, Oshanin G and Redner S, eds., 2014 First-passage phenomena and their applications, (World Scientific Publishers: Singapore)
  * [30] Lindenberg K, Metzler R and Oshanin G, eds., 2019 Chemical Kinetics: Beyond the Textbook, (World Scientific Publishers Europe: London)
  * [31] Grebenkov DS 2019 Spectral theory of imperfect diffusion-controlled reactions on heterogeneous catalytic surfaces, J. Chem. Phys. 151 104108
  * [32] Grebenkov DS and Traytak SD 2019 Semi-analytical computation of Laplacian Green functions in three-dimensional domains with disconnected spherical boundaries J. Comput. Phys. 379 91
  * [33] Grebenkov DS 2020 Diffusion toward non-overlapping partially reactive spherical traps: fresh insights onto classic problems J. Chem. Phys. 152 244108
  * [34] Sokolowski TR, Paijmans J, Bossen L, Miedema T, Wehrens M, Becker NB, Kaizu K, Takahashi K, Dogterom M, and Rein ten Wolde P 2019 eGFRD in all dimensions J. Chem. Phys. 150 054108
  * [35] Shoup D, Lipari G and Szabo A 1981 Diffusion-controlled bimolecular reaction rates. The effect of rotational diffusion and orientation constraints Biophys. J. 36 697
  * [36] Oshanin G, Popescu MN and Dietrich S 2017 Active colloids in the context of chemical kinetics J. Phys. A: Math. Theor. 50 134001
  * [37] Grebenkov DS, Metzler R and Oshanin G 2017 Effects of the target aspect ratio and intrinsic reactivity onto diffusive search in bounded domains New J. Phys. 19 103025
  * [38] Grebenkov DS, Metzler R, Oshanin G, Dagdug L, Berezhkovskii AM and Skvortsov AT 2019 Trapping of diffusing particles by periodic absorbing rings on a cylindrical tube J. Chem. Phys. 150 206101
  * [39] Grebenkov D, Metzler R and Oshanin G 2019 Full distribution of first exit times in the narrow escape problem New J. Phys. 21 122001
  * [40] Grebenkov D, Metzler R and Oshanin 2018 Towards a full quantitative description of single-molecule reaction kinetics in biological cells Phys. Chem. Chem. Phys. 20 16393
  * [41] Bénichou O, Moreau M and Oshanin G 2000 Kinetics of stochastically gated diffusion-limited reactions and geometry of random walk trajectories Phys. Rev. E 61 3388
  * [42] Collins FC and Kimball GE 1949 Diffusion-controlled reaction rates J. Coll. Sci. 4 425
  * [43] Lawley SD and Miles CE 2019 Diffusive search for diffusing targets with fluctuating diffusivity and gating J. Nonlinear Sci. 29 2955
  * [44] Lawley SD and Keener JP 2015 A New Derivation of Robin Boundary Conditions through Homogenization of a Stochastically Switching Boundary SIAM J. Appl. Dyn. Sys. 14 1845-1867
  * [45] Grebenkov DS Imperfect diffusion-controlled reactions, in Ref. [30], Ch. $8$
  * [46] Grebenkov DS 2020 Paradigm Shift in Diffusion-Mediated Surface Phenomena Phys. Rev. Lett. 125 078102
  * [47] v Smoluchowski M 1916 Drei Vorträge über Diffusion, Brownsche Molekularbewegung und Koagulation von Kolloidteilchen, Z. Phys. 17, 557
  * [48] Grebenkov DS, Metzler R and Oshanin G 2021 A molecular relay race: sequential first-passage events to the terminal reaction centre in a cascade of diffusion controlled processes New J. Phys. 23 093004
  * [49] Sapoval B (1994) General Formulation of Laplacian Transfer Across Irregular Surfaces Phys. Rev. Lett. 73 3314
  * [50] Sapoval B, Filoche M, and Weibel E (2002) Smaller is better – but not too small: A physical scale for the design of the mammalian pulmonary acinus Proc. Nat. Acad. Sci. USA 99 10411
  * [51] Grebenkov DS, Filoche M, and Sapoval B (2003) Spectral Properties of the Brownian Self-Transport Operator Eur. Phys. J. B 36 221-231
  * [52] Sapoval B, Andrade JS Jr, Baldassari A, Desolneux A, Devreux F, Filoche M, Grebenkov DS, Russ S (2005) New Simple Properties of a Few Irregular Systems Physica A 357 1-17
  * [53] Grebenkov DS (2006) Scaling Properties of the Spread Harmonic Measures Fractals 14 231-243
  * [54] Grebenkov DS (2015) Analytical representations of the spread harmonic measure Phys. Rev. E 91 052108
  * [55] Arendt W A. ter Elst AFM, Kennedy JB and Sauter M 2014 The Dirichlet-to-Neumann operator via hidden compactness J. Funct. Anal. 266 1757-1786
  * [56] Daners D 2014 Non-positivity of the semigroup generated by the Dirichlet-to-Neumann operator Positivity 18 235-256
  * [57] Arendt W and ter Elst AFM 2015 The Dirichlet-to-Neumann Operator on Exterior Domains Potential Anal. 43 313-340
  * [58] Hassell A and Ivrii V 2017 Spectral asymptotics for the semiclassical Dirichlet to Neumann operator, J. Spectr. Theory 7 881-905
  * [59] Girouard A and Polterovich I 2017 Spectral geometry of the Steklov problem J. Spectr. Theory 7 321-359
  * [60] Grebenkov DS 2020 Surface Hopping Propagator: An Alternative Approach to Diffusion-Influenced Reactions, Phys. Rev. E 102 032125
  * [61] Grebenkov DS 2020 Joint distribution of multiple boundary local times and related first-passage time problems with multiple targets J. Stat. Mech. 103205
  * [62] Grebenkov D, Metzler R and Oshanin G 2018 Strong defocusing of molecular reaction times results from an interplay of geometry and reaction control Comm. Chem. 1 96
  * [63] Godec A and Metzler R 2015 _Optimization and universality of Brownian search in quenched heterogeneous media_ Phys. Rev. E 91 052134
  * [64] Godec A and Metzler R 2016 _Universal proximity effect in target search kinetics in the few encounter limit_ Phys. Rev. X 6 041037
  * [65] Godec A and Metzler R 2017 _First passage time statistics for two-channel diffusion_ , J. Phys. A 50 084001
  * [66] Bénichou O, Chevalier C, Klafter J, Meyer B, and Voituriez R (2010) Geometry-controlled kinetics, Nature Chem. 2 472-477
  * [67] Reva M, DiGregorio DA, and Grebenkov DS 2021 _A first-passage approach to diffusion-influenced reversible binding: insights into nanoscale signaling at the presynapse_ , Sci. Rep. 11 5377
  * [68] Singer A, Schuss Z, Holcman D and Eisenberg B 2006 Narrow Escape I J. Stat. Phys. 122 437
  * [69] Singer A, Schuss Z and Holcman D 2006 Narrow Escape II J. Stat. Phys. 122 465
  * [70] Singer A, Schuss Z and Holcman D 2006 Narrow Escape III J. Stat. Phys. 122 491
  * [71] Grebenkov DS and Rupprecht J-F 2017 _The escape problem for mortal walkers_ , J. Chem. Phys. 146 084106
  * [72] Yu J, Xiao J, Ren X, Lao K and Xie XS 2006 Probing Gene Expression in Live Cells, One Protein Molecule at a Time Science 311 1600
  * [73] Elf J, Li G-W and Xie XS 2007 Probing Transcription Factor Dynamics at the Single-Molecule Level in a Living Cell, Science 316, 1192
  * [74] Wang M, Zhang J, Xu H and Golding I 2019 Measuring transcription at a single gene copy reveals hidden drivers of bacterial individuality, Nature Microbiol. 4 2118
  * [75] Weigel AV, Tamkun MM and Krapf D 2013 Quantifying the dynamic interactions between a clathrin-coated pit and cargo molecules, Proc. Natl. Acad. Sci. USA 110 E4591
  * [76] Képès F 2004 Periodic Transcriptional Organization of the E. coli Genome, J. Mol. Biol. 340
  * [77] Kolesov G, Wunderlich Z, Laikova ON, Gelfand MS and Mirny LA 2007 How gene order is influenced by the biophysics of transcription regulation, Proc. Natl. Acad. Sci. USA 104 13948
  * [78] Pulkkinen O and Metzler R 2013 Distance matters: the impact of gene proximity in bacterial gene regulation, Phys. Rev. Lett. 110 198101
  * [79] Pulkkinen O and Metzler R 2015 _Variance-corrected Michaelis-Menten equation predicts transient rates of single-enzyme reactions and response times in bacterial gene regulation_ , Sci. Rep. 5 17820
  * [80] Kühn T, Ihalainen TO, Hyväluoma J, Dross N, Willman SF, Langowski J, Vihinen-Rante M and Timonen J 2011 Protein Diffusion in Mammalian Cell Cytoplasm, PLoS ONE 6 e22962
  * [81] Han D, Korabel N, Chen R, Johnston M, Gavrilova A, Allan VJ, Fedotov S and Wigh TA 2020 Deciphering anomalous heterogeneous intracellular transport with neural networks, eLife 9, e52224
  * [82] Ma J, Do M, Le Gros MA, Peskin CS, Larabell CA, Mori Y and Isascson SA 2020 Strong intracellular signal inactivation produces sharper and more robust signaling from cell membrane to nucleus, PLoS Comp. Biol. 16 e1008356
  * [83] Witzel P, Götz M, Lanoiselée Y, Franosch T, Grebenkov DS, and Heinrich D (2019) Heterogeneities Shape Passive Intracellular Transport Biophys. J. 117 203-213
  * [84] Grebenkov D, Metzler R and Oshanin G 2020 From single-particle stochastic kinetics to macroscopic reaction rates: fastest first-passage time of $N$ random walkers New J. Phys. 22 103004
  * [85] Bray AJ and Blythe RA 2002 Exact Asymptotics for One-Dimensional Diffusion with Mobile Traps Phys. Rev. Lett. 89 150601
  * [86] Oshanin G, Bénichou O, Coppey M and Moreau M 2002 Trapping reactions with randomly moving traps: Exact asymptotic results for compact exploration Phys. Rev. E 66 060101(R)
  * [87] Moreau M, Oshanin G, Bénichou O and Coppey M 2003 Pascal principle for diffusion-controlled trapping reactions Phys. Rev. E 67 045104(R); 2004 Lattice theory of trapping reactions with mobile species Phys. Rev. E 69 046101
  * [88] Bray AJ, Majumdar SN and Blythe RA 2003 Formal solution of a class of reaction-diffusion models: Reduction to a single-particle problem Phys. Rev. E 67 060102(R)
  * [89] Yuste SB, Oshanin G, Lindenberg K, Bénichou O and Klafter J 2008 Survival probability of a particle in a sea of mobile traps: A tale of tails Phys. Rev. E 78 021105
  * [90] Le Vot F, Yuste SB, Abad E and Grebenkov DS 2020 First-encounter time of two diffusing particles in confinement Phys. Rev. E 102 032118
  * [91] Meerson B and Redner S 2015 Mortality, redundancy, and diversity in stochastic search Phys. Rev. Lett. 114 198101
  * [92] Reynaud K, Schuss Z, Rouach N and Holcman D 2015 Why so many sperm cells? Commun. Integr. Biol. 8 e1017156
  * [93] Lawley S D and Madrid J B 2020 A probabilistic approach to extreme statistics of Brownian escape times in dimensions $1$, $2$, and $3$ J. Nonlinear Sci. 30 1207
  * [94] Lawley S D 2020 Distribution of extreme first passage times of diffusion J. Math. Biol. 80 2301
  * [95] Grebenkov DS 2019 _A unifying approach to first-passage time distributions in diffusing diffusivity and switching diffusion models_ , J. Phys. A 52 174001
  * [96] Lanoiselée Y, Moutal N, and Grebenkov DS (2018) _Diffusion-limited reactions in dynamic heterogeneous media_ , Nature Commun. 9 4398
  * [97] Grebenkov DS 2019 Probability distribution of the boundary local time of reflected Brownian motion in Euclidean domains Phys. Rev. E 100 062110