# Reversible Target-Binding Kinetics of Multiple Impatient Particles

Denis S. Grebenkov  Laboratoire de Physique de la Matière Condensée (UMR
7643),
CNRS – Ecole Polytechnique, IP Paris, 91120 Palaiseau, France
<EMAIL_ADDRESS>Aanjaneya Kumar  Department of Physics,
Indian Institute of Science Education and Research,
Dr. Homi Bhabha Road, Pune 411008, India
<EMAIL_ADDRESS>
###### Abstract

Certain biochemical reactions can only be triggered after binding of a
sufficient number of particles to a specific target region such as an enzyme
or a protein sensor. We investigate the distribution of the reaction time,
i.e., the first instance when all independently diffusing particles are bound
to the target. When each particle binds irreversibly, this is equivalent to
the first-passage time of the slowest (last) particle. In turn, reversible
binding to the target renders the problem much more challenging and
drastically changes the distribution of the reaction time. We derive the exact
solution of this problem and investigate the short-time and long-time
asymptotic behaviors of the reaction time probability density. We also analyze
how the mean reaction time depends on the unbinding rate and the number of
particles. Our exact and asymptotic solutions are compared to Monte Carlo
simulations.

diffusion-controlled reactions, reversible binding, first-passage times,
extreme statistics

###### pacs:

02.50.-r, 05.60.-k, 05.10.-a, 02.70.Rr

##  I Introduction

Diffusion-controlled reactions play an important role in many chemical and
biological processes. In a typical scenario, particles diffuse in a confining
domain towards a specific target region to react or to trigger a biological
event Lauffenburger ; Alberts ; Redner ; Schuss ; Metzler ; Oshanin ;
Grebenkov07 ; Benichou14 ; Holcman14 . Various aspects of such diffusive
search processes have been investigated, with the particular emphasis on
first-passage times (FPTs) that characterize how fast a single particle finds
a single target. The distribution of the first-passage time $\tau$ is usually
described by the survival probability, $S(t)=\mathbb{P}\\{\tau>t\\}$, or,
equivalently, by the probability density function $H(t)=-dS(t)/dt$. The
distribution and, particularly, the mean value $\langle\tau\rangle$ and the
associate reaction rate, have been thoroughly analyzed Grigoriev02 ; Singer06a
; Singer06b ; Singer06c ; Condamin07 ; Benichou08 ; Benichou10 ; Benichou10b ;
Pillay10 ; Cheviakov10 ; Grebenkov10 ; Cheviakov12 ; Caginalp12 ; Mattos12 ;
Berezhkovsky12 ; Rupprecht15 ; Godec16 ; Godec16b ; Grebenkov16 ; Marshall16 ;
Grebenkov17b ; Lanoiselee18 ; Grebenkov18a ; Grebenkov18 ; Sposini19 ;
Grebenkov19 ; Levernier19 ; Grebenkov19b ; Hartich19 ; Hartich19b ;
Grebenkov20a ; Grebenkov20b .

As the diffusive search is typically long, many independent searchers are
generally involved to speed up this process. In this setting, the arrival of
the fastest particle among $N$ particles can trigger the reaction. If
$\tau_{1},\ldots,\tau_{N}$ denote the FPTs of these particles, the fastest
first-passage time (fFPT) is
$\mathcal{T}^{0}_{1,N}=\min\\{\tau_{1},\ldots,\tau_{N}\\}$. As the particles
search independently, the distribution of the fFPT is simply

$\mathbb{P}\\{\mathcal{T}^{0}_{1,N}>t\\}=\mathbb{P}\\{\tau_{1}>t\\}\cdots\mathbb{P}\\{\tau_{N}>t\\}=[S(t)]^{N}.$
(1)

Similarly, the first-passage time $\mathcal{T}^{0}_{K,N}$ of the $K$-th
fastest particle to arrive onto the target is governed by the law

$\mathbb{P}\\{\mathcal{T}^{0}_{K,N}>t\\}=\sum\limits_{j=0}^{K-1}{N\choose
j}[S(t)]^{N-j}[1-S(t)]^{j},$ (2)

where ${N\choose j}$ is the binomial coefficient. The associated probability
density follows immediately:

$\displaystyle H_{K,N}^{0}(t)$ $\displaystyle=$
$\displaystyle-\frac{d\mathbb{P}\\{\mathcal{T}^{0}_{K,N}>t\\}}{dt}$ (3)
$\displaystyle=$ $\displaystyle K{N\choose K}[S(t)]^{N-K}[1-S(t)]^{K-1}H(t).$

While these expressions fully characterize the random variable
$\mathcal{T}^{0}_{K,N}$, finding the large-$N$ asymptotic behavior of its
moments, $\langle[\mathcal{T}^{0}_{K,N}]^{m}\rangle$, is a difficult problem.
More generally, random variables $\mathcal{T}^{0}_{K,N}$ present an example of
extreme value statistics Majumdar20 .

This problem was first studied by Weiss et al. who showed by analyzing the
exact form of $S(t)$ for one-dimensional diffusion on an interval that the
mean $\langle\mathcal{T}^{0}_{K,N}\rangle$ decreases logarithmically slowly
with $N$: $\langle\mathcal{T}^{0}_{K,N}\rangle\propto 1/\ln N$ as $N\to\infty$
Weiss83 . They also briefly considered higher-order moments and argued the
universality of the logarithmic decay for other diffusive processes. This
seminal work was further extended by several authors Basnayake18 ; Basnayake19
; Lawley20a ; Lawley20b ; Reynaud15 ; Schuss19 ; Madrid20 ; Grebenkov20 . For
instance, Basnayake et al. as well as Lawley and Madrid gave rigorous
mathematical proofs for the asymptotic behavior of these moments Basnayake18 ;
Basnayake19 ; Lawley20a  (see also Appendix A.1 for new results concerning the
behavior of the mean of the slowest FPT $\mathcal{T}^{0}_{N,N}$). In addition,
Lawley found the parameters of the asymptotic Gumbel distribution of
$\mathcal{T}^{0}_{1,N}$ for a large class of diffusion processes Lawley20b .
Moreover, this result was extended to the $K$-th fastest FPT
$\mathcal{T}^{0}_{K,N}$ and the asymptotic form of the joint distribution of
$\\{\mathcal{T}^{0}_{1,N},\ldots,\mathcal{T}^{0}_{K,N}\\}$ was derived. The
logarithmic scaling of the mean fFPT was evoked to rationalize the redundancy
in the number of searchers in some biological systems, such as the large
number of sperm cells Reynaud15 ; Schuss19 . We stress, however, that a
logarithmic speed up of the search process is too costly from a practical
point of view; for instance, a tenfold reduction of the mean time would
require more than twenty thousands of particles. Note that if the starting
positions of $N$ particles are uniformly distributed in the domain, the
logarithmic decay is replaced by a much faster power law decay:
$\langle\mathcal{T}^{0}_{1,N}\rangle\propto 1/N$ for partially reactive
targets, and $\langle\mathcal{T}^{0}_{1,N}\rangle\propto 1/N^{2}$ for
perfectly reactive targets Grebenkov20 . More generally, the power law decay
$1/N$ was shown to emerge as a transient regime of moderately large $N$ if the
target is small or if finding the target requires escaping a potential well
Madrid20 . Moreover, as the mean value is not always the relevant time scale
of the process Grebenkov17b ; Grebenkov18 ; Reva21 , evolutionary optimization
of diffusive search does not necessarily aim at reducing the mean value of the
fFPT.

In the above discussion, each particle that arrived onto the target was
supposed to remain on it forever. In particular, the number $\mathcal{N}(t)$
of particles bound to the target at time $t$, is a non-decreasing stochastic
process that increases by $1$ at the arrival of each new particle. As a
consequence, the $K$-th fastest FPT $\mathcal{T}^{0}_{K,N}$ is equal to the
first instance $\mathcal{T}_{K,N}=\inf\\{t>0~{}:~{}\mathcal{N}(t)=K\\}$ when
$K$ particles among $N$ are bound to the target. In many chemical and
biological settings, however, binding to the target is a reversible process,
i.e., each particle remains on the target for some waiting time, unbinds from
it and resumes its bulk diffusion. The waiting time is usually considered to
be an independent random variable obeying an exponential law with the rate
$k_{\rm off}$. As each unbinding event diminishes $\mathcal{N}(t)$ by $1$, the
number of bound particles is no longer a non-decreasing process (Fig. 1). Even
though the dynamics of all particles is Markovian (i.e., their positions and
states at time $t$ fully determine the probabilities of their positions and
states in the future), the number $\mathcal{N}(t)$ of bound particles is a
non-Markovian process. The reversible binding does not affect the statistics
of the first instance $\mathcal{T}_{1,N}$ when one particle (the fastest one)
is bound to the target, i.e., Eq. (1) governs the probability law for
$\mathcal{T}_{1,N}=\mathcal{T}^{0}_{1,N}$. In contrast, the first instance
$\mathcal{T}_{K,N}$ for $K$ particles to be bound to the target is no longer
equal to the $K$-th fastest FPT $\mathcal{T}^{0}_{K,N}$. Indeed, before the
binding of the $K$-th particle, some of the previously bound particles can
unbind, and thus $\mathcal{T}_{K,N}\geq\mathcal{T}^{0}_{K,N}$ (the superscript
$0$ highlights that the first-passage times $\mathcal{T}^{0}_{K,N}$ correspond
to irreversible binding with $k_{\rm off}=0$). Even though the particles are
independent, random waiting times spent by these particles on the target
render the characterization of the reaction times $\mathcal{T}_{K,N}$ much
more challenging than that of $\mathcal{T}^{0}_{K,N}$. As reversible binding
allows for some particles to leave the target before the arrival of the
others, they were termed impatient Grebenkov17 . In Ref. Grebenkov17 , the
problem of two impatient particles diffusing on an interval was mapped onto
intermittent diffusion on a square. Solving the latter problem, the mean
$\langle\mathcal{T}_{2,2}\rangle$ was obtained, and the effect of reversible
binding was analyzed. Even this basic case with two particles required
sophisticated analysis.

Figure 1:  Illustration of a simulated process $\mathcal{N}(t)$ that counts
the number of bound particles at time $t$, with $N=3$. Three filled circles
indicate the first-crossing times $\mathcal{T}_{K,N}$ for $K=1,2,3$. On the
bottom, there is a schematic presentation of the state of each of three
particles: free state (thin dashed line) vs bound state (thick solid line).
Three filled circles indicate the first-passage times $\mathcal{T}^{0}_{K,N}$
for the first, the second, and the third particles (in the order of arrival).
While $\mathcal{T}_{1,N}=\mathcal{T}^{0}_{1,N}$, unbinding events imply that
$\mathcal{T}_{K,N}\geq\mathcal{T}^{0}_{K,N}$ for any $K>1$.

Lawley and Madrid proposed a remarkable approximation to the general problem
Lawley19 . Assuming that the first-binding time and the rebinding time (i.e.,
the random time between unbinding of a particle from the target and its next
rebinding) can be approximated by an exponential random variable with some
rate $\nu$, the number $\mathcal{N}(t)$ can be modeled by a Markovian birth-
death process $\bar{\mathcal{N}}(t)$ between $N+1$ states of $0,1,2,\ldots,N$
bound particles:

${0}$${1}$${2\quad\cdots\quad
N-1}$${N}$$\scriptstyle{N\nu}$$\scriptstyle{k_{\rm
off}}$$\scriptstyle{(N-1)\nu}$$\scriptstyle{2k_{\rm
off}}$$\scriptstyle{\nu}$$\scriptstyle{Nk_{\rm off}}$ (4)

(here and throughout the text, bar denotes the quantities corresponding to the
Lawley-Madrid approximation). Introducing an $(N+1)\times(N+1)$-dimensional
matrix $W$ with zero elements except for

$W_{i,i+1}=ik_{\rm off},\quad W_{i+1,i}=(N+1-i)\nu\quad(i=1,2,\ldots,N),$

and $W_{i,i}$ are chosen so that $W$ has zero column sums, the distribution of
the first-crossing time
$\bar{\mathcal{T}}_{K,N}=\inf\\{t>0~{}:~{}\bar{\mathcal{N}}(t)=K\\}$ can be
written as Lawley19

$\mathbb{P}\\{\bar{\mathcal{T}}_{K,N}>t\\}=\sum\limits_{j=1}^{K}\bigl{[}\exp(W^{(K)}t)\bigr{]}_{j,1}\,,$
(5)

where $W^{(K)}$ is the $K\times K$ matrix obtained by retaining the first $K$
columns and $K$ rows from $W$ and discarding everything else. Here the initial
state was assumed to be $0$, i.e., no bound particles. In other words, the
distribution is expressed in terms of the matrix exponential of $W^{(K)}$. The
probability density of $\bar{\mathcal{T}}_{K,N}$ is even simpler:

$\bar{H}_{K,N}(t)=\nu(N-K+1)\bigl{[}\exp(W^{(K)}t)\bigr{]}_{K,1}\,.$ (6)

Finally, the mean time is fully explicit:

$\langle\bar{\mathcal{T}}_{K,N}\rangle=\frac{1}{\nu}\sum\limits_{m=1}^{K}\biggl{(}\frac{1}{b_{m}}+\sum\limits_{j=m+1}^{K}\frac{(k_{\rm
off}/\nu)^{j-m}}{b_{j}}\prod\limits_{i=m}^{j-1}\frac{d_{i}}{b_{i}}\biggr{)},$
(7)

with $b_{m}=N-K+m$ and $d_{m}=K-m$. Lawley and Madrid proved that
$\bar{\mathcal{N}}(t)$ and $\bar{\mathcal{T}}_{K,N}$ are universal bounds to
$\mathcal{N}(t)$ and $\mathcal{T}_{K,N}$:

$\mathcal{N}(t)\geq\bar{\mathcal{N}}(t)~{}~{}\textrm{for
all}~{}t,\qquad\mathcal{T}_{K,N}\leq\bar{T}_{K,N}~{}~{}\textrm{for all}~{}K,$
(8)

which actually means that

$\displaystyle\mathbb{P}\\{\mathcal{N}(t)>K\\}$ $\displaystyle\geq$
$\displaystyle\mathbb{P}\\{\bar{\mathcal{N}}(t)>K\\},$ (9)
$\displaystyle\mathbb{P}\\{\mathcal{T}_{K,N}>t\\}$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\\{\bar{\mathcal{T}}_{K,N}>t\\}.$ (10)

Moreover, when the target region is getting smaller and smaller, these bounds
become more and more accurate.

The Lawley-Madrid approximation (LMA) opens a way to investigate in detail the
role of reversible binding onto the statistics of sophisticated biochemical
processes involving the arrival of several molecules onto the target region.
The prominent example is the signalling process between neurons when the
fusion of a neurotransmitters vesicle with the presynaptic bouton membrane is
triggered by the arrival of five calcium ions onto the sensor protein Alberts
; Berridge03 ; Eggermann12 ; Nakamura15 ; Dittrich13 ; Guerrier16 ; Reva21 .
It was recently shown by extensive simulations that unbinding events
considerably affect the fusion probability Reva21 .

In spite of numerous advantages of the LMA, it relies on a rough assumption
that both first-binding and rebinding times can be modeled by an exponential
random variable. However, the probability density of the rebinding time is in
general more sophisticated; for instance, in the case of a spherical target,
it diverges at short times as $H(t)\propto t^{-1/2}$ (see Appendix B), in
sharp contrast to the assumed exponential density $\nu e^{-\nu t}$ that
behaves as $\nu+O(t)$ as $t\to 0$. This observation suggests that the LMA does
not capture correctly the short-time behavior that can be relevant for some
applications. In this paper, we undertake a systematic study of the problem of
impatient particles in the case $K=N$. Using the renewal approach, we derive
the exact solution for this problem. We deduce the short-time and long-time
behavior of this solution and compare it with the LMA predictions, as well as
to Monte Carlo simulations. The short-time asymptotic analysis is also
extended to arbitrary $K$. We show that the LMA captures the qualitative
behavior of the reaction time distribution at moderate and long times.
However, the LMA overestimates the mean reaction time and the decay time, and
fails at short times. From the practical point of view, the LMA can thus be
used for qualitative estimations but further improvements are necessary for
getting more accurate results.

The paper is organized as follows. In Sec. II, we start with the mathematical
model of impatient particles, derive the exact form of the probability density
of the reaction time $\mathcal{T}_{N,N}$, and analyze its short-time and long-
time asymptotic behavior. We also obtain the mean reaction time
$\langle\mathcal{T}_{N,N}\rangle$ and the decay time $T_{N}$ of the
exponential decrease of the probability density at long times. We illustrate
the obtained results for a relevant example of restricted diffusion towards a
spherical target. Section III is devoted to a systematic comparison of the
exact solution with two reference solutions: the irreversible binding case and
the LMA. In Sec. IV, we summarize our findings, discuss eventual applications,
and provide final remarks and perspectives. Appendices contain technical
details of the asymptotic analysis (Sec. A), summary of formulas for
restricted diffusion between two concentric spheres (Sec. B), numerical
implementation of the exact solution (Sec. C) and description of Monte Carlo
simulations (Sec. D).

##  II Main results

###  II.1 Mathematical model of impatient particles

We consider $N$ independent indistinguishable point-like particles diffusing
with diffusion coefficient $D$ inside a bounded Euclidean domain
$\Omega\subset\mathbb{R}^{d}$. The boundary $\partial\Omega$ of $\Omega$ is
reflecting everywhere, except for a partially reactive target region $\Gamma$.
After hitting the target, a particle can bind to it with some probability
controlled by the reactivity $\kappa$ Collins49 ; Sano79 ; Shoup82 ; Zwanzig90
; Sapoval94 ; Filoche99 ; Benichou00 ; Grebenkov03 ; Berezhkovskii04 ;
Grebenkov06 ; Grebenkov06a ; Reingruber09 ; Lawley15 ; Grebenkov17b ;
Bernoff18b ; Grebenkov19b ; Grebenkov20 . This binding is reversible, i.e.,
the bound particle stays on the target for an independent random waiting time
distributed according to the exponential law with the rate $k_{\rm off}$.
After unbinding from the target, the particle resumes its bulk diffusion from
a random uniformly distributed location on the target boundary, until the next
time it binds to the target. We are interested in computing the probability
density $H_{N}(t|\bm{x}_{0})=H_{N,N}(t|\bm{x}_{0})$ of the first instance
$\mathcal{T}_{N}=\mathcal{T}_{N,N}$ when all of the $N$ particles are bound to
the target (i.e., when the process $\mathcal{N}(t)$ crosses the level $N$ for
the first time). As the arrival of $N$ particles to the target is supposed to
trigger some reaction event, the first-crossing time $\mathcal{T}_{N}$ is
called the reaction time. For simplicity, we assume that all the particles are
initially free (not bound to the target) and start from the same initial
position $\bm{x}_{0}\in\Omega$. These starting assumptions can be easily
relaxed.

###  II.2 Exact solution

To proceed, we introduce the probability $\mathcal{P}_{t}(n|m)$ that starting
from $m$ particles bound to the target at time $0$, there are $n$ particles
bound to the target at time $t$. This probability is hard to compute in
general due to unbinding events. However, there are two particular cases for
which $\mathcal{P}_{t}(n|m)$ can be expressed in terms of a single particle
dynamics. Let $P(t|\bm{x}_{0})$ denote the occupancy probability that an
initially free particle that started from a point $\bm{x}_{0}$ is bound to the
target at time $t$. Since all particles are independent, the probability of
finding $m$ bound particles on the target at time $t$ is

$\mathcal{P}_{t}(m|0)={N\choose
m}[P(t|\bm{x}_{0})]^{m}\big{[}1-P(t|\bm{x}_{0})\big{]}^{N-m},$ (11)

where we have chosen the initial condition that all particles are initially
free. Similarly, if $Q(t)$ denotes the probability that, starting from the
bound state at time $0$, the particle is bound to the target at time $t$, then

$\mathcal{P}_{t}(m|N)={N\choose m}[Q(t)]^{m}\big{[}1-Q(t)\big{]}^{N-m}.$ (12)

The probability density of the reaction time $\mathcal{T}_{N}$ can then be
obtained from a standard renewal equation:

$\mathcal{P}_{t}(N|0)=\int_{0}^{t}dt^{\prime}\,H_{N}(t^{\prime})\,\mathcal{P}_{t-t^{\prime}}(N|N).$
(13)

Switching to Laplace space allows us to get

$\tilde{H}_{N}(p|\bm{x}_{0})=\frac{\mathcal{L}\\{[P(t|\bm{x}_{0})]^{N}\\}}{\mathcal{L}\\{[Q(t)]^{N}\\}}\,,$
(14)

where both $\mathcal{L}$ and tilde denote the Laplace transform, e.g.

$\tilde{f}(p)=\mathcal{L}[f(t)]=\int\limits_{0}^{\infty}dt\,e^{-pt}\,f(t).$

The inversion of the Laplace transform gives the probability density in time
domain:

$H_{N}(t|\bm{x}_{0})=\mathcal{L}^{-1}\biggr{\\{}\frac{\mathcal{L}\\{[P(t|\bm{x}_{0})]^{N}\\}}{\mathcal{L}\\{[Q(t)]^{N}\\}}\biggr{\\}}.$
(15)

The last step consists in relating the probabilities $P(t|\bm{x}_{0})$ and
$Q(t)$ to the first-passage time statistics of a single particle. This can be
done in a standard way by summing contributions according to the number of
unbinding events (see, e.g., Reva21 ). For instance, one finds

$Q(t)=\Psi(t)+\int\limits_{0}^{t}dt_{1}\,\psi(t_{1})\int\limits_{t_{1}}^{t}dt^{\prime}_{1}\,H(t^{\prime}_{1}-t_{1})\,\Psi(t-t^{\prime}_{1})+\ldots,$

where $\Psi(t)=e^{-k_{\rm off}t}$ is the probability of staying in the bound
state up to time $t$, $\psi(t)=-d\Psi(t)/dt=k_{\rm off}e^{-k_{\rm off}t}$ is
the probability density of the associated waiting time, and $H(t)$ is the
probability density of the rebinding time. The first term in the above
equation is the contribution without unbinding. In the second term, the
particle unbinds at time $t_{1}$, diffuses in the bulk until the next
rebinding at time $t^{\prime}_{1}$, and remains bound until time $t$. The
third and next terms correspond to 2, 3, etc. unbinding events. In Laplace
domain, one simply gets

$\displaystyle\tilde{Q}(p)$ $\displaystyle=$
$\displaystyle\tilde{\Psi}(p)+\tilde{\psi}(p)\tilde{H}(p)\tilde{\Psi}(p)+\ldots=\frac{\tilde{\Psi}(p)}{1-\tilde{\psi}(p)\tilde{H}(p)}$
(16) $\displaystyle=$ $\displaystyle\frac{1}{p+k_{\rm
off}(1-\tilde{H}(p))}\,.$

In turn, the occupancy probability $P(t|\bm{x}_{0})$ includes an additional
step of the first-passage to the target that yields:

$\displaystyle\tilde{P}(p|\bm{x}_{0})$ $\displaystyle=$
$\displaystyle\tilde{H}(p|\bm{x}_{0})\,\tilde{Q}(p)$ (17) $\displaystyle=$
$\displaystyle\frac{\tilde{H}(p|\bm{x}_{0})}{p+k_{\rm
off}(1-\tilde{H}(p))}\,,$ (18)

where $\tilde{H}(p|\bm{x}_{0})$ is the Laplace transform of the probability
density $H(t|\bm{x}_{0})$ of the first-passage time to the target when the
particle started from a point $\bm{x}_{0}$. Note that as the particle is
released after unbinding from a uniformly distributed point on the target
boundary $\Gamma$, one also gets

$H(t)=\frac{1}{|\Gamma|}\int\limits_{\Gamma}d\bm{x}_{0}\,H(t|\bm{x}_{0}),$
(19)

where $|\Gamma|$ is the Lebesgue measure of $\Gamma$ (e.g., the area of
$\Gamma$ in the three-dimensional case). In this way, both probabilities
$P(t|\bm{x}_{0})$ and $Q(t)$ are expressed in terms of the first-passage time
probability density $H(t|\bm{x}_{0})$ for a single particle. In the case
$N=1$, comparison of Eqs. (14, 17) yields immediately that
$H_{1}(t|\bm{x}_{0})=H(t|\bm{x}_{0})$, as expected.

###  II.3 Spectral decompositions

As we deal with restricted diffusion in a bounded domain, the probabilities
$P(t|\bm{x}_{0})$ and $Q(t)$ can be formally deduced from their Laplace
transforms by applying the residue theorem. Let $\\{p_{n}\\}$ be the poles of
$\tilde{P}(p|\bm{x}_{0})$ that lie on the negative real axis:
$0=p_{0}>p_{1}\geq p_{2}\geq\searrow-\infty$. According to Eq. (18), these
poles satisfy the equation:

$p_{n}+k_{\rm off}(1-\tilde{H}(p_{n}))=0.$ (20)

Note that since the poles of $\tilde{H}(p|\bm{x}_{0})$ and $\tilde{H}(p)$ are
the same, they cancel each other in Eq. (18) and thus are not included in the
set of poles of $\tilde{P}(p|\bm{x}_{0})$. If all the poles are simple, the
inverse Laplace transform yields

$P(t|\bm{x}_{0})=P_{\infty}+\sum\limits_{n=1}^{\infty}v_{n}(\bm{x}_{0})\,e^{p_{n}t},$
(21)

where $v_{n}(\bm{x}_{0})$ is the residue of $\tilde{P}(p|\bm{x}_{0})$
evaluated at the pole $p_{n}$. The steady-state limit $P_{\infty}$ corresponds
to the pole at $0$, which can be obtained by using the Taylor expansion

$\tilde{H}(p|\bm{x}_{0})=\langle
e^{-p\tau}\rangle_{\bm{x}_{0}}=1-p\langle\tau\rangle_{\bm{x}_{0}}+O(p^{2}),$
(22)

where $\langle\tau\rangle_{\bm{x}_{0}}$ is the mean FPT to the target for a
single particle started from $\bm{x}_{0}$. Similarly,
$\tilde{H}(p)=1-p\langle\tau\rangle+O(p^{2})$, where $\langle\tau\rangle$ is
the mean rebinding time:

$\langle\tau\rangle=\frac{1}{|\Gamma|}\int\limits_{\Gamma}d\bm{x}_{0}\,\langle\tau\rangle_{\bm{x}_{0}}.$
(23)

As a consequence, Eq. (18) implies that

$P_{\infty}=\frac{1}{1+k_{\rm off}\langle\tau\rangle}\,.$ (24)

The mean rebinding time $\langle\tau\rangle$ can be found explicitly by
writing the boundary value problem for the mean FPT:

$\left\\{\begin{array}[]{r l
l}D\Delta\langle\tau\rangle_{\bm{x}_{0}}&=-1&(\bm{x}_{0}\in\Omega),\\\
-D\partial_{n}\langle\tau\rangle_{\bm{x}_{0}}&=\kappa\mathbf{1}_{\Gamma}(\bm{x}_{0})\langle\tau\rangle_{\bm{x}_{0}}&(\bm{x}_{0}\in\partial\Omega),\\\
\end{array}\right.$ (25)

where $\mathbf{1}_{\Gamma}(\bm{x}_{0})$ is the indicator function of $\Gamma$:
$\mathbf{1}_{\Gamma}(\bm{x}_{0})=1$ if $\bm{x}_{0}\in\Gamma$, and $0$
otherwise. Integrating the first relation over $\bm{x}_{0}\in\Omega$ and
applying the Green’s formula, one gets

$\displaystyle-|\Omega|$ $\displaystyle=$
$\displaystyle\int\limits_{\Omega}d\bm{x}_{0}\,D\Delta\langle\tau\rangle_{\bm{x}_{0}}=\int\limits_{\partial\Omega}d\bm{x}_{0}\,D\partial_{n}\langle\tau\rangle_{\bm{x}_{0}}$
$\displaystyle=$
$\displaystyle-\int\limits_{\Gamma}d\bm{x}_{0}\,\kappa\langle\tau\rangle_{\bm{x}_{0}}=-\kappa|\Gamma|\langle\tau\rangle,$

from which

$\langle\tau\rangle=\frac{|\Omega|}{\kappa|\Gamma|}\,,$ (26)

where $|\Omega|$ is the volume of the domain. For a spherical target of radius
$\rho$, the reactivity can be expressed in terms of the forward constant
$k_{\rm on}=\kappa(4\pi\rho^{2}N_{A})$ (with $N_{A}\approx 6.02\cdot
10^{23}~{}\rm{mol}^{-1}$ being the Avogadro number) Shoup82 ; Lauffenburger
so that the mean rebinding time also reads as
$\langle\tau\rangle=N_{A}|\Omega|/k_{\rm on}$. Defining the dimensionless
quantity

$\eta=k_{\rm off}\langle\tau\rangle=\frac{k_{\rm
off}|\Omega|}{\kappa|\Gamma|}=\frac{k_{\rm off}|\Omega|N_{A}}{k_{\rm on}}\,,$
(27)

we simply get $P_{\infty}=1/(1+\eta)$.

In general, the poles are not necessarily simple. In particular, if the
unbinding rate $k_{\rm off}$ is such that $\tilde{H}(-k_{\rm off})=0$, then
$-k_{\rm off}$ is the pole of $\tilde{P}(p|\bm{x}_{0})$ of higher order than
$1$. For instance, if $-k_{\rm off}$ is the pole of order $2$, the
corresponding term in the spectral expansion (21) is of the form $te^{-k_{\rm
off}t}$. As the set of zeros of the function $\tilde{H}(p)$ is discrete, we
will ignore such specific values of the unbinding rate $k_{\rm off}$.

Introducing

$\tilde{P}(p)=\frac{1}{|\Gamma|}\int\limits_{\Gamma}d\bm{x}_{0}\,\tilde{P}(p|\bm{x}_{0})=\tilde{H}(p)\,\tilde{Q}(p),$
(28)

we can express $\tilde{Q}(p)$ from Eq. (16) as

$\tilde{Q}(p)=\frac{1}{p+k_{\rm off}}+\frac{k_{\rm off}}{p+k_{\rm
off}}\tilde{P}(p),$ (29)

which in time domain reads

$Q(t)=e^{-k_{\rm off}t}+k_{\rm off}\int\limits_{0}^{t}dt^{\prime}\,e^{-k_{\rm
off}t^{\prime}}\,P(t-t^{\prime}).$ (30)

This relation implies that $Q(t)$ monotonously decreases from $Q(0)=1$ to
$Q(\infty)=P_{\infty}$ (see Appendix A.4). Substituting Eq. (21) into this
relation, we get

$Q(t)=P_{\infty}+Q_{0}e^{-k_{\rm
off}t}+\sum\limits_{n=1}^{\infty}q_{n}e^{p_{n}t},$ (31)

where

$\displaystyle Q_{0}$ $\displaystyle=$ $\displaystyle
1-P_{\infty}-\sum\limits_{n=1}^{\infty}q_{n}\,,$ $\displaystyle q_{n}$
$\displaystyle=$ $\displaystyle\frac{v_{n}}{1+p_{n}/k_{\rm
off}}\quad(n=1,2,\ldots),$ $\displaystyle v_{n}$ $\displaystyle=$
$\displaystyle\frac{1}{|\Gamma|}\int\limits_{\Gamma}d\bm{x}_{0}v_{n}(\bm{x}_{0}),$

and we assumed that $-k_{\rm off}$ is not the pole. On one hand, evaluating
$\tilde{P}(p)$ at $p=-k_{\rm off}$, one finds

$\tilde{P}(-k_{\rm off})=\frac{P_{\infty}}{-k_{\rm
off}}-\sum\limits_{n=1}^{\infty}\frac{v_{n}}{k_{\rm off}+p_{n}}\,.$

On the other hand, Eq. (18) implies

$\tilde{P}(-k_{\rm off})=\frac{\tilde{H}(-k_{\rm off})}{-k_{\rm off}+k_{\rm
off}(1-\tilde{H}(-k_{\rm off}))}=-\frac{1}{k_{\rm off}}\,,$

yielding $Q_{0}\equiv 0$, and thus

$Q(t)=P_{\infty}+\sum\limits_{n=1}^{\infty}q_{n}e^{p_{n}t}.$ (32)

In summary, Eqs. (15, 21, 32) fully determine the exact form of the
probability density $H_{N}(t|\bm{x}_{0})$ in terms of the first-passage time
statistics $H(t|\bm{x}_{0})$ of a single particle. Even if $H(t|\bm{x}_{0})$
is known explicitly (see an example in Appendix B), a numerical implementation
of this exact solution remains challenging because it involves: finding zeros
$\\{p_{n}\\}$ of Eq. (20), evaluation of the residues at these poles,
computation of spectral expansions (21, 32) and finally the inverse Laplace
transform in Eq. (15). The practical details of this computation are discussed
in Appendix C. At the same time, our exact solution opens a way to investigate
the asymptotic behavior of the exact probability density $H_{N}(t|\bm{x}_{0})$
in a rather general setting. Before turning to this analysis, we discuss the
mean reaction time.

###  II.4 Mean reaction time

The relation (15) allows one to access the moments of the reaction time:

$\langle\mathcal{T}_{N}^{k}\rangle=(-1)^{k}\lim_{p\to
0}\frac{\partial^{k}}{\partial
p^{k}}\,\frac{\mathcal{L}\\{[P(t|\bm{x}_{0})]^{N}\\}}{\mathcal{L}\\{[Q(t)]^{N}\\}}\,.$
(33)

In particular, the mean reaction time is

$\displaystyle\langle\mathcal{T}_{N}\rangle$ $\displaystyle=$
$\displaystyle\lim\limits_{p\to
0}\biggl{(}\frac{\mathcal{L}\\{t[P(t|\bm{x}_{0})]^{N}\\}}{\mathcal{L}\\{[Q(t)]^{N}\\}}$
$\displaystyle-$
$\displaystyle\frac{\mathcal{L}\\{[P(t|\bm{x}_{0})]^{N}\\}\,\mathcal{L}\\{t[Q(t)]^{N}\\}}{(\mathcal{L}\\{[Q(t)]^{N}\\})^{2}}\biggr{)}.$

As both $P(t|\bm{x}_{0})$ and $Q(t)$ tend to $P_{\infty}$ in the long-time
limit, setting $p=0$ in the above Laplace transforms would yield divergence.
To overcome this issue, one can add and subtract the term $P_{\infty}^{N}$ to
each Laplace transform, e.g.,

$\displaystyle\mathcal{L}\\{t[P(t|\bm{x}_{0})]^{N}\\}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\infty}dt\,t\,e^{-pt}\bigl{(}[P(t|\bm{x}_{0})]^{N}-P_{\infty}^{N}+P_{\infty}^{N}\bigr{)}$
$\displaystyle=$ $\displaystyle P_{\infty}^{N}/p^{2}+a_{1}+o(1)\qquad(p\to
0),$

where

$a_{k}=\int\limits_{0}^{\infty}dt\,t^{k}\,\bigl{(}[P(t|\bm{x}_{0})]^{N}-P_{\infty}^{N}\bigr{)}.$
(34)

Introducing also

$b_{k}=\int\limits_{0}^{\infty}dt\,t^{k}\,\bigl{(}[Q(t)]^{N}-P_{\infty}^{N}\bigr{)},$
(35)

we compute the above limit as

$\langle\mathcal{T}_{N}\rangle=\frac{b_{0}-a_{0}}{P_{\infty}^{N}}=\int\limits_{0}^{\infty}dt\,\frac{[Q(t)]^{N}-[P(t|\bm{x}_{0})]^{N}}{P_{\infty}^{N}}\,.$
(36)

Higher-order moments of $\mathcal{T}_{N}$ can be expressed in a similar way.
For $N=1$, this relation implies that
$\langle\mathcal{T}_{1}\rangle=\langle\tau\rangle_{\bm{x}_{0}}$, as expected.
Equation (36) is a generalization of the expression for the mean slowest FPT
governed by the probability density in Eq. (3):

$\langle\mathcal{T}^{0}_{N,N}\rangle=\int\limits_{0}^{\infty}dt\bigl{(}1-[1-S(t)]^{N}\bigr{)}.$
(37)

In fact, if there is no unbinding ($k_{\rm off}=0$), one gets $Q(t)=1$,
$P_{\infty}=1$, and $P(t|\bm{x}_{0})=1-S(t|\bm{x}_{0})$. In Appendix A.1, we
derive the large-$N$ asymptotic behavior of this mean time:

$\langle\mathcal{T}^{0}_{N,N}\rangle\propto\langle\tau\rangle\ln
N\qquad(N\to\infty).$ (38)

The unbinding mechanism drastically changes this asymptotic behavior into

$\langle\mathcal{T}_{N}\rangle\propto\frac{(1+k_{\rm
off}\langle\tau\rangle)^{N}}{k_{\rm off}N}\qquad(N\to\infty),$ (39)

i.e., a very slow logarithmic increase turns into exponential growth
controlled by the unbinding rate $k_{\rm off}$ (see Appendix A.2). As a
consequence, when many particles are needed to trigger the reaction, even a
small unbinding rate can considerably alter predictions of the irreversible
setting. We emphasize however, that Eq. (39) captures only the large-$N$
asymptotic behavior and is not applicable at small $N$. In particular, an non-
monotonous dependence of the right-hand side of Eq. (39) on $k_{\rm off}$ and
$N$ is not reproduced for the mean reaction time (see further discussion in
Appendix A.2).

###  II.5 Long-time behavior

The probability density $H_{N}(t|\bm{x}_{0})$ can be formally obtained via the
inverse Laplace transform in Eq. (14) by finding the poles $p_{n,N}$ of the
function $\tilde{H}_{N}(p|\bm{x}_{0})$ in the complex plane and applying the
residue theorem. This is a difficult task, even numerically, especially for
large $N$. We focus therefore on the pole $p_{1,N}$ with the smallest absolute
value that determines the decay time $T_{N}=1/|p_{1,N}|$ of the probability
density at long times:

$H_{N}(t|\bm{x}_{0})\propto e^{-t/T_{N}}\qquad(t\to\infty).$ (40)

As $P(t|\bm{x}_{0})$ admits the spectral decomposition (21) with the poles
$p_{n}$, the poles of the numerator $\mathcal{L}\\{[P(t|\bm{x}_{0})]^{N}\\}$
of Eq. (14) are obtained as all linear combinations of the form
$p_{n_{1}}+p_{n_{2}}+\ldots+p_{n_{N}}$. In particular, the pole with the
smallest absolute value is still $p_{1}$ (apart from the pole at $0$). The
situation is more difficult for the denominator $\mathcal{L}\\{[Q(t)]^{N}\\}$,
for which we are looking not for its poles, but for zeros. Let us search for a
zero of this function:

$\mathcal{L}\\{[Q(t)]^{N}\\}(p)=\sum\limits_{n_{1}=0}^{\infty}\cdots\sum\limits_{n_{N}=0}^{\infty}\frac{q_{n_{1}}\cdots
q_{n_{N}}}{p-p_{n_{1}}-\ldots-p_{n_{N}}}\,,$ (41)

where we included the pole at $0$ by setting $p_{0}=0$ and $q_{0}=P_{\infty}$.
As $p\to 0$, the leading term of this expression is $q_{0}^{N}/p$, which can
be separated from the other terms. In the leading-order approximation, one can
set $p=0$ in the remaining terms:

$0=\mathcal{L}\\{[Q(t)]^{N}\\}(p_{1,N})\approx\frac{q_{0}^{N}}{p_{1,N}}-\sum\limits_{\begin{subarray}{c}n_{1}=0,\ldots\\\
\ldots,n_{N}=0,\\\
n_{1}+\ldots+n_{N}>0\end{subarray}}^{\infty}\frac{q_{n_{1}}\cdots
q_{n_{N}}}{p_{n_{1}}+\ldots+p_{n_{N}}}.$ (42)

The multiple sum, from which the term with $n_{1}=n_{2}=\ldots=n_{N}=0$ was
subtracted, can be expressed in terms of an integral, yielding an
approximation for the pole $p_{1,N}$:

$p_{1,N}\approx-
P_{\infty}^{N}\biggl{(}\int\limits_{0}^{\infty}dt\bigl{(}[Q(t)]^{N}-P_{\infty}^{N}\bigr{)}\biggr{)}^{-1}.$
(43)

As a consequence, the decay rate is

$T_{N}\approx
P_{\infty}^{-N}\int\limits_{0}^{\infty}dt\bigl{(}[Q(t)]^{N}-P_{\infty}^{N}\bigr{)}.$
(44)

Curiously, this expression is very similar to the expression (36) for the mean
reaction time.

The accuracy of this approximation depends on various parameters such as
$k_{\rm off}$ and $N$. In fact, in order to get the sum in Eq. (42), we
neglected $p_{1,N}$ under the assumption that $|p_{1,N}|$ is much smaller than
$|p_{1}|$. As $N$ increases, the reaction event occurs at longer times, i.e.,
the decay time increases, and the approximation gets more accurate. In turn,
the case $N=1$ is the worst for this approximation (see discussion in Appendix
A.3). Similarly, as $k_{\rm off}$ increases, the particles unbind more often,
the decay time increases, yielding a more accurate approximation. Note that
the approximation (43) can be improved by accounting perturbatively for next-
order corrections.

###  II.6 Short-time behavior

At short times, the main contribution to the probability density of the first-
passage time comes from the particles that follow almost “direct trajectories”
to the target Godec16b ; Basnayake18 ; Grebenkov18 . As a consequence, the
short-time behavior is generally

$H(t|\bm{x}_{0})\approx
C_{\bm{x}_{0}}\,t^{\alpha}\,e^{-\delta^{2}/(4Dt)}\qquad(t\to 0),$ (45)

where $\delta$ is the distance between the starting point $\bm{x}_{0}$ and the
target region $\Gamma$, $t^{\alpha}$ is a power-law correction, and
$C_{\bm{x}_{0}}$ is the prefactor depending on the starting point, the shape
of the domain, and the reactivity $\kappa$. Note that $\delta$ is either the
Euclidean distance (i.e., the length of the shortest interval connecting
$\bm{x}_{0}$ and $\Gamma$), or the geodesic distance along the shortest
curvilinear path from $\bm{x}_{0}$ and $\Gamma$ that bypasses eventual
obstacles. As a rigorous derivation of this relation is beyond the scope of
the paper, we use it as an assumption, under which the following results are
valid (see an example in Appendix B.3).

When the starting point $\bm{x}_{0}$ lies on the target, $\delta=0$ and Eq.
(45) implies $H(t)\approx Ct^{\alpha}$ with
$C=(1/|\Gamma|)\int\nolimits_{\Gamma}d\bm{x}_{0}\,C_{\bm{x}_{0}}$, from which
$\tilde{H}(p)\approx C\Gamma(\alpha+1)p^{-1-\alpha}$ as $p\to\infty$. If
$\alpha>-2$, the leading term in the denominator of Eq. (18) is $p$ that
implies for any $\bm{x}_{0}$:

$\tilde{P}(p|\bm{x}_{0})\approx\frac{\tilde{H}(p|\bm{x}_{0})}{p}\,,$ (46)

which in time domain gives

$P(t|\bm{x}_{0})\approx\int_{0}^{t}dt^{\prime}\,H(t^{\prime}|\bm{x}_{0})=1-S(t|\bm{x}_{0}).$
(47)

When $\bm{x}_{0}\notin\Gamma$, the integral of Eq. (45) yields in the lowest
order

$P(t|\bm{x}_{0})\approx\frac{4DC_{\bm{x}_{0}}}{\delta^{2}}\,t^{\alpha+2}\,e^{-\delta^{2}/(4Dt)}\,,$
(48)

from which Eq. (11) implies

$\mathcal{P}_{t}(N|0)\approx\left(\frac{4DC_{\bm{x}_{0}}}{\delta^{2}}\right)^{N}\,t^{N(\alpha+2)}\,e^{-N\delta^{2}/(4Dt)}.$
(49)

At short times, one has $Q(t)\approx 1$, and thus $\mathcal{P}_{t}(N|N)\approx
1$, so that

$\tilde{H}_{N}(p|\bm{x}_{0})\approx
p\mathcal{L}\big{\\{}\mathcal{P}_{t}(N|0)\big{\\}}\,,$ (50)

which in turn gives us the short-time behavior

$\displaystyle
H_{N}(t|\bm{x}_{0})\approx\frac{\partial{\mathcal{P}_{t}(N|0)}}{\partial{t}}$
$\displaystyle\approx\frac{N\delta^{2}}{4D}\left(\frac{4DC_{\bm{x}_{0}}}{\delta^{2}}\right)^{N}\,t^{N(\alpha+2)-2}\,e^{-N\delta^{2}/(4Dt)}.$
(51)

This leading-order asymptotic relation can be improved by computing the next-
order term in the integral of Eq. (45) that yields the correction $O(t)$ to
Eq. (51), which is still independent of the unbinding rate $k_{\rm off}$. In
turn, $k_{\rm off}$ appears in the correction $O(t^{2})$ by using $Q(t)\approx
1-k_{\rm off}t$ instead of $Q(t)\approx 1$ in the above derivation. The
integral of this expression yields, in the leading order:

$S_{N}(t|\bm{x}_{0})\approx
1-\left(\frac{4DC_{\bm{x}_{0}}}{\delta^{2}}\right)^{N}t^{N(\alpha+2)}\,e^{-N\delta^{2}/(4Dt)}.$
(52)

It is easy to see that for $N=1$, the leading term of the short-time behavior
in Eq. (45) is recovered. We note that the short-time behavior of the reaction
time density $H_{N}(t|\bm{x}_{0})$ is identical to that of the probability
density $H^{0}_{N,N}(t)$ of the first-passage time $\mathcal{T}^{0}_{N,N}$.
This result is independent of the unbinding rate $k_{\rm off}$ because the
probability of an unbinding event is small at times $t\ll k_{\rm off}^{-1}$.
As a consequence, for any $K$, one can approximate the short-time behavior of
$\mathcal{T}_{K,N}$ by that of $\mathcal{T}^{0}_{K,N}$, for which the
probability density is given explicitly by Eq. (3). Substituting here the
short-time asymptotic relations (45, 47, 48), we get then

$\displaystyle H_{K,N}(t)\approx H^{0}_{K,N}(t)$ (53)
$\displaystyle\approx{N\choose
K}\frac{K\delta^{2}}{4D}\left(\frac{4DC_{\bm{x}_{0}}}{\delta^{2}}\right)^{K}t^{K(\alpha+2)-2}\,e^{-K\delta^{2}/(4Dt)},$

which generalizes Eq. (51).

##  III Discussion

In order to illustrate our general results, we consider a relevant example of
restricted diffusion inside a reflecting sphere of radius $R$ towards a
partially reactive spherical target of radius $\rho$ located at the origin.
This geometrical setting is a simplified model of passive diffusion inside the
cytoplasm towards the nucleus. It was also employed to model diffusion of
calcium ions inside a presynaptic bouton towards a calcium-sensing protein
Reva21 . The distribution of the first-passage time of a single particle was
investigated in Grebenkov18 . In the presence of unbinding events, the exact
spectral decompositions for both probabilities $P(t|\bm{x}_{0})$ and $Q(t)$
were derived in Reva21 . Appendix B summarizes former results needed for
studying the problem of impatient particles.

The numerical method for evaluating the probability density
$H_{N}(t|\bm{x}_{0})$ in Eq. (15) is described in Appendix C. To validate the
accuracy of this exact solution, we also performed Monte Carlo simulations, as
described in Appendix D. In the following, we set $\rho=1$ and $D=1$ to fix
the units of length and time. The radius of outer reflecting sphere is set as
$R=10$ so that the target is relatively small. All the particles start from a
fixed point $\bm{x}_{0}$ such that either $|\bm{x}_{0}|=5$ (relatively far
from the target), or $|\bm{x}_{0}|=2$ (relatively close to the target). To
analyze the effect of unbinding events, we fix the reactivity $\kappa=1$ (and
thus the forward constant $k_{\rm on}$) and vary the unbinding rate $k_{\rm
off}$. For $\kappa=1$, the mean rebinding time $\langle\tau\rangle$ in Eq.
(26) is equal to $333$. Setting $k_{\rm off}=0.003$ or $k_{\rm off}=0.03$, we
can thus examine two settings of moderate ($\eta=1$) and strong ($\eta=10$)
unbinding kinetics, respectively. We will compare our exact solution in Eq.
(15) with Monte Carlo simulations, the LMA, the irreversible binding solution,
and the short-time asymptotic relation.

Table 1 presents the mean reaction time for $N=2$ and $N=3$ with two unbinding
rates $k_{\rm off}$, showing an excellent agreement between Eq. (36) and Monte
Carlo simulations. In turn, the LMA overestimates the mean reaction time, the
largest deviation corresponding to stronger unbinding $k_{\rm off}$ and larger
$N$. In addition, Table 1 presents the decay time in the same setting.
Expectedly, our approximation (44) is least accurate for $N=2$ and the small
unbinding rate $k_{\rm off}=0.003$ (see Sec. II.5). At $N=3$, the agreement is
better. Moreover, for faster unbinding with $k_{\rm off}=0.03$, the
approximation (44) is in excellent agreement with the exact values for both
$N=2$ and $N=3$. In contrast, the LMA predictions are much less accurate.

| Mean $\langle\mathcal{T}_{N}\rangle(\times 10^{3})$ | Decay time $T_{N}(\times 10^{3})$
---|---|---
$N$ | $k_{\rm off}$ | Theory | MC | LMA | Theory | Approx. | LMA
$2$ | 0.003 | 1.20 | 1.21 | 1.47 | 1.04 | 0.76 | 1.33
0.03 | 3.90 | 3.94 | 6.45 | 3.87 | 3.83 | 6.42
$3$ | 0.003 | 2.01 | 2.03 | 3.08 | 1.71 | 1.46 | 2.83
0.03 | 27.9 | 28.3 | 81.3 | 28.0 | 27.8 | 81.3
Table 1:  Mean reaction time $\langle\mathcal{T}_{N}\rangle$ and the decay
time $T_{N}$ for restricted diffusion towards a spherical target of radius
$\rho=1$ and reactivity $\kappa=1$, surrounded by a reflecting concentric
sphere of radius $R=10$, for $N$ particles started from $|\bm{x}_{0}|=5$ with
$D=1$ (see Appendix B for details). Monte Carlo (MC) values of
$\langle\mathcal{T}_{N}\rangle$ were estimated from $10^{6}$ realizations (see
Appendix D), its theoretical values were obtained by numerical integration of
Eq. (36), while LMA values were given by Eq. (7). Theoretical values of the
decay time $T_{N}$ were estimated by fitting
$S_{N}(t|\bm{x}_{0})/H_{N}(t|\bm{x}_{0})$ over a selected range of times,
approximate values were obtained by numerical integration of Eq. (44), while
the LMA values were deduced from the inverse of the smallest eigenvalue of the
matrix $-W^{(N)}$. Note that all times in the table should be multiplied by
$10^{3}$.

Figure 2:  Probability density $H_{N}(t|\bm{x}_{0})$ of the reaction time
$\mathcal{T}_{N}$ for restricted diffusion between concentric spheres of radii
$\rho=1$ and $R=10$, with $D=1$, $\kappa=1$, three values of $k_{\rm off}$
(see legend), two starting positions $|\bm{x}_{0}|=5$ (a,b) and
$|\bm{x}_{0}|=2$ (c,d), and two values of $N$: $N=2$ (a,c) and $N=3$ (b,d).
Symbols show empirical histograms from Monte Carlo simulations with $10^{6}$
particles. Thick lines indicate our exact solution (15) evaluated numerically
as described in Appendix C, whereas thin lines show the Lawley-Madrid
approximation (6). Thin gray solid line presents the short-time asymptotic
behavior (51).

###  III.1 Comparison with irreversible binding case

First, we note that the limit of irreversible binding can be achieved by
setting either $k_{\rm off}=0$ or $\kappa=\infty$. In fact, in the latter
case, any particle that unbinds from the target, immediately rebinds and thus
never leaves the target. As a consequence, the natural parameter
characterizing the unbinding kinetics is the dimensionless quantity $\eta$
defined by Eq. (27). When $\eta$ is small, unbinding kinetics is usually
considered as irrelevant. In the following, we consider the irreversible
binding limit by keeping $\kappa$ fixed and setting $k_{\rm off}\to 0$.

For irreversible binding, the short-time behavior of the probability density
$H_{K,N}(t|\bm{x}_{0})$ is given by Eq. (53). In turn, the long-time behavior
follows from the spectral expansion of the probability density
$H(t|\bm{x}_{0})$ of the first-passage time. In fact, as restricted diffusion
occurs in a bounded domain, the governing Laplace operator, $-\Delta$, has a
discrete spectrum, i.e., a countable set of eigenvalues
$0<\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\leq\ldots\nearrow\infty$ that are
associated to $L_{2}(\Omega)$-normalized eigenfunctions $\\{u_{n}(\bm{x})\\}$
forming a complete orthonormal basis in $L_{2}(\Omega)$ Grebenkov13 . As a
consequence, the survival probability admits the standard spectral
decomposition Gardiner

$S(t|\bm{x}_{0})=\sum\limits_{n=1}^{\infty}c_{n}\,u_{n}(\bm{x}_{0})e^{-Dt\lambda_{n}},\quad
c_{n}=\int\limits_{\Omega}d\bm{x}\,u_{n}(\bm{x}),$ (54)

from which Eq. (3) implies

$H^{0}_{N,N}(t|\bm{x}_{0})\approx
N\,D\lambda_{1}\,c_{1}u_{1}(\bm{x}_{0})e^{-Dt\lambda_{1}}\,.$ (55)

When the target is small, one has

$\lambda_{1}\approx\frac{\kappa|\Gamma|}{D|\Omega|}=\frac{1}{D\langle\tau\rangle}\,,$
(56)

where we used Eq. (26) for the mean rebinding time. We get therefore

$H^{0}_{N,N}(t|\bm{x}_{0})\propto
e^{-t/\langle\tau\rangle}\qquad(t\to\infty).$ (57)

One sees that the decay time here, $\langle\tau\rangle$, does not depend on
$N$, in sharp contrast to the exponential growth of $T_{N}$ in Eq. (44) for
reversible binding.

Lawley found that the mean of the FPT $\mathcal{T}^{0}_{K,N}$ was determined
for any fixed $K$ as Lawley20b

$\langle\mathcal{T}^{0}_{K,N}\rangle\approx\frac{C}{\ln N}\qquad(N\to\infty),$
(58)

with some constant $C$ and higher-order corrections $1/(\ln N)^{2}$ depending
on $K$. However, this behavior cannot be applied to $K=N$. In Appendix A.1, we
show that

$\langle\mathcal{T}^{0}_{N,N}\rangle\approx C^{\prime}\ln
N\qquad(N\to\infty),$ (59)

with another constant $C^{\prime}$ determined by the decay time of the
survival probability for a single particle. Even though the mean arrival time
of the slowest particle differs by a factor $(\ln N)^{2}$ from that of the
fastest particle, the need for $N$ particles to trigger the reaction event
does not considerably slow down the irreversible reaction kinetics. This
observation is totally different in the case of reversible binding, for which
the mean reaction time $\langle\mathcal{T}_{N}\rangle$ in Eq. (39) exhibits an
exponential growth with $N$.

###  III.2 Comparison with the LMA

Now, we compare our exact results to the Lawley-Madrid approximation. This
approximation was designed under assumption that the rebinding time
distribution can be approximated by an exponential law:
$S(t)\approx\bar{S}(t)=e^{-\nu t}$, with an appropriate rate $\nu$. There are
two natural choices for this rate. In order to get the correct long-time
behavior of the survival probability, one can set $\nu=D\lambda_{1}$ to match
the leading term of the exact spectral expansion (54). Alternatively, as the
rebinding time $\tau$ is approximated by an exponential law, one can set
$\nu=1/\langle\tau\rangle$. When the target is small and weakly reactive, Eq.
(56) indicates that $1/\langle\tau\rangle$ is close to $D\lambda_{1}$, and
both choices yield the same result. One sees that the approximate equality
$1/\langle\tau\rangle\approx D\lambda_{1}$ ensures the self-consistence of the
Lawley-Madrid approximation and can thus serve as a practical indicator of its
validity. As a consequence, the LMA is expected to capture the long-time
behavior of the probability density $H_{N}(t|\bm{x}_{0})$ in the limit of
small targets. In the remaining part of this section, we assume that the
validity conditions of the LMA are fulfilled and set
$\nu=1/\langle\tau\rangle$.

First, we look at the mean reaction time. Lawley and Madrid analyzed the
asymptotic behavior of their Eq. (7) in two limits: (i) when $K$ is fixed and
$N\to\infty$, in which case
$\langle\bar{\mathcal{T}}_{K,N}\rangle\sim\langle\tau\rangle K/N$, i.e.,
essentially a linear growth with $K$; (ii) when $1/(1+k_{\rm
off}\langle\tau\rangle)<K/N<1$ is fixed, in which case
$\langle\bar{\mathcal{T}}_{K,N}\rangle$ exhibits a very rapid growth Lawley19
. While the limiting case $K=N$ was not discussed, we deduced the asymptotic
behavior of Eq. (7) by using similar tools:

$\langle\bar{\mathcal{T}}_{N,N}\rangle\approx\frac{(1+k_{\rm
off}\langle\tau\rangle)^{N}}{k_{\rm off}N}\quad(N\gg 1).$ (60)

This expression coincides with Eq. (39) that we obtained from the exact
solution (36). This highlights that the LMA captures qualitatively the long-
time behavior. In turn, as discussed earlier and illustrated in Table 1, both
Eq. (7) and its asymptotic form (60) overestimate the mean reaction time.

Let us now turn to the approximation (6) of the probability density. Denoting
by $0>\nu_{1}\geq\ldots\geq\nu_{N}$ the negative eigenvalues of the matrix
$W^{(N)}$, one sees that the long-time asymptotic behavior is determined by
the largest eigenvalue $\nu_{1}$,

$\bar{H}_{N,N}(t)\propto e^{-t/\bar{T}_{N}}\qquad(t\to\infty),$ (61)

with $\bar{T}_{N}=-1/\nu_{1}$. In turn, the short-time approximation reads

$\bar{H}_{N,N}(t)\approx\frac{N}{\langle\tau\rangle}\,(t/\langle\tau\rangle)^{N-1}+O(t^{N})\quad(t\to
0),$ (62)

where the lower-order powers of $t$ vanish because of the three-diagonal
structure of the matrix $W^{(N)}$, whereas the prefactor in front of the
leading term $t^{N-1}$ is
$\frac{[(W^{(N)})^{N-1}]_{N,1}}{(N-1)!}=N/\langle\tau\rangle^{N-1}$.
Expectedly, this asymptotic behavior is different from relation (51) derived
from our exact solution. Note that in the limit $k_{\rm off}\to 0$, the
assumed exponential law for the rebinding time implies

$\lim\limits_{k_{\rm off}\to
0}\bar{H}_{N,N}(t)=\frac{N}{\langle\tau\rangle}e^{-t/\langle\tau\rangle}(1-e^{-t/\langle\tau\rangle})^{N-1}.$
(63)

Figure 2 illustrates the behavior of the probability density
$H_{N}(t|\bm{x}_{0})$ for two values $N=2$ and $N=3$, as well as two starting
positions $|\bm{x}_{0}|=5$ and $|\bm{x}_{0}|=2$. First of all, we note that
the probability density is broad, spanning over 4 to 6 orders of magnitude in
time. At short times, the probability density $H_{N}(t|\bm{x}_{0})$ does not
depend on the unbinding rate $k_{\rm off}$, yielding the universal behavior of
the left tail of the distribution given by Eq. (3). Note that the short-time
asymptotic relation (51) is not accurate on the considered range of times but
captures correctly the leading-order term. This relation can be improved by
including next-order corrections. At the timescale $1/k_{\rm off}$, the
unbinding mechanism starts to play a role, yielding deviations from the
irreversible binding case. These deviations are actually visible already at
$t\gtrsim 5$ for $k_{\rm off}=0.03$ and $t\gtrsim 50$ for $k_{\rm off}=0.003$.
As unbinding events slow down the reaction, the right tail of the distribution
is shifted towards longer times as $k_{\rm off}$ increases. In fact, the long-
time decay (40) is determined by the exponential function with the decay time
$T_{N}$ increasing with $k_{\rm off}$. Note that Monte Carlo simulations are
in perfect agreement with the exact solution.

The comparison with the LMA reveals its advantages and limitations. The LMA
captures correctly the behavior of the probability density for moderate and
long times, the agreement being better as $k_{\rm off}$ is smaller. One sees
that the LMA systematically overestimates the decay time that controls the
long-time behavior (see Table 1). Deviations become larger as $N$ and $k_{\rm
off}$ increase. Expectedly, the LMA totally fails at short times. Deviations
are stronger when the particles start closer to the target. In fact, when
$|\bm{x}_{0}|=2$, there is a notable maximum around $t\sim 1$ that is not
captured by the LMA. The most probable time determining the position of this
maximum is several orders of magnitude smaller than the mean reaction time.
This maximum can be relevant for applications when the source of particles is
close to the target (see Reva21  for further discussions). Nevertheless, the
explicit character of the Lawley-Madrid approximation and a much simpler
computation of the probability density via Eq. (6) make it a valuable tool for
a first-step analysis of reversible reactions with multiple particles. Further
improvements of the LMA present an important perspective.

Figure 3 presents a complementary view onto the behavior of the reaction time
$\mathcal{T}_{N}$ by showing its cumulative distribution function
$1-S_{N}(t|\bm{x}_{0})$.

Figure 3:  Cumulative distribution function
$\mathbb{P}\\{\mathcal{T}_{N}<t\\}=1-S_{N}(t|\bm{x}_{0})$ of the reaction time
$\mathcal{T}_{N}$ for restricted diffusion between concentric spheres of radii
$\rho=1$ and $R=10$, with $D=1$, $\kappa=1$, three values of $k_{\rm off}$
(see legend), two starting positions $|\bm{x}_{0}|=5$ (a,b) and
$|\bm{x}_{0}|=2$ (c,d), and two values of $N$: $N=2$ (a,c) and $N=3$ (b,d).
Symbols show the empirical cumulative distribution function of reaction times
from Monte Carlo simulations with $10^{6}$ particles. Thick lines indicate the
integral of our exact solution (15) evaluated numerically as described in
Appendix C, whereas thin lines show the Lawley-Madrid approximation (5). Thin
gray solid line presents the short-time asymptotic behavior (52).

##  IV Conclusions and perspectives

Diffusion-controlled reactions involving multiple particles are abundant and
particularly relevant in biochemistry. The need for a sufficient number of
bound particles can be considered as a sort of protection mechanism against
spontaneous triggering, as well as a mean for reliable control of reactions.
The overwhelming majority of former studies in this field were focused on
first-passage times of a single particle, with a straightforward extension to
the extreme statistics of many particles with irreversible binding to the
target. In turn, the problem of impatient particles with reversible binding
seems to remain unnoticed, in spite of its practical relevance Grebenkov17 .
For instance, five calcium ions have to bind to a calcium-sensing protein to
initiate the release of neurotransmitters for signalling between neurons
Berridge03 ; Eggermann12 ; Nakamura15 ; Dittrich13 ; Guerrier16 ; Reva21 . To
outline the role of unbinding kinetics onto this process, we take the
following estimates from Ref. Reva21 : $R=300$ nm, $k_{\rm on}=6.35\cdot
10^{8}~{}\rm{M}^{-1}~{}\rm{s}^{-1}=6.35\cdot
10^{5}~{}\rm{mol}~{}\rm{m}^{-3}~{}\rm{s}^{-1}$ and $k_{\rm off}=1.57\cdot
10^{4}~{}\rm{s}^{-1}$, from which the mean rebinding time is
$\langle\tau\rangle\approx 0.1$ s, see Eqs. (26, 27). As a consequence,
$\eta=k_{\rm off}\langle\tau\rangle\approx 1.57\cdot 10^{3}\gg 1$, so that one
cannot simply ignore reversible binding that drastically changes the
distribution of the reaction time.

Even if the particles diffuse independently, their randomly “asynchronized”
waiting times on the target render the problem of exact characterization of
the reaction times $\mathcal{T}_{K,N}$ mathematically challenging. The
remarkable work by Lawley and Madrid brought an elegant approximate solution
to this problem Lawley19 . The good accuracy of this approximation, as
reported by its authors, might seem to suggest that this challenging problem
is fully solved. In this paper, we showed that this is far from being the end
of the story.

We focused on the particular case of the first time
$\mathcal{T}_{N}=\mathcal{T}_{N,N}$ when all $N$ particles are bound to the
target. This choice allowed us to derive, for the first time, the exact
complete solution of the problem of impatient particles, i.e., to express the
probability density of the random variable $\mathcal{T}_{N}$ in terms of the
first-passage time distribution of a single particle. This exact solution
revealed some limitations and deficiencies of the LMA. In particular, we
showed that the approximate solution captures the qualitative behavior at
moderate and long times but fails at short times. Moreover, the LMA
overestimates the mean reaction time and the decay time so that its
predictions are inaccurate in some settings. At the same time, the complexity
of the exact solution for $\mathcal{T}_{N,N}$ and yet a fully open problem of
finding the exact solution in the general case $\mathcal{T}_{K,N}$ make the
LMA a valuable tool for the qualitative analysis and preliminary estimations.
Moreover, the accuracy of the LMA is expected to be much higher in the limit
of very small targets. We believe that further improvements of the LMA or
development of alternative methods can bring important insights on the problem
of impatient particles. This is an interesting perspective of the present
work.

We also emphasize that impatient particles offer an excellent example of a
physical problem, for which standard numerical methods may be insufficient for
getting the whole picture. In particular, as the mean reaction time and the
decay time grow exponentially fast with the number of particles, getting the
whole distribution of the reaction time $\mathcal{T}_{N}$ was not possible
even for moderate $N$. For instance, a Monte Carlo simulation with $10^{6}$
realizations used to plot the empirical probability density in Fig. 2 took one
day on a laptop. However, this simulation allowed to get the behavior of
$H_{N}(t|\bm{x}_{0})$ only for a limited range of time scales (e.g., from
$10^{2}$ to $10^{4}$ for $N=2$). Even though parallelization can easily
increase the number of realizations (say, by a factor 100 or 1000), it would
not be enough to get the short-time behavior. As the computational time
explodes with $N$, we could not complete Monte Carlo simulations even for
moderate values of $N$ such as $N=5$ or $N=10$. Here, analytical tools and
approximations are indispensable.

While we focused on the setting when all particles start from the same fixed
point $\bm{x}_{0}$, our exact solution can be easily extended to a more
general case with distinct starting points. Moreover, the starting point of
each particle can also be random. In the case of a uniform distribution of the
starting points, the properties of the fastest FPT $\mathcal{T}^{0}_{1,N}$
were studied in Grebenkov20 . An extension to the reaction time
$\mathcal{T}_{N}$ is straightforward.

The exact expression (36) for the mean reaction time
$\langle\mathcal{T}_{N}\rangle$ opens a way to investigate the role of
different parameters onto the reaction kinetics. A rough approximation allowed
us to access the large-$N$ asymptotic behavior of this quantity. However, the
asymptotic formula (39) lacks an exact prefactor and also fails at small $N$.
More accurate analysis of Eq. (36) could hopefully improve this formula to get
a quantitatively accurate description of the mean reaction time. Its extension
to other reaction times $\langle\mathcal{T}_{K,N}\rangle$ presents an exciting
perspective.

###### Acknowledgements.

DG acknowledges the Alexander von Humboldt Foundation for support within a
Bessel Prize award. AK was supported by the Prime Minister’s Research
Fellowship (PMRF) of the Government of India.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Mathematical details

In this Appendix, we discuss some asymptotic relations and derivations.

###  A.1 Mean reaction time for irreversible binding

The mean fastest FPT $\mathcal{T}^{0}_{1,N}$ and, more generally, the mean
$K$-th fastest FPT $\mathcal{T}^{0}_{K,N}$, were thoroughly investigated in
the irreversible binding case Weiss83 ; Basnayake18 ; Basnayake19 ; Lawley20a
; Lawley20b . For any fixed $K$, the mean value
$\langle\mathcal{T}^{0}_{K,N}\rangle$ behaves universally as $1/\ln N$ in the
large $N$ limit, whereas the higher-order corrections $O(1/(\ln N)^{2})$
depend on $K$. In turn, the asymptotic behavior of the slowest FPT
$\mathcal{T}^{0}_{N,N}$ was not discussed, to our knowledge. In particular,
the former result for any fixed $K$ cannot be applied to the case $K=N$. Here,
we sketch the main steps of this analysis, more rigorous derivations being
beyond the scope of this paper.

As the mean time $\langle\mathcal{T}^{0}_{N,N}\rangle$ is given by Eq. (37),
its asymptotic analysis is reduced to that of the survival probability
$S(t|\bm{x}_{0})$ for a single particle. It is easy to check that the function
$f(t)=1-(1-S(t|\bm{x}_{0}))^{N}$ monotonously decreases from $1$ at $t=0$ to
$0$ as $t\to\infty$. The integral in Eq. (37) can be evaluated by
approximating $f(t)$ by the Heaviside step function $\Theta(t_{N}-t)$, where
$t_{N}$ is chosen by setting $f(t_{N})=\zeta$, with $\zeta$ being around $1/2$
(see below). This equation yields $S(t_{N}|\bm{x}_{0})=1-(1-\zeta)^{1/N}$.
When $N$ is large, the right-hand side of this relation is close to $0$. In
other words, the limit $N\to\infty$ corresponds to large $t_{N}$, for which
the spectral expansion (54) can be truncated to a single term,
$S(t_{N}|\bm{x}_{0})\approx c_{1}u_{1}(\bm{x}_{0})e^{-D\lambda_{1}t_{N}}$. As
a consequence, one gets

$\displaystyle\langle\mathcal{T}^{0}_{N,N}\rangle$ $\displaystyle\approx$
$\displaystyle
t_{N}\approx-\frac{\ln(1-(1-\zeta)^{1/N})-\ln(c_{1}u_{1}(\bm{x}_{0}))}{D\lambda_{1}}$
(64) $\displaystyle\approx$ $\displaystyle\frac{\ln
N+\ln(c_{1}u_{1}(\bm{x}_{0}))-\ln\ln\frac{1}{1-\xi}}{D\lambda_{1}}\,.$

Even though this approximate relation depends on a somewhat arbitrary choice
of $\zeta$ around $1/2$, this dependence is weak and corresponds to the sub-
leading (constant) term, as compared to the leading term $\ln N$. Note that
$1/(D\lambda_{1})$ is the decay time for a single particle which determines
the natural timescale of the problem. Figure 4 illustrates the dependence of
$\langle\mathcal{T}^{0}_{N,N}\rangle$ on $N$ and its large-$N$ asymptotic
behavior (64).

Figure 4:  Mean slowest FPT $\langle\mathcal{T}^{0}_{N,N}\rangle$ for
restricted diffusion between concentric spheres of radii $\rho=1$ and $R=10$,
with $|\bm{x}_{0}|=5$, $D=1$, $\kappa=1$, and $k_{\rm off}=0$ (irreversible
binding). Empty circles show the results of a numerical integration of Eq.
(37), dashed line presents Eq. (64) with $\zeta=0.5$, while solid line
illustrates Eq. (64) with $\zeta=0.428$, which was selected to get the best
agreement.

###  A.2 Mean reaction time

The mean reaction time $\langle\mathcal{T}_{N}\rangle$ is determined by Eq.
(36). We note that the function $[Q(t)]^{N}-[P(t|\bm{x}_{0})]^{N}$
monotonously decreases from $1$ at $t=0$ to $0$ as $t\to\infty$. We also
checked that, for large $N$, this function decreases fast enough to allow for
truncation of the integral at some finite time $t_{N}$, whereas the term
$[P(t|\bm{x}_{0})]^{N}$ is small for $t<t_{N}$ and can be omitted.

One can thus apply the same approximation as in the previous subsection. In
fact, we aim at evaluating $t_{N}$ at which $[Q(t_{N})]^{N}=\zeta$ or,
equivalently, $Q(t_{N})=\zeta^{1/N}$, with some $\zeta$ around $1/2$. As
$\zeta^{1/N}$ is close to $1$, one considers the short-time approximation, for
which $Q(t)=1-k_{\rm off}t+O(t^{2})$ (see Eq. (30) and Appendix A.4). We get
thus $t_{N}=(1-\zeta^{1/N})/k_{\rm off}\approx\ln(1/\zeta)/(Nk_{\rm off})$ as
$N\to\infty$, from which

$\langle\mathcal{T}_{N}\rangle\approx\frac{(1+k_{\rm
off}\langle\tau\rangle)^{N}\ln(1/\zeta)}{k_{\rm off}N}\qquad(N\to\infty),$
(65)

where we used Eq. (24). We stress that the prefactor $\ln(1/\zeta)$ stands
here in front of the leading term, whereas in Sec. A.1, an arbitrary parameter
$\zeta$ appeared only in the sub-leading term in Eq. (64), while the leading
term was universal. This feature highlights the deficiency of the
approximation (65). Figure 5 compares the exact mean reaction time
$\langle\mathcal{T}_{N}\rangle$ and the asymptotic behavior
$P_{\infty}^{-N}/(k_{\rm off}N)$ for several values of $N$. While the overall
behavior is correctly captured, deviations are considerable and depend on the
parameters. The curves shown in Fig. 5 can be interpreted as the dependence of
$\ln(1/\zeta)$ on $k_{\rm off}$ and $N$. Further improvements of this
approximation present an interesting perspective. Note also that the
asymptotic behavior (65) with $\ln(1/\zeta)=1$ is identical to Eq. (60) from
the LMA.

Figure 5:  The rescaled mean reaction time
$\langle\mathcal{T}_{N}\rangle\,k_{\rm off}NP_{\infty}^{N}$ (lines) and the
rescaled decay time $T_{N}\,k_{\rm off}NP_{\infty}^{N}$ (symbols) as functions
of $k_{\rm off}$ for restricted diffusion between concentric spheres of radii
$\rho=1$ and $R=10$, with $|\bm{x}_{0}|=5$, $D=1$, $\kappa=1$, and $N=2$
(circles and solid line), $N=3$ (squares and dashed line), $N=4$ (diamonds and
dash-dotted line) and $N=5$ (triangles and dotted line). Theoretical values of
$\langle\mathcal{T}_{N}\rangle$ and $T_{N}$ were obtained by numerical
integration of Eqs. (36) and (44), respectively.

In a first approximation, one may attempt to set the factor $\ln(1/\zeta)$ to
$1$, as in Eq. (60). The non-monotonous dependence of the right-hand side of
the asymptotic form (65) on $k_{\rm off}$ and $N$ may suggest that the mean
reaction time can be optimized with respect to these parameters. In fact, its
derivative with respect to $N$ vanishes at

$N_{c}=\frac{1}{\ln(1+k_{\rm off}\langle\tau\rangle)}\,,$ (66)

suggesting that $\langle\mathcal{T}_{N}\rangle$ can be minimized with respect
to $N$ when $k_{\rm off}\langle\tau\rangle$ is small enough. Similarly, the
derivative with respect to $k_{\rm off}$ vanishes at

$k_{\rm off,c}=\frac{1}{(N-1)\langle\tau\rangle}\,,$ (67)

suggesting a minimum of $\langle\mathcal{T}_{N}\rangle$. However, this
fictitious optimality results from a rough asymptotic formula (65) and does
not occur when the exact solution (36) is considered (see Fig. 6). This
example illustrates danger of relying on approximate solutions and urges for a
more elaborate analysis of the exact solution.

Figure 6:  Mean reaction time $\langle\mathcal{T}_{N}\rangle$ (symbols) and
its asymptotic form (39) as functions of $N$ for restricted diffusion between
concentric spheres of radii $\rho=1$ and $R=10$, with $|\bm{x}_{0}|=5$, $D=1$,
$\kappa=1$, and two values of $k_{\rm off}$. Theoretical values of
$\langle\mathcal{T}_{N}\rangle$ were obtained by numerical integration of Eq.
(36).

###  A.3 The decay time for $N=1$

As mentioned in Sec. II.5, the approximation (44) of the decay time $T_{N}$ is
least accurate in the case $N=1$. To illustrate this point, we note that the
integral in Eq. (44) can be found explicitly for $N=1$. In fact, one has

$\displaystyle\int\limits_{0}^{\infty}dt\bigl{(}Q(t)-P_{\infty}\bigr{)}=\lim\limits_{p\to
0}\biggl{(}\tilde{Q}(p)-\frac{P_{\infty}}{p}\biggr{)}$
$\displaystyle=\lim\limits_{p\to 0}\biggl{(}\frac{1}{p+k_{\rm
off}(1-\tilde{H}(p))}-\frac{P_{\infty}}{p}\biggr{)}$
$\displaystyle=\lim\limits_{p\to 0}\biggl{(}\frac{1}{p+k_{\rm
off}(p\langle\tau\rangle-p^{2}\langle\tau^{2}/2\rangle+O(p^{3}))}-\frac{P_{\infty}}{p}\biggr{)}$
$\displaystyle=\frac{1}{2}P_{\infty}^{2}k_{\rm off}\langle\tau^{2}\rangle,$

where we used Eq. (24). As a consequence, the approximation (44) reads

$p_{1,1}\approx-\frac{2}{P_{\infty}k_{\rm off}\langle\tau^{2}\rangle}\,.$ (68)

Note that for restricted diffusion between concentric spheres, the second
moment of the rebinding time is known explicitly, see Eq. (76). However, a
similar approximation can be used to estimate the first pole of
$\tilde{Q}(p)$, by expanding $\tilde{H}(p)$ in Eq. (16) up to the second order
in $p$, from which $p_{1}\approx-2/(P_{\infty}k_{\rm
off}\langle\tau^{2}\rangle)$. In other words, $p_{1,1}$ turns out to be
identical to $p_{1}$, thus invalidating the approximation in the case $N=1$.

###  A.4 Behavior of the function $Q(t)$

According to Eq. (30), $Q(t)$ is a monotonously decreasing function. In fact,
the time derivative of Eq. (30) reads

$Q^{\prime}(t)=-k_{\rm off}\bigl{(}Q(t)-P(t)\bigr{)}.$ (69)

Comparing the probabilities $P(t)$ and $Q(t)$ of finding the particle bound to
the target, one realizes that the former includes an additional step of
binding to the target and thus $P(t)\leq Q(t)$, implying $Q^{\prime}(t)\leq
0$. This property can also be deduced in a more formal way. In fact, as $P(t)$
is the convolution of $Q(t)$ and $H(t)$ (see Eq. (28)), its integration by
parts yields

$\displaystyle P(t)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{t}dt^{\prime}\,Q(t-t^{\prime})H(t^{\prime})$
$\displaystyle=$
$\displaystyle-\bigl{[}Q(0)S(t)-S(0)Q(t)\bigr{]}+\int\limits_{0}^{t}dt^{\prime}\,Q^{\prime}(t-t^{\prime})S(t^{\prime}).$

Since $S(0)=Q(0)=1$, one deduces

$Q^{\prime}(t)=-k_{\rm
off}\biggl{(}S(t)-\int\limits_{0}^{t}dt^{\prime}\,S(t^{\prime})Q^{\prime}(t-t^{\prime})\biggr{)},$
(70)

where $S(t)\geq 0$ is the survival probability. Applying a sort of induction
argument, one can check that the right-hand side is negative. Note also that
this relation implies $Q^{\prime}(0)=-k_{\rm off}$ and thus $Q(t)\approx
1-k_{\rm off}t+O(t^{2})$ as $t\to 0$.

##  Appendix B Diffusion between concentric spheres

In this Appendix, we summarize former results needed for evaluating the
probability density of the reaction time $\mathcal{T}_{N}$ for the practically
relevant scenario of particles diffusing in a shell-like domain
$\Omega=\\{\bm{x}\in\mathbb{R}^{3}:\rho<|\bm{x}|<R\\}$ bounded between two
concentric spheres of radii $\rho$ and $R$. The inner sphere is a partially
reactive target with reactivity $\kappa$, whereas the outer sphere is
reflecting. The rotational symmetry of the problem allows for an explicit
solution by separation of variables Redner ; Gardiner ; Carslaw . The first-
passage time distribution was discussed in Grebenkov18 , whereas the exact
solution for the probability $P(t|\bm{x}_{0})$ was given in Reva21 .

###  B.1 First-passage time density

The probability density of the first-passage time can be found by separation
of variables in a standard way (see Grebenkov18  for details). The rotational
symmetry implies that $H(t|\bm{x}_{0})$ and other related quantities depend
only on time $t$ and the radial coordinate $r=|\bm{x}_{0}|$. In Laplace
domain, one has

$\tilde{H}(p|\bm{x}_{0})=\frac{g(r)}{g(\rho)-g^{\prime}(\rho)\frac{D}{\kappa}}\,,$
(71)

where

$g(r)=\frac{R\sqrt{p/D}\cosh\xi-\sinh\xi}{r\sqrt{p/D}}\,,$ (72)

with $\xi=(R-r)\sqrt{p/D}$, and $g^{\prime}(r)$ is given by

$g^{\prime}(r)=\frac{(1-Rrp/D)\sinh\xi-\xi\cosh\xi}{r^{2}\sqrt{p/D}}\,.$ (73)

The moments of the first-passage time can be found as

$\langle\tau^{k}\rangle_{\bm{x}_{0}}=(-1)^{k}\lim\limits_{p\to
0}\frac{\partial^{k}\tilde{H}(p|\bm{x}_{0})}{\partial p^{k}}\,.$ (74)

Setting $|\bm{x}_{0}|=\rho$, one also determines the moments of the rebinding
time, e.g.,

$\langle\tau\rangle=\frac{R^{3}-\rho^{3}}{3\kappa\rho^{2}}$ (75)

and

$\displaystyle\langle\tau^{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{2(R^{3}-\rho^{3})^{2}}{9\kappa^{2}\rho^{4}}$ (76)
$\displaystyle+$
$\displaystyle\frac{2(5R^{6}-9R^{5}\rho+5R^{3}\rho^{3}-\rho^{6})}{45D\kappa\rho^{3}}\,.$

The inversion of the Laplace transform in Eq. (71) by means of the residue
theorem yields:

$H(t|\bm{x}_{0})=\frac{D}{\rho^{2}}\sum\limits_{n=1}^{\infty}\hat{\alpha}_{n}^{2}\,\hat{c}_{n}\,u\bigl{(}\hat{\alpha}_{n},|\bm{x}_{0}|\bigr{)}\,e^{-Dt\hat{\alpha}_{n}^{2}/\rho^{2}},$
(77)

with

$\displaystyle
u(\alpha,r)=\frac{\rho\,\sin\bigl{(}\alpha\frac{R-r}{\rho}\bigr{)}-R\alpha\cos\bigl{(}\alpha\frac{R-r}{\rho}\bigr{)}}{r}\,,$
(78)
$\displaystyle\hat{c}_{n}=-\frac{2\mu\rho^{2}}{\hat{\alpha}_{n}}\biggl{[}(\mu
R(R-\rho)+R^{2}+\rho^{2})$
$\displaystyle\times\hat{\alpha}_{n}\sin(\hat{\alpha}_{n}\beta)+(R(R-\rho)\hat{\alpha}_{n}^{2}-\mu\rho^{2})\cos(\hat{\alpha}_{n}\beta)\biggr{]}^{-1}\,,$

and $\hat{\alpha}_{n}$ (with $n=1,2,\ldots$) denoting the positive solutions
of the trigonometric equation

$\tan(\alpha\beta)=\frac{\alpha\bigl{(}\beta+(1+\beta)\mu\bigr{)}}{1+\mu+(1+\beta)\alpha^{2}}\,,$
(79)

with

$\mu=\kappa\rho/D,\qquad\beta=(R-\rho)/\rho.$ (80)

Note that the survival probability is obtained by integrating Eq. (77):

$S(t|\bm{x}_{0})=\sum\limits_{n=1}^{\infty}\hat{c}_{n}\,u\bigl{(}\hat{\alpha}_{n},|\bm{x}_{0}|\bigr{)}\,e^{-Dt\hat{\alpha}_{n}^{2}/\rho^{2}}.$
(81)

###  B.2 The occupancy probability

In turn, the spectral expansion (21) of the occupancy probability
$P(t|\bm{x}_{0})$ was derived in Reva21 , with
$v_{n}(\bm{x}_{0})=c_{n}u(\alpha_{n},|\bm{x}_{0}|)$, where $u(\alpha,r)$ is
given by Eq. (78), $p_{n}=-\alpha_{n}^{2}D/\rho^{2}$, and

$\displaystyle c_{n}$ $\displaystyle=$
$\displaystyle\frac{2\mu}{\sin(\alpha_{n}\beta)\bigl{(}\alpha_{n}^{2}w_{1}+w_{2}\bigr{)}+\alpha_{n}\cos(\alpha_{n}\beta)\bigl{(}\alpha_{n}^{2}w_{3}+w_{4}\bigr{)}},$

with

$\displaystyle w_{1}$ $\displaystyle=4(1+\beta)+\beta(\beta+\mu(1+\beta)),$
(82a) $\displaystyle w_{2}$
$\displaystyle=2(1+\mu-\lambda(1+\beta))-\lambda\beta^{2},$ (82b)
$\displaystyle w_{3}$ $\displaystyle=\beta(1+\beta),$ (82c) $\displaystyle
w_{4}$ $\displaystyle=\beta(1+\mu-\lambda(1+\beta))-3(\beta+\mu(1+\beta)),$
(82d)

$\lambda=k_{\rm off}\rho^{2}/D$, and $\alpha_{n}$ are strictly positive
solutions of the trigonometric equation

$\sin(\alpha_{n}\beta)=\frac{\bigl{[}\alpha_{n}^{2}(\beta+\mu(1+\beta))-\lambda\beta\bigr{]}\alpha_{n}\cos(\alpha_{n}\beta)}{\alpha_{n}^{4}(1+\beta)+\alpha_{n}^{2}(1+\mu-\lambda(1+\beta))-\lambda}\,.$
(83)

enumerated by $n=1,2,\ldots$. Note that the coefficients $v_{n}$ determining
$Q(t)$ in Eq. (32) are simply $v_{n}=c_{n}\,u(\alpha_{n},\rho)$.

###  B.3 Short-time asymptotic behavior

Here we focus on the short-time behavior of $H(t|\bm{x}_{0})$. As the solution
in Eq. (77) depends only on the radial coordinate $r=|\bm{x}_{0}|$, we replace
$\bm{x}_{0}$ by $r$ in the following expressions. Setting $s=(R-\rho)^{2}p/D$
and $\nu=D/(\kappa\rho)$, we can rewrite $\tilde{H}(p|r)$ as

$\displaystyle\tilde{H}(p|r)=\frac{\beta\sqrt{s}\cosh{\big{(}\sqrt{s}\frac{R-r}{R-\rho}\big{)}}-\sinh{\big{(}\sqrt{s}\frac{R-r}{R-\rho}\big{)}}}{r/\rho}$
(84)
$\displaystyle\times\bigg{(}(\nu+\beta)\sqrt{s}\cosh(\sqrt{s})-(1+\nu-\gamma
s)\sinh(\sqrt{s})\bigg{)}^{-1},$

where $\gamma=\nu\beta^{2}\rho/R$. We first study the case $r=\rho$, for which

$\tilde{H}(p|\rho)=\frac{1-\frac{1}{\beta\sqrt{s}}\tanh(\sqrt{s})}{1+\nu/\beta-\frac{1+\nu-\gamma
s}{\beta\sqrt{s}}\tanh(\sqrt{s})}\,.$ (85)

For large $s$, $\tanh(\sqrt{s})=1+O(e^{-2\sqrt{s}})$, which further implies

$\displaystyle\tilde{H}(p|\rho)$ $\displaystyle\approx$
$\displaystyle\frac{\beta\sqrt{s}-1}{(\nu+\beta)\sqrt{s}-(1+\nu)+\gamma s}$
$\displaystyle=$
$\displaystyle\frac{1}{\gamma}\biggr{(}\frac{A_{1}}{\sqrt{s}+x_{1}}+\frac{A_{2}}{\sqrt{s}+x_{2}}\bigg{)},$

where

$x_{1,2}=\frac{(\nu+\beta)\pm\sqrt{(\nu+\beta)^{2}+4\gamma(1+\nu)}}{2\gamma}\,,$
(86)

and

$A_{1}=\frac{1+\beta x_{1}}{x_{1}-x_{2}}\,,\qquad A_{2}=\frac{1+\beta
x_{2}}{x_{2}-x_{1}}\,.$ (87)

Using the following inverse Laplace transform,

$\mathcal{L}\biggr{\\{}\frac{1}{\sqrt{s}+a}\biggr{\\}}=\frac{1}{\sqrt{\pi\tau}}-ae^{a^{2}\tau}\text{erfc}(a\sqrt{\tau})=f_{a}(\tau)$
(88)

(where $\textrm{erfc}(z)$ is the complementary error function), we find the
short-time approximation

$H(t|\rho)\approx
D\frac{A_{1}f_{x_{1}}\big{(}\frac{Dt}{(R-\rho)^{2}}\big{)}+A_{2}f_{x_{2}}\big{(}\frac{Dt}{(R-\rho)^{2}}\big{)}}{(R-\rho)^{2}\gamma}\,.$
(89)

Using $f_{a}(\tau)\approx 1/\sqrt{\pi\tau}-a+2a^{2}\sqrt{\tau/\pi}+O(\tau)$,
one gets

$\displaystyle H(t|\rho)$ $\displaystyle\approx$
$\displaystyle\frac{\kappa}{\sqrt{\pi Dt}}+\kappa(1/\rho+\kappa/D)$ (90)
$\displaystyle+$ $\displaystyle
2\frac{\kappa(1/\rho+\kappa/D)^{2}}{\sqrt{\pi}}\sqrt{Dt}+O(t).$

For the case of $\rho<r<R$, we have

$\tilde{H}(p|r)\approx\frac{\rho}{r}e^{-\sqrt{s}(r-\rho)/(R-\rho)}\tilde{H}(p|\rho),$
(91)

where we neglected the terms of the order $e^{-2\sqrt{s}}$ and
$e^{-2\sqrt{s}\delta}$, with $\delta=(R-r)/(R-\rho)$. Note that if $r$ is
close to $R$ (i.e., if $\delta$ is very small), the above approximation would
be slightly modified. Since the Laplace transform is expressed as a product of
two terms, the inverse Laplace transform yields the convolution

$H(t|r)\approx\frac{\rho}{r}\int_{0}^{t}dt^{\prime}H(t-t^{\prime}|\rho)\frac{(r-\rho)e^{-(r-\rho)^{2}/(4Dt^{\prime})}}{\sqrt{4\pi
Dt^{\prime 3}}}\,,$ (92)

which can be evaluated using the asymptotic expression for $H(t|\rho)$ to give

$H(t|r)=\frac{\rho\kappa e^{-(r-\rho)^{2}/(4Dt)}}{r\sqrt{\pi
Dt}}\bigg{(}1+\frac{2Dt(1+\kappa\rho/D)}{\rho(r-\rho)}+\cdots\bigg{)}.$ (93)

One can recognize Eq. (45) in the leading term, with $\alpha=-1/2$ and

$C_{\bm{x}_{0}}=\frac{\rho\kappa}{r\sqrt{\pi D}}\,.$ (94)

##  Appendix C Numerical implementation

Our central formula (15) expresses the probability density
$H_{N}(t|\bm{x}_{0})$ in terms of the accessible probabilities
$P(t|\bm{x}_{0})$ and $Q(t)$. However, its practical implementation requires
the computation of two Laplace transforms,
$\mathcal{L}\\{[P(t|\bm{x}_{0})]^{N}\\}$ and $\mathcal{L}\\{[Q(t)]^{N}\\}$,
and then the evaluation of the inverse Laplace transform of their ratio. Since
both $P(t|\bm{x}_{0})$ and $Q(t)$ are given as spectral expansions, such a
computation becomes numerically difficult, especially at small and large times
when $H_{N}(t|\bm{x}_{0})$ rapidly decays. We also attempted a direct solution
of the related deconvolution problem:

$[P(t|\bm{x}_{0})]^{N}=\int\limits_{0}^{t}dt^{\prime}\,H_{N}(t^{\prime}|\bm{x}_{0})\,[Q(t-t^{\prime})]^{N},$
(95)

but it was numerically unstable.

To resolve this difficulty, one can integrate Eq. (95) by parts to transform
it into an integral equation on the survival probability
$S_{N}(t|\bm{x}_{0})=\mathbb{P}\\{\mathcal{T}_{N}>t\\}$:

$\displaystyle[Q(t)]^{N}-[P(t|\bm{x}_{0})]^{N}$ $\displaystyle=$
$\displaystyle S_{N}(t|\bm{x}_{0})$ $\displaystyle-$
$\displaystyle\int\limits_{0}^{t}dt^{\prime}\,S_{N}(t^{\prime}|\bm{x}_{0})\,q_{N}(t-t^{\prime}),$

where we used $Q(0)=S_{N}(0|\bm{x}_{0})=1$, and defined

$q_{N}(t)=-\frac{d}{dt}[Q(t)]^{N}=-N[Q(t)]^{N-1}\frac{dQ}{dt}\,.$ (97)

Considering the last term in Eq. (C) as the application of an integral
operator $\mathcal{Q}$ to the function $S_{N}(t|\bm{x}_{0})$, one can formally
invert this relation to get

$S_{N}(t|\bm{x}_{0})=(I-\mathcal{Q})^{-1}\bigl{(}[Q(t)]^{N}-[P(t|\bm{x}_{0})]^{N}\bigr{)},$
(98)

where $I$ is the identity operator. Expanding the operator
$(I-\mathcal{Q})^{-1}$ into the geometric series, one finally expresses the
survival probability as

$S_{N}(t|\bm{x}_{0})=\int\limits_{0}^{t}dt^{\prime}\bigl{(}[Q(t^{\prime})]^{N}-[P(t^{\prime}|\bm{x}_{0})]^{N}\bigr{)}R(t-t^{\prime}),$
(99)

where

$\displaystyle R(t)$ $\displaystyle=$
$\displaystyle(I-\mathcal{Q})^{-1}\delta(t)$ $\displaystyle=$
$\displaystyle\delta(t)+q_{N}(t)+\int\limits_{0}^{t}dt_{1}\,q_{N}(t_{1})\,q_{N}(t-t_{1})+\ldots,$

i.e., the sum of convolutions of $q_{N}(t)$ with itself of all orders.

In practice, we compute both $[Q(t)]^{N}-[P(t|\bm{x}_{0})]^{N}$ and $q_{N}(t)$
over a linear grid of points $\\{0,\delta,2\delta,\ldots,\delta(K-1)\\}$ and
then evaluate convolutions by fast Fourier transform (FFT). In this way, one
gets the survival probability evaluated at grid points:

$S_{N}(j\delta|\bm{x}_{0})\approx\mathcal{F}^{-1}\left\\{\frac{\mathcal{F}\\{p_{k}\\}}{1-\mathcal{F}\\{q_{k}\\}}\right\\}\,,$
(101)

with $j=0,1,\ldots,K-1$, where
$p_{k}=c_{k}\bigl{(}[Q(k\delta)]^{N}-[P(k\delta|\bm{x}_{0})]^{N}\bigr{)}$,
$q_{k}=\delta c_{k}q_{N}(k\delta)$, $c_{0}=1/2$, $c_{k}=1$ for $0<k<K$, and
$c_{k}=0$ for $K\leq k<2K$. Here the coefficient $c_{0}$ accounts for the
integration weight $1/2$ of the first point, whereas $c_{k}$ for $k\geq K$
allow one to pad the vectors by $0$ for the proper computation of linear
convolutions via direct ($\mathcal{F}$) and inverse ($\mathcal{F}^{-1}$) FFTs
applied to vectors of length $2K$. Note that the probability density
$H_{N}(t|\bm{x}_{0})$ can also be found via FFT as

$H_{N}(j\delta|\bm{x}_{0})\approx\frac{1}{\delta}\mathcal{F}^{-1}\left\\{\bigl{(}e^{2\pi
ik/(2K)}-1\bigr{)}\frac{\mathcal{F}\\{p_{k}\\}}{1-\mathcal{F}\\{q_{k}\\}}\right\\}\,,$
(102)

in analogy with the evaluation of a derivative via standard Fourier transform:
$f^{\prime}(x)=\mathcal{F}^{-1}\\{ik\mathcal{F}\\{f(x)\\}\\}$.

The time step $\delta$ sets the minimal time at which both
$S_{N}(t|\bm{x}_{0})$ and $H_{N}(t|\bm{x}_{0})$ are available, and controls
the accuracy of the whole computation. In fact, it determines how accurately
discrete sums approximate convolution integrals. This is particularly
important for the evaluation of $\mathcal{F}\\{q_{k}\\}$, whose maximal value
is achieved at the zero frequency:

$\displaystyle\max\\{|\mathcal{F}\\{q_{k}\\}|\\}$ $\displaystyle=$
$\displaystyle\mathcal{F}_{0}\\{q_{k}\\}=\sum\limits_{k=0}^{2K-1}q_{k}\approx\int\limits_{0}^{t_{\rm
max}}dt\,q_{N}(t)$ $\displaystyle=$ $\displaystyle 1-[Q(t_{\rm
max})]^{N}\approx 1-P_{\infty}^{N}<1.$

As a consequence, $1/(1-\mathcal{F}\\{q_{k}\\})$ is well defined. However,
when $N$ or $k_{\rm off}$ increase, the maximum approaches to $1$. If $\delta$
is not small enough, inaccurate discretization may result in
$\mathcal{F}_{0}\\{q_{k}\\}$ exceeding $1$ and thus strong instabilities in
the above computation. For the computation of theoretical curves in Fig. 2, we
used $\delta=0.01$ in all cases, except for the case $k_{\rm off}=0.03$ and
$N=3$, for which $\delta=0.005$ was needed.

##  Appendix D Monte Carlo simulations

Monte Carlo simulations were realized via a standard event-driven scheme. Each
particle was equipped by its internal “clock” $t_{i}$ and the binary state
variable $s_{i}$ indicating whether the particle is bound or not. At time $0$,
all particles are free ($s_{i}=1$) and released from a fixed point
$\bm{x}_{0}$, with their clocks being set to $0$. The particles diffuse
independently and bind the target at random times sampled from the probability
density $H(t|\bm{x}_{0})$. The internal clock of each particle is thus set to
its (individual) first-binding time, while their states are set to $0$
(bound). We emphasize that these FPTs account for partial reactivity of the
target, i.e., for eventual failed binding attempts and reflections from the
target, until the successful binding. Selecting the particle with the minimal
internal clock (say, $t_{i}$), one updates this clock by adding a random
waiting time $\delta_{i}$ generated from the exponential law with the rate
$k_{\rm off}$, and sets its state variable $s_{i}$ to $1$ (free). In other
words, $t_{i}$ is replaced by $t_{i}+\delta_{i}$, which is the instance when
the $i$-th particle unbinds from the target and resumes its diffusion. From
now on, the following step is repeated: one selects the particle with the
minimal internal time (say, $t_{j}$); if $s_{j}=0$ (i.e., at the instance
$t_{j}$ the particle binds to the target), we evaluate the number of bound
particles at time $t_{j}$, and the simulation is stopped if all particles are
bound; if the simulation is not stopped, the clock $t_{j}$ is updated by
adding either a newly generated random waiting time $\delta_{j}$ (if
$s_{j}=0$), or a random rebinding time $\tau$ sampled from the probability
density $H(t)$ (if $s_{j}=1$). This step is repeated until the simulation is
stopped (see Fig. 1).

The first-binding times are generated from the known probability density
$H(t|\bm{x}_{0})$ given by Eq. (77). To sample from a broad distribution
$H(t|\bm{x}_{0})$ spanning several orders of magnitude in time, we first
perform a change of variable $\zeta=\ln t$, and obtain the associated
probability density $H_{1}(\zeta|\bm{x}_{0})$. Prior to running simulations,
we create a linear grid of possible values $\zeta_{k}$, ranging from
$\zeta_{\min}$ to $\zeta_{\rm max}$, with a step $d\zeta=0.01$, and a grid
containing the probability weight $H_{1}(\zeta_{k}|\bm{x}_{0})d\zeta$ of each
value $\zeta_{k}$. Using these probability weights, a (pseudo)-random value of
$\zeta$ is generated by using the Matlab function `randsample`, and the
corresponding first-binding time is obtained as $e^{\zeta}$. The same method
is used for generating rebinding times from the known probability density
$H(t)$.

For the considered example of restricted diffusion between two spheres, the
explicit form of the survival probability $H(t|\bm{x}_{0})$ is provided in
Appendix B. The spectral decomposition (81) was truncated at a large order
$n=10000$ in order to access accurately the short-time behavior of
$H(t|\bm{x}_{0})$. The zeros $\hat{\alpha}_{n}$ of Eq. (79) were found by the
bisection method (see Grebenkov18 ; Reva21  for details). The grid bounds
$\zeta_{\rm min}$ and $\zeta_{\rm max}$ depend on the parameters and were
chosen manually to cover a broad range of times whose probability density is
not negligible (e.g., we used $\zeta_{\rm min}=-15$ and $\zeta_{\rm max}=10$
for computing $H(t)$ for $\kappa=1$).

## References

  * (1) D. A. Lauffenburger and J. Linderman, Receptors: Models for Binding, Trafficking, and Signaling (Oxford University Press, Oxford, 1993).
  * (2) Alberts, B. et al., Molecular Biology of the Cell 5th edn. (Garland Science, Taylor & Francis Group, New York, 2008).
  * (3) S. Redner, A Guide to First Passage Processes (Cambridge: Cambridge University press, 2001).
  * (4) Z. Schuss, Brownian Dynamics at Boundaries and Interfaces in Physics, Chemistry and Biology (Springer, New York, 2013).
  * (5) R. Metzler, G. Oshanin, and S. Redner (Eds.) First-Passage Phenomena and Their Applications (Singapore: World Scientific, 2014).
  * (6) G. Oshanin, R. Metzler, K. Lindenberg (Eds.) Chemical Kinetics: Beyond the Textbook (New Jersey: World Scientific, 2019).
  * (7) D. S. Grebenkov, “NMR Survey of Reflected Brownian Motion”, Rev. Mod. Phys. 79, 1077-1137 (2007).
  * (8) O. Bénichou and R. Voituriez, “From first-passage times of random walks in confinement to geometry-controlled kinetics”, Phys. Rep. 539, 225-284 (2014).
  * (9) D. Holcman and Z. Schuss, “The Narrow Escape Problem”, SIAM Rev. 56, 213-257 (2014).
  * (10) I. V. Grigoriev, Y. A. Makhnovskii, A. M. Berezhkovskii, and V. Y. Zitserman, “Kinetics of escape through a small hole”, J. Chem. Phys. 116, 9574 (2002).
  * (11) A. Singer, Z. Schuss, D. Holcman, and R. S. Eisenberg, “Narrow Escape, Part I”, J. Stat. Phys. 122, 437-463 (2006).
  * (12) A. Singer, Z. Schuss, and D. Holcman, “Narrow Escape, Part II The circular disk”, J. Stat. Phys. 122, 465 (2006).
  * (13) A. Singer, Z. Schuss, and D. Holcman, “Narrow Escape, Part III Riemann surfaces and non-smooth domains”, J. Stat. Phys. 122, 491 (2006).
  * (14) S. Condamin, O. Bénichou, V. Tejedor, R. Voituriez, and J. Klafter, “First-passage time in complex scale-invariant media”, Nature 450, 77 (2007).
  * (15) O. Bénichou and R. Voituriez, “Narrow-Escape Time Problem: Time Needed for a Particle to Exit a Confining Domain through a Small Window”, Phys. Rev. Lett. 100, 168105 (2008).
  * (16) O. Bénichou, D. S. Grebenkov, P. Levitz, C. Loverdo, and R. Voituriez, “Optimal Reaction Time for Surface-Mediated Diffusion”, Phys. Rev. Lett. 105, 150606 (2010).
  * (17) O. Bénichou, C. Chevalier, J. Klafter, B. Meyer, and R. Voituriez, “Geometry-controlled kinetics”, Nature Chem. 2, 472-477 (2010).
  * (18) S. Pillay, M. J. Ward, A. Peirce, and T. Kolokolnikov, “An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part I: Two-Dimensional Domains”, SIAM Multi. Model. Simul. 8, 803-835 (2010).
  * (19) A. F. Cheviakov, M. J. Ward, and R. Straube, “An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part II: The Sphere”, SIAM Multi. Model. Simul. 8, 836-870 (2010).
  * (20) D. S. Grebenkov, “Searching for partially reactive sites: Analytical results for spherical targets”, J. Chem. Phys. 132, 034104 (2010).
  * (21) A. F. Cheviakov, A. S. Reimer, and M. J. Ward, “Mathematical modeling and numerical computation of narrow escape problems”, Phys. Rev. E 85, 021131 (2012).
  * (22) C. Caginalp and X. Chen, “Analytical and Numerical Results for an Escape Problem”, Arch. Rational. Mech. Anal. 203, 329-342 (2012).
  * (23) T. G. Mattos, C. Mejia-Monasterio, R. Metzler, and G. Oshanin, “First passages in bounded domains: When is the mean first passage time meaningful”, Phys. Rev. E 86, 031143 (2012).
  * (24) A. M. Berezhkovsky and L. Dagdug, “Effect of Binding on Escape from Cavity through Narrow Tunnel”, J. Chem. Phys. 136, 124110 (2012).
  * (25) J.-F. Rupprecht, O. Benichou, D. S. Grebenkov, and R. Voituriez, “Exit time distribution in spherically symmetric two-dimensional domains”, J. Stat. Phys. 158, 192-230 (2015).
  * (26) A. Godec and R. Metzler, “First passage time distribution in heterogeneity controlled kinetics: going beyond the mean first passage time”, Sci. Rep. 6, 20349 (2016).
  * (27) A. Godec and R. Metzler, “Universal Proximity Effect in Target Search Kinetics in the Few-Encounter Limit”, Phys. Rev. X 6, 041037 (2016).
  * (28) D. S. Grebenkov, “Universal formula for the mean first passage time in planar domains”, Phys. Rev. Lett. 117, 260201 (2016).
  * (29) J. S. Marshall, “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 (2016).
  * (30) D. S. Grebenkov and G. Oshanin, “Diffusive escape through a narrow opening: new insights into a classic problem”, Phys. Chem. Chem. Phys. 19, 2723-2739 (2017).
  * (31) Y. Lanoiselée, N. Moutal, and D. S. Grebenkov, “Diffusion-limited reactions in dynamic heterogeneous media”, Nature Commun. 9, 4398 (2018).
  * (32) D. S. Grebenkov, R. Metzler, and G. Oshanin, “Towards a full quantitative description of single-molecule reaction kinetics in biological cells”, Phys. Chem. Chem. Phys. 20, 16393-16401 (2018).
  * (33) D. S. Grebenkov, R. Metzler, and G. Oshanin, “Strong defocusing of molecular reaction times results from an interplay of geometry and reaction control”, Commun. Chem. 1, 96 (2018).
  * (34) V. Sposini, A. Chechkin, and R. Metzler, “First passage statistics for diffusing diffusivity”, J. Phys. A: Math. Theor. 52, 04LT01 (2019).
  * (35) D. S. Grebenkov, R. Metzler, and G. Oshanin, “Full distribution of first exit times in the narrow escape problem”, New J. Phys. 21, 122001 (2019).
  * (36) N. Levernier, M. Dolgushev, O. Bénichou, R. Voituriez, and T. Guérin, “Survival probability of stochastic processes beyond persistence exponents”, Nature Comm. 10, 2990 (2019).
  * (37) D. S. Grebenkov, “Spectral theory of imperfect diffusion-controlled reactions on heterogeneous catalytic surfaces”, J. Chem. Phys. 151, 104108 (2019).
  * (38) D. Hartich, and A. Godec, “Extreme value statistics of ergodic Markov processes from first passage times in the large deviation limit”, J. Phys. A: Math. Theor. 52, 244001 (2019).
  * (39) D. Hartich, and A. Godec, “Interlacing relaxation and first-passage phenomena in reversible discrete and continuous space Markovian dynamics”, J. Stat. Mech. 024002 (2019).
  * (40) D. S. Grebenkov, “Paradigm Shift in Diffusion-Mediated Surface Phenomena”, Phys. Rev. Lett. 125, 078102 (2020).
  * (41) D. S. Grebenkov, “Diffusion toward non-overlapping partially reactive spherical traps: fresh insights onto classic problems”, J. Chem. Phys. 152, 244108 (2020).
  * (42) S. N. Majumdar, A. Pal, and G. Schehr, “Extreme value statistics of correlated random variables: a pedagogical review”, Phys. Rep. 840, 1-32 (2020).
  * (43) G. H. Weiss, K. E. Shuler, and K. Lindenberg, “Order Statistics for First Passage Times in Diffusion Processes”, J. Stat. Phys. 31, 255-278 (1983).
  * (44) K. Basnayake, A. Hubl, Z. Schuss, and D. Holcman, “Extreme narrow escape: Shortest paths for the first particles among n to reach a target window”, Phys. Lett. A 382, 3449-3454 (2018).
  * (45) K. Basnayake, Z. Schuss, and D. Holcman, “Asymptotic formulas for extreme statistics of escape times in 1, 2 and 3-dimensions”, J. Nonlinear Sci. 29, 461-499 (2019).
  * (46) D. S. Lawley and J. B. Madrid, “A Probabilistic Approach to Extreme Statistics of Brownian Escape Times in Dimensions 1, 2, and 3,” J. Nonlinear Sci. 30, 1207-1227 (2020).
  * (47) S. D. Lawley, “Distribution of extreme first passage times of diffusion”, J. Math. Biol. 80, 2301-2325 (2020).
  * (48) D. S. Grebenkov, R. Metzler, and G. Oshanin, “From single-particle stochastic kinetics to macroscopic reaction rates: fastest first-passage time of N random walkers”, New J. Phys. 22, 103004 (2020).
  * (49) K. Reynaud, Z. Schuss, N. Rouach, and D. Holcman, “Why so many sperm cells”, Commun. Integr. Biol. 8, e1017156 (2015).
  * (50) Z. Schuss, K. Basnayake, and D. Holcman, “Redundancy principle and the role of extreme statistics in molecular and cellular biology”, Phys. Life Rev. 28, 52-79 (2019).
  * (51) J. Madrid and S. D. Lawley, “Competition between slow and fast regimes for extreme first passage times of diffusion”, J. Phys. A: Math. Theor. 53, 335002 (2020).
  * (52) M. Reva, D. A. DiGregorio, and D. S. Grebenkov, “A first-passage approach to diffusion-influenced reversible binding: insights into nanoscale signaling at the presynapse”, Sci. Rep. 11, 5377 (2021).
  * (53) D. S. Grebenkov, “First passage times for multiple particles with reversible target-binding kinetics”, J. Chem. Phys. 147, 134112 (2017).
  * (54) S. D. Lawley and J. B. Madrid, “First passage time distribution of multiple impatient particles with reversible binding”, J. Chem. Phys. 150, 214113 (2019).
  * (55) M. J. Berridge, M. D. Bootman, and H. L. Roderick, “Calcium: calcium signalling: dynamics, homeostasis and remodelling”, Nat. Rev. Mol. Cell Biol. 4, 517 (2003).
  * (56) E. Eggermann, I. Bucurenciu, S. P. Goswami, and P. Jonas, “Nanodomain coupling between Ca2+ channels and sensors of exocytosis at fast mammalian synapses”, Nat. Rev. Neurosci. 13, 7-21 (2012).
  * (57) M. Dittrich et al. “An excess-calcium-binding-site model predicts neurotransmitter release at the neuromuscular junction”, Biophys. J. 104, 2751-2763 (2013).
  * (58) Y. Nakamura et al. “Nanoscale distribution of presynaptic Ca2+ channels and its impact on vesicular release during development”, Neuron 85, 145-158 (2015).
  * (59) C. Guerrier and D Holcman, “Hybrid Markov-mass action law model for cell activation by rare binding events: application to calcium induced vesicular release at neuronal synapses”, Sci. Rep. 6, 1-10 (2016).
  * (60) F. C. Collins and G. E. Kimball, “Diffusion-controlled reaction rates”, J. Coll. Sci. 4, 425 (1949).
  * (61) H. Sano and M. Tachiya, “Partially diffusion-controlled recombination”, J. Chem. Phys. 71, 1276-1282 (1979).
  * (62) D. Shoup and A. Szabo, “Role of diffusion in ligand binding to macromolecules and cell-bound receptors”, Biophys. J. 40, 33 (1982).
  * (63) R. Zwanzig, “Diffusion-controlled ligand binding to spheres partially covered by receptors: an effective medium treatment”, Proc. Natl. Acad. Sci. USA 87, 5856 (1990).
  * (64) B. Sapoval, “General Formulation of Laplacian Transfer Across Irregular Surfaces”, Phys. Rev. Lett. 73, 3314-3317 (1994).
  * (65) M. Filoche and B. Sapoval, “Can One Hear the Shape of an Electrode? II. Theoretical Study of the Laplacian Transfer”, Eur. Phys. J. B 9, 755-763 (1999).
  * (66) O. Bénichou, M. Moreau, and G. Oshanin, “Kinetics of stochastically gated diffusion-limited reactions and geometry of random walk trajectories”, Phys. Rev. E 61, 3388-3406 (2000).
  * (67) D. S. Grebenkov, M. Filoche, and B. Sapoval, “Spectral Properties of the Brownian Self-Transport Operator”, Eur. Phys. J. B 36, 221-231 (2003).
  * (68) A. Berezhkovskii, Y. Makhnovskii, M. Monine, V. Zitserman, and S. Shvartsman, “Boundary homogenization for trapping by patchy surfaces”, J. Chem. Phys. 121, 11390 (2004).
  * (69) D. S. Grebenkov, Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena, in “Focus on Probability Theory”, Ed. L. R. Velle, pp. 135-169 (Nova Science Publishers, 2006).
  * (70) D. S. Grebenkov, M. Filoche, and B. Sapoval, “Mathematical Basis for a General Theory of Laplacian Transport towards Irregular Interfaces”, Phys. Rev. E 73, 021103 (2006).
  * (71) J. Reingruber and D. Holcman, “Gated Narrow Escape Time for Molecular Signaling”, Phys. Rev. Lett. 103 148102 (2009).
  * (72) S. Lawley and J. P. Keener, “A New Derivation of Robin Boundary Conditions through Homogenization of a Stochastically Switching Boundary”, SIAM J. Appl. Dyn. Sys. 14, 1845-1867 (2015).
  * (73) A. Bernoff, A. Lindsay, and D. Schmidt “Boundary Homogenization and Capture Time Distributions of Semipermeable Membranes with Periodic Patterns of Reactive Sites”, Multiscale Model. Simul. 16, 1411-1447 (2018).
  * (74) D. S. Grebenkov and B.-T. Nguyen, “Geometrical structure of Laplacian eigenfunctions”, SIAM Rev. 55, 601-667 (2013).
  * (75) C. W. Gardiner, Handbook of stochastic methods for physics, chemistry and the natural sciences (Berlin, Springer, 1985).
  * (76) H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, 2nd Ed. (Oxford University Press, 1959).