Reversible reactions controlled by surface diffusion on a sphere
Denis S. Grebenkov





 denis.grebenkov@polytechnique.edu

Laboratoire de Physique de la Matière Condensée (UMR 7643), 
 
CNRS – Ecole Polytechnique, IP Paris, 91128 Palaiseau, France







We study diffusion of particles on the surface of a sphere toward a
partially reactive circular target with partly reversible binding
kinetics.  We solve the coupled diffusion-reaction equations and
obtain the exact expressions for the time-dependent concentration of
particles and the total diffusive flux.  Explicit asymptotic formulas
are derived in the small target limit.  This study reveals the strong
effects of reversible binding kinetics onto diffusion-mediated
reactions that may be relevant for many biochemical reactions on cell
membranes.



02.50.-r, 05.40.-a, 02.70.Rr, 05.10.Gg














Reversible reactions controlled by surface diffusion on a sphere
    Denis S. Grebenkov
    December 30, 2023
================================================================





§ INTRODUCTION


Various vital processes occur on the cell membrane, examples ranging
from cell signaling to ionic and molecular transfer <cit.>.  In
some cases, proteins diffuse on the cell membrane and search for
specific receptors to trigger signal transduction via a cascade of
biochemical reactions.  For instance, G-protein-coupled receptors
transmit their signals to the cell interior by interacting with G
proteins and thus mediate the biological effects of many hormones and
neurotransmitters <cit.>.  A comprehensive theory of
these crucial phenomena requires a better understanding of diffusive
processes on two-dimensional surfaces and, in particular, on a
spherical surface.  In this direction, Bloomfield and Prager studied
the impact of rotational diffusion on tail fiber attachment on
bacteriophage T4 <cit.>.  Assuming that the fiber tip
diffuses freely over the surface of the sphere until it reaches a
baseplate site, they computed the first two moments of the
distribution of reaction times.  Chao et al.  investigated
surface diffusion into a trap as a promising passive mechanism of
localization of cell membrane components <cit.> (see also
<cit.>).  They derived the concentration of components by
solving the time-dependent diffusion equation in the presence of a
perfectly absorbing sink.  Linderman and Lauffenburger analyzed
trapping of diffusing receptors into tubules as a mechanism of
intracellular receptor/ligand sorting <cit.>.  Sano and
Tachiya treated a more general case of imperfect sinks <cit.>,
for which they obtained the mean reaction time and discussed an
exponential approximation for the survival probability.  Berg
considered the effect of weak binding of a diffusing particle to the
nonspecific surface around a target on the sphere <cit.>.  The
resulting surface-mediated diffusion, in which three-dimensional
diffusion in the bulk is coupled to two-dimensional diffusion on the
surface, facilitates the target search.  The optimality of such
intermittent processes was further investigated
<cit.>.
Berg and Purcell estimated the reaction rate controlled by bulk
diffusion toward multiple small patches evenly distributed over the
surface of a sphere <cit.>.  This seminal paper and the
associated homogenization concept were further extended by many
authors
<cit.>.
Shoup et al. devised an efficient approximation to compute the
reaction rate from the bulk to an active circular region on the
surface of a sphere, even in the presence of rotational motion
<cit.>.  This approximation was further extended to compute
the mean first-passage time (MFPT) from the bulk <cit.>.
Singer et al.  generalized the small-target asymptotic behavior
of the MFPT to two-dimensional manifolds <cit.>.  In turn,
Coombs et al.  derived the asymptotic behavior of the MFPT for a
diffusing particle confined to the surface of a sphere, in the
presence of many partially absorbing traps of small radii
<cit.>.  Prüstel and Tachiya considered the problem of
reversible binding reactions on several two-dimensional domains,
including the spherical case <cit.>.  By employing
convolution relations between the survival probabilities for
reversible and irreversible reactions, they derived an exponential
approximation for these probabilities.

In this paper, we further extend these former works by considering
surface diffusion of particles on a sphere toward a partially reactive
circular target with partly reversible binding kinetics
(Sec. <ref>).  By generalizing the approach by Chao et
al. <cit.>, we provide in Sec. <ref> the exact
solution for the concentration of diffusing particles in both Laplace
domain and time domain.  In Sec. <ref>, we illustrate
the properties of our solution, reveal the strong effects of
reversible binding kinetics, and discuss implications for biochemical
reactions on cell membranes.  Explicit asymptotic formulas in the
small target limit are also presented.  Section <ref>
concludes the paper with the summary of results and perspectives for
future research.  Technical derivations are reported in Appendices.







§ MODEL



We consider the following model: a single particle A of radius
ρ is fixed at the South pole of the sphere of radius R whereas
point-like particles B diffuse on the surface of that sphere with
the diffusion coefficient D (Fig. <ref>).  In what
follows, we call the particle A a target (or a sink, or a trap).
This model can also describe the situation when the particle A
diffuses with the coefficient coefficient D_A or/and when the
particles B have a finite radius ρ_B.  In this case, D is
replaced by D + D_A while ρ is replaced by ρ + ρ_B.
However, the particles B are still assumed to be independent,
without any mutual interaction (even without the excluded volume
constraint).  In practice, it means that the concentration of
particles B should be rather dilute.




We study a common diffusion-influenced reaction


    [every arrow/.append style=shift left]
     A + B r    AB lλμ⟶ A + C .

Once a particle B comes in contact with the target particle A, it
can either bind with the target, or continue its diffusive motion on
the surface.  The probability of the binding event is determined by
the reactivity κ (in units m/s) in a standard way
<cit.>.  The
latter can also be related to the conventional association rate
constant  (in units m^2/s/mol) as


    k_ on = κ N_A (2πρ) ,

where N_A ≃ 6· 10^23 1/mol is the Avogadro number, and
2πρ is the target perimeter <cit.>
(another notation, k_d = k_ on/N_A = 2πρκ, was
employed in <cit.>).  A finite reactivity κ
mimics conformational or energy activation barrier that the particles
A and B have to overcome to form a metastable complex AB.  Two
rates λ and μ allow us to distinguish three scenarios of
irreversible (λ = 0), reversible (λ >0, μ = 0),
and partly reversible (λ > 0, μ > 0) binding.  When the
dissociation (or desorption) rate λ is zero, the particle B
never resumes its diffusion, being either chemically transformed into
another particle C with the rate μ, or staying bound forever to
the target (if μ = 0).  In the former case, the target A plays
the role of a catalyst for the transformation of B into C.  When
λ > 0 and μ = 0, the formation of the complex AB is
fully reversible, i.e., the particle B remains bound for a random
exponentially distributed time, determined by the dissociation rate
λ, and then resumes its surface diffusion until the next
binding.  This mechanism accounts for reversible reactions 
<cit.>.  In turn, the last case
λ > 0 and μ > 0 describes partly reversible scenario when
the metastable complex AB can dissociate either into A + B with
the rate λ (in which case B resumes its diffusion until the
next binding), or into A + C with the rate μ.  Note that the
rate λ was also denoted as k_ off or k_a.

In spherical angular coordinates (θ,φ), the diffusion
domain on the surface of the sphere is Ω = { (θ,φ)
 :  0 ≤θ < π - ,  0 ≤φ < 2π}, whereas the
target surface is = { (θ,φ)  : θ = π - ,
 0 ≤φ < 2π}, where  is the angular size of the
target: ρ = R sin.  We introduce the surface concentration
of particles B in a point (θ,φ), c(θ,φ,t)
(in mol/m^2) that obeys the diffusion equation


    ∂_t c = D Δ c    in Ω,

where Δ is the Laplace-Beltrami operator on the sphere:


    Δ = 1/R^2(1/sinθ∂_θsinθ∂_θ 
    + 1/sin^2θ∂^2_φ).

Even though a Dirac point-like source would yield a more general
description, we restrict the analysis to the uniform initial
concentration of particles B:


    c(θ,φ,t=0) = C_0.

This assumption, which is rather common for applications, simplifies
technical computations.  In fact, the axial symmetry of the domain and
of the initial condition ensures that the solution
c(θ,φ,t) does not depend on the azimuthal angle φ
that we drop in the following.  It is worth noting that the developed
method can be extended to a general setting without axial symmetry.

The diffusion equation (<ref>) is completed by the Robin
boundary condition on the target boundary  


    - D (∂_n c)_| = κ  c_| - λn_T/2πρ ,

where ∂_n = (1/R) ∂_θ is the normal derivative
directed toward the target, and n_T(t) (in mol) is the amount of
particles B bound to the target at time t.  Here, n_T/(2πρ) can be understood as the density of bound particles along the
target boundary of length 2πρ.  In chemical physics, the
relation (<ref>) is called “back-reaction”, “generalized
radiation”, or “generalized Collins-Kimball” boundary condition
<cit.>.  It states that the net
diffusive flux density of particles B toward the target (the
left-hand side) is equal to the reactive flux density, κ
c_|, minus the flux density related to the dissociation of the
metastable complex AB (the last term).  The reactivity κ
characterizes the difficulty for a particle B to overcome an energy
activation barrier for binding to the target
<cit.> (see also a recent overview in
<cit.>).  The limit κ = ∞
describes an immediate reaction upon the first encounter (no
activation barrier), whereas κ = 0 corresponds to no reaction
(infinite activation barrier) on reflecting boundary.  The
concentration is also required to be regular (not diverging) at the
North pole that will serve as the second boundary condition at θ
= 0, see below.

In turn, the amount of particles B bound to the target, n_T(t),
obeys the first-order differential equation:


    ∂_t n_T = J(t) - μ  n_T,

with the initial condition n_T(t=0) = 0, where


    J(t) = ∫_ dl  (-D (∂_n c)_|) = -2πρ D (∂_n c)_|

is the total diffusive flux onto the target, and the second term
accounts for irreversible transformation of the metastable complex
AB into A and C.  Once J(t) is found, one gets immediately


    n_T(t) = ∫_0^t dt'   e^-μ (t-t')  J(t').

In the case μ > 0, the amount n_C(t) of produced particles C
obeys an ordinary differential equation,


    ∂_t n_C(t) = μ  n_T(t)  ,

with the initial condition n_C(0) = 0.  Once n_T(t) is found,
n_C(t) follows immediately by integration.








§ SOLUTION



In this section, we solve the coupled partial differential equations
introduced in Sec. <ref>.  As it is common for diffusion
problems, we first obtain the solution in Laplace domain (for
Laplace-transformed quantities) and then proceed to the solution in
time domain via Laplace transform inversion.





 §.§ Solution in Laplace domain


Setting


    x = cosθ,

the Laplace-Beltrami operator can be written as


    Δ = Ł/R^2,     Ł = ∂_x (1-x^2)∂_x.

In order to get rid of the derivative with respect to time t, the
Laplace transform is performed,



    c̃(x,s)    =   ∫_0^∞ dt   e^-ts  c(x,t) , 
    ñ_T(s)    =   ∫_0^∞ dt   e^-ts  n_T(t)


(throughout the paper, tilde denotes Laplace-transformed quantities).
The new functions c̃(x,s) and ñ_T(s) satisfy



    s  c̃(x,s) - D/R^2 Ł c̃(x,s)    =    C_0 , 
    
    (s + λ + μ)  ñ_T(s)/2πρ   =   κ c̃(a,s) , 
    D/R√(1-a^2)(∂_x c̃(x,s))_x=a = κ c̃(a,s)    -   λñ_T(s)/2πρ ,


where we used ∂_n = - (1/R) sinθ ∂_x, the initial
condition n_T(t=0) = 0, and


    a = cos (π-) = - cos = - √(1 - (ρ/R)^2)

is the location of the target surface .  As Eq. (<ref>)
relates ñ_T(s) and c̃(a,s), one can rewrite
Eq. (<ref>) as Robin boundary condition


    (∂_x c̃(x,s))_x=a = q_s  c̃(a,s) ,

with


    q_s = q  s+μ/s+λ + μ ,

where


    q = /√(1-a^2)≥ 0  ,      = κ R/D

(throughout the paper, bar denotes dimensionless quantities rescaled
by R and D).  In other words, the above coupled differential
equations for c and n_T are decoupled in Laplace domain and
reduced to a single differential equation


    (R^2 s/D - Ł) c̃ = R^2 C_0/D     (a < x < 1) ,

subject to the Robin boundary condition (<ref>) at x = a
(the target), and the regularity condition at x = 1 (the North
pole).  This implies that the reversible kinetics can be described in
Laplace domain as irreversible one, with a “frequency”-dependent
reactivity


    κ(s) = κ s+μ/s+λ+μ .

In time domain, this would result in a convolution-type Robin boundary
condition with an “effective” time-dependent reactivity.  We note
that Eqs. (<ref>, <ref>) present an extension of the
well-known back-reaction boundary condition in Laplace domain (see,
e.g., <cit.>) by accounting for the rate μ.


The solution of the above equation can be obtained in terms of the
Legendre functions of the first kind, P_ν(x), that satisfy the
Legendre equation on a < x < 1,


    ∂_x (1-x^2) ∂_x_Ł P_ν(x) + ν(ν+1) P_ν(x) = 0  ,

and are regular at x = 1 (see Appendix <ref> for some
properties of these functions).  The Robin boundary condition at x =
a reads


    P'_ν(a) = q_s   P_ν(a) ,

where the prime denotes the derivative with respect to x:
P'_ν(a) = (∂_x P_ν(x))_x=a.  As a and q_s are
fixed, this is an equation on the degree ν that has infinitely
many real-valued solutions, denoted as ν_n^s, with n =
0,1,2,… (see Appendix <ref>):


    0 ≤ν_0^s ≤ν_1^s ≤… ↗ +∞

(negative solutions are not considered because of the symmetry:
P_ν(x) = P_-ν-1(x)).  The eigenvalues of the operator -Ł
are then ν_n^s (ν_n^s+1).  We emphasize that ν_n^s and thus
the eigenfunctions P_ν_n^s(x) depend on s through the
parameter q_s in the Robin boundary condition (<ref>).

As the operator Ł is self-adjoint, its eigenfunctions
P_ν_n^s(x) form a complete basis in the space of
square-integrable functions on (a,1).  As a consequence, one can
search for a solution of Eq. (<ref>) in the form


    c̃(x,s) = ∑_n=0^∞ a_n(s)   b_n   P_ν_n^s(x) ,

where the constants b_n ensure the L_2-normalization of
P_ν_n^s(x),


    b_n = (∫_a^1 dx   [P_ν_n^s(x)]^2 )^-1/2

(see the identity (<ref>) for computing these constants).
The unknown coefficients a_n(s) can be found by substituting the
anzats (<ref>) into Eq. (<ref>), multiplying
by P_ν_m^s(x), integrating from a to 1, and using the
orthogonality of P_ν_n^s(x) (see Appendix
<ref>):


    a_n(s) = C_0   b_n/s + D ν_n^s(ν_n^s+1)/R^2∫_a^1 dx   P_ν_n^s(x)  .
 
The integral of the Legendre function can be easily calculated via
Eq. (<ref>), yielding


    c̃(x,s) = C_0 (1-a^2) ∑_n=0^∞b_n^2   P_ν_n^s(x)   P'_ν_n^s(a)/ν_n^s(ν_n^s+1)(s + Dν_n^s(ν_n^s+1)/R^2)  .

While we used here the uniform initial condition for c(x,t), one can
substitute it by any given (axially symmetric) initial condition,
including the Dirac distribution δ(x-x_0), in which case the
integral of P_ν_n^s(x) in Eq. (<ref>) would simply yield
P_ν_n^s(x_0).

From Eqs. (<ref>, <ref>), one deduces the Laplace
transform of the total diffusive flux onto the target:


    J̃(s) = 2π C_0 D (1 - a^2)^2  ∑_n=0^∞b_n^2   [P'_ν_n^s(a)]^2/ν_n^s(ν_n^s+1)(s + D/R^2ν_n^s(ν_n^s+1)) .

Eqs. (<ref>, <ref>) constitute the exact solution
of the problem in Laplace domain.










 §.§ Solution in time domain


The Laplace transform inversion can be performed with the aid of the
residue theorem that requires finding the poles of
Eq. (<ref>) in the complex plane of s.  We recall that
ν_n^s and thus b_n and P_ν_n^s(x) depend on s through
the parameter q_s in the Robin boundary condition
(<ref>).  Since we deal with a diffusion problem, the poles
should lie exclusively on the negative real axis.  We expect that
P_ν_n^s(x) as a function of s does not have poles (see a short
comment in Appendix <ref>, even though a rigorous
proof of this claim is beyond the scope of this paper) so that the
poles come from the zeros of the denominator in Eq. (<ref>).
In the particular case λ > 0 and μ = 0, the solution
reaches a nontrivial steady-state limit with ν = 0 which is
determined by the pole s = 0 that we will treat separately (in other
cases, this pole does not appear).  The other poles are strictly
negative and determined as the zeros of the functions


    f_n(s) ≡ s + D/R^2ν_n^s(ν_n^s + 1) = 0.

Since ν_n^s as solutions of Eq. (<ref>) depend on s,
one has to search for a pair (s,ν) as a solution of the system of
two nonlinear equations (<ref>, <ref>).
Expressing s from Eq. (<ref>) in terms of ν_n^s, one
reduces this system to a single equation


    ( +  - ν(ν+1)) P'_ν(a) = q( - ν(ν+1)) P_ν(a) ,

with dimensionless rates


    = λ R^2/D,         = μ R^2/D ,

and q is given by Eq. (<ref>).  In Appendix
<ref>, we prove that all the solutions of this
equation are real; moreover, negative solutions are ignored due to the
symmetry P_ν(x) = P_-ν-1(x).  We denote the positive solutions
as


    0 ≤ν_0 ≤ν_1 ≤… ↗ +∞

(without the superscript s) to distinguish them from ν_n^s.  The
poles are simply


    s_n = - D ν_n(ν_n+1)/R^2.

As ν_n obeys the boundary condition (<ref>) for such
s that satisfies Eq. (<ref>), one can formally write ν_n
= ν_n^s_n.

Quite unusually, the numerical computation in time domain turns out to
be simpler than that in the Laplace domain, because the set of
solutions of Eq. (<ref>) has to be computed only once,
whereas solutions of Eq. (<ref>) had to be found for each
value of s.

If all poles s_n are simple, then the inverse Laplace transform of
the concentration reads


    c(x,t)     = C_0 ∑_n=0^∞ c_n   e^- Dt ν_n(ν_n+1)/R^2  P_ν_n(x) , 
    
    c_n     = (1-a^2) b_n^2    P'_ν_n(a)/ν_n(ν_n+1) (∂_s f_n)(s_n)  ,

and that of the total flux is


    J(t)     = C_0 D ∑_n=0^∞ J_n   e^- Dtν_n(ν_n+1)/R^2 ,  
    
    J_n     = 2π (1-a^2)^2  b_n^2   [P'_ν_n(a)]^2 /ν_n(ν_n+1) (∂_s f_n)(s_n) .

In Appendix <ref>, we show how the contribution
(∂_s f_n)(s_n) to each residue can be computed in practice,
whereas Appendix <ref> presents an example of such
computation.

Finally, we also get from Eqs. (<ref>, <ref>,
<ref>):


    n_T(t) = C_0 D ∑_n=0^∞ J_n  e^-μ t - e^-Dtν_n(ν_n+1)/R^2/Dν_n(ν_n+1)/R^2 - μ

(note that if μ coincides with Dν_n(ν_n+1)/R^2 for some n,
the corresponding term would be J_n   t  e^-μ t; similarly, a
specific term appears in the case when μ = 0 and ν_0 = 0) and


    n_C(t)     = μ C_0 D ∑_n=0^∞J_n/Dν_n(ν_n+1)/R^2 - μ
       ×(1-e^-μ t/μ - 1-e^-Dtν_n(ν_n+1)/R^2/Dν_n(ν_n+1)/R^2)

(again, a specific term appears if μ coincides with
Dν_n(ν_n+1)/R^2 for some n).  With these expressions, we
completed obtaining the exact solution of the our model of
diffusion-influenced reactions on the sphere.








§ DISCUSSION



In this section, we discuss the main features of the derived
solutions.




 §.§ Conservation of particles


When μ = 0, the total number of particles B is preserved.  In
fact, integrating Eq. (<ref>) over the domain Ω, we
get an equation for the total number of particles diffusing on the
sphere:


    ∂_t N(t) = D ∫_Ω ds  Δ c =  ∫_ dl   (D ∂_n c)_| = - ∂_t n_T(t),

with the last equality coming from Eq. (<ref>).
The total number of particles, N(t) + n_T(t), is thus preserved:


    N(t) + n_T(t) = N(0) + n_T(0) = 2π (1-a) R^2 C_0 ,

where we used the initial conditions c(x,t=0) = C_0 and n_T(t=0) =
0.

In the steady-state limit (as t→∞), the populations of
particles on the surface and on the target are equilibrated.  As the
steady-state flux is zero, Eq. (<ref>) implies


    c_∞ = c(x,t=∞) = λ/2πρκ  n_T(t=∞) ,

from which Eq. (<ref>) yields


    c_∞ = C_0 (1 + κ/λ R√(1+a/1-a))^-1 .

This expression results from the detailed balance between free and
bound particles B, which is controlled by the reactivity κ
and the dissociation rate λ.








 §.§ Survival probabilities


While the above discussion dealt with the concentration c(x,t) of
particles on the sphere, it is instructive to provide its
probabilistic interpretation.  Given that the initial concentration
was set to be a constant C_0, the dimensionless ratio c(x,t)/C_0
can be interpreted as the probability of finding a particle B, that
started from a point x = cosθ at time 0, in the free state
(not bound to the target) at a later time t.  For irreversible
reactions, it is called the survival probability <cit.> because
irreversible binding to the target can be understood as “killing” of
a free particle B.  The same term was adopted by Agmon et al.
for reversible reactions in <cit.>, in which
c(x,t)/C_0 was denoted as S_ rev(t|x).

Moreover, Agmon et al. also considered the probability S_
rev(t|*) for a particle, that was initially in the bound state, to
be unbound in a later time t.  This probability was expressed in
terms of S_ rev(t|a) (see <cit.>, Eq. (3.15)).
Comparing that relation with Eq. (<ref>) at μ = 0, one
can identify


    S_ rev(t|*) = λ/2πρκ C_0  n_T(t) .

The explicit solution (<ref>) yields thus


    S_ rev(t|*) = λ R^2/2πρκ∑_n=0^∞ J_n  1 - e^-Dtν_n(ν_n+1)/R^2/ν_n(ν_n+1) .

In the long-time limit, this probability approaches a steady-state
limit, which can be either obtained from the above spectral
decomposition, or from Eqs. (<ref>, <ref>,
<ref>):


    S_ rev(t=∞|*) = (1 + κ/λ R√(1+a/1-a))^-1 .

This probability is smaller than 1, in contrast to the case of
diffusion in unbounded domains outside a disk or a sphere, studied in
<cit.>, for which S_
rev(t=∞|*) = 1.

We also note that Eq. (<ref>) with μ = 0 implies that the
time derivative of S_ rev(t|*) gives back the total flux
J(t), so that one recovers the relation (3.17) from <cit.>
for the time-dependent rate coefficient:


    k_ rev(t) = J(t)/C_0 = ∂_t n_T(t)/C_0 = 2πρκ/λ ∂_t S_ rev(t|*).

In the following, we focus on the concentration c(x,t) and the total
flux J(t), bearing in mind that other important quantities can be
directly accessed. 






 §.§ Limiting cases



Our solution generalizes previous works on diffusion-influenced
reactions on the spherical surface.  We discuss several limiting cases
which yield simpler formulas.






  §.§.§ Irreversible reactions


For irreversible reactions (λ = 0), a particle B that is
bound to the target, is never released.  As a consequence, the
boundary condition (<ref>) for the concentration c is
decoupled from n_T, so that one can solve first the boundary value
problem for c with Robin boundary condition, evaluate the diffusive
flux J(t), and then determine n_T from Eq. (<ref>).  In
fact, Eq. (<ref>) is reduced to q_s = q, so that Robin
boundary condition (<ref>) reads


    P'_ν(a) = q P_ν(a).

Its solutions ν_n^s do not depend on s and actually coincide
with the solutions of ν_n of Eq. (<ref>).  As a
consequence, 


    (∂_s f_n)(s_n) = 1,

and Eqs. (<ref>, <ref>) are simplified.  In the following,
we refer to the corresponding concentration and total flux as
c_0(x,t) and J_0(t), where subscript 0 highlights that this
solution corresponds to λ = 0.  This is a generalization of
the solution presented in Ref. <cit.> to the case of a
partially absorbing sink.

Whenever a < 1 and κ > 0 (i.e., a nontrivial partially
reactive target is present), it is easy to check that the first
solution ν_0 of Eq. (<ref>) is strictly positive (in
fact, if ν_0 was equal to 0, then P_ν_0(x) = 1 and thus
Eq. (<ref>) cannot be satisfied).  As a consequence, the
concentration of particles c_0(x,t) and the flux J_0(t) vanish as
t→∞.





  §.§.§ Perfect sink


For a perfect sink (κ = ∞), the Robin boundary condition
(<ref>) is reduced to the Dirichlet condition (also known
as totally absorbing or Smoluchowski boundary condition) c|_ = 0
that reads as


    P_ν(a) = 0 ,

which does not depend on s.  As a consequence, Eq. (<ref>)
holds again, and the concentration c(x,t) and the diffusive flux
J(t) are still determined by Eqs. (<ref>, <ref>), in
which ν_n are the solutions of Eq. (<ref>).  This is
precisely the solution presented in Ref. <cit.>.  Note that
when the sink is perfect, it does not matter whether the reaction is
reversible or not: any released particle will be immediately
re-adsorbed by the target.





  §.§.§ Steady-state regime


For reversible reactions with μ = 0, a steady-state solution is
established at long times so that there should exist an eigenmode with
ν = 0.  First, using the properties of Legendre functions, we
rewrite (<ref>) as


    [a (ν+1) - (1-a^2) q_s ] P_ν(a) = (ν+1) P_ν+1(a) .

Searching for ν_0 in the form η s^α (with some constant
η and degree α) and substituting it into
Eq. (<ref>), we get with the aid of
Eqs. (<ref>)


    ν_0 ≃η  s + O(s^2),    with  η = κ R/Dλ√(1+a/1-a) ,

from which P_ν_0(x) ≃ 1 + O(s), P'_ν_0(a) = q_s +
O(s^2), b_n^2 ≃ (1-a)^-1, and thus the leading term of
Eq. (<ref>) behaves asymptotically as


    c̃(x,s) ≃C_0/(1 + Dη/R^2)s     (s→ 0) .

Inverting the Laplace transform, we retrieve a uniform steady-state
concentration in Eq. (<ref>).







 §.§ Numerical implementation



The concentration c(x,t) and the total flux J(t) are given via
exact spectral decompositions (<ref>, <ref>).  Since the
solutions ν_n grow up to infinity as n→∞, both
decompositions converge exponentially fast for any fixed t > 0.  As
time t stands in front of ν_n(ν_n+1) in the exponential
function, larger Dt/R^2 yields faster convergence.  In turn, the
effect of other parameters on the convergence speed is more subtle and
requires a systematic analysis that is beyond the scope of this paper.

In practice, the spectral decompositions have to be truncated up to
some order , which can be selected according to the minimal
time at which c(x,t) and J(t) need to be calculated.  The
properties of the Legendre function P_ν(x) summarized in Appendix
<ref> were thoroughly used for efficient numerical
computations.  In particular, P_ν(x) and its derivatives with
respect to x and ν were accurately and rapidly computed via
Eqs. (<ref>, <ref>, <ref>).  Using
bisection or Newton's method, one computes numerically the needed
number of the solutions ν_n of Eq. (<ref>).  These
solutions then determine the poles s_n, the normalization constants
b_n and the contribution (∂_s f_n)(s_n) to each residue
(see Appendix <ref> for details), from which the
coefficients c_n and J_n are computed via Eqs. (<ref>,
<ref>).  Note that if some poles are not simple, their
contribution to the spectral decompositions (<ref>, <ref>)
should be modified according the residue theorem; in practice, the
poles were simple in all considered examples.  In Appendix
<ref>, we illustrate this numerical scheme by computing
the total flux.








 §.§ Comparison with Tachiya's approach



In his seminal work <cit.>, Tachiya proposed an elegant
approach to relate the solutions for reversible and irreversible
reactions (see also <cit.>).  His straightforward
probabilistic arguments couple c(x,t) for reversible reactions (with
λ > 0) to c_0(x,t) for irreversible reactions (with λ
= 0) via convolution relations that can be easily solved in the
Laplace domain.  In our notations, one gets


    c̃(x,s) = 1/s + c̃_0(x,s) - 1/s/1 + λ c̃_0(a,s)   s/(s+μ) ,

where we also included the rate μ.  One can easily check that the
function c̃(x,s) satisfies Eqs. (<ref>,
<ref>) with q_s from Eq. (<ref>), given that
c̃_0(x,s) satisfies Eqs. (<ref>,
<ref>) with q instead of q_s.  This is a significant
simplification because Eq. (<ref>) for c̃_0(x,s)
involves ν_n^s that do not depend on s and thus have to be
computed only once.  However, the simple and instructive relation
(<ref>) in the Laplace domain is not much helpful for the
Laplace inversion.  In contrast, our approach allows one to get the
exact solutions in time domain so that two approaches are
complementary to each other.





We also note that our description with μ = 0 is equivalent to that
given by Prüstel and Tachiya <cit.>.  In fact, the
boundary condition (<ref>) can be re-written in the form of
Eq. (16) from Ref. <cit.>:


    2πρD/Rsin∂_x c_| = k_a c_| - k_d (1 - 2π R^2 ∫_a^1 dx'   c(x',t)) ,

with k_a = 2πρκ, k_d = λ, and


    n_T(t) = 1 - 2π R^2 ∫_a^1 dx'   c(x',t).

The time derivative of the last equation yields, after
simplifications, Eq. (<ref>).

Instead of solving exactly the problem for irreversible kinetics,
Prüstel and Tachiya used an exponential approximation for the
survival probability,


    S_0(x,t) = c_0(x,t)/C_0≈exp(- t/ τ(x)),

where τ(x) is the mean first passage time in the case of
irreversible reaction (with λ = 0):


    τ(x) = ∫_0^∞ dt  S_0(x,t)  ,

which was found in <cit.> (see also Appendix
<ref>):


    τ(x) = R^2/Dln(1+x/1+a) + R/κ√(1-a/1+a)  .

With this approximation, the Laplace inversion of
Eq. (<ref>) yielded an approximate formula for the
concentration c(x,t) which reads in our notations as


    c(x,t)/C_0≃λ/λ + 1/τ(a) 
    + e^-t/τ(x)1/τ(x) - 1/τ(a)/1/τ(x) - 1/τ(a) - λ
           - e^-t(λ + 1/τ(a))[λ/λ + 1/τ(a) 
    + λ/1/τ(x) - 1/τ(a) - λ]  .


Since the validity of this approximation was not tested in
<cit.>, we undertake this study here.  Figure <ref>
illustrates the quality of the exponential approximation
(<ref>) for irreversible reactions.  When the
particle B starts far from a small perfectly absorbing target (x =
0, κ = ∞, Fig. <ref>a), the diffusive search
takes long time that is enough for a particle to “forget” about its
starting point.  In this setting, Eq.  (<ref>)
accurately approximates the survival probability.  In contrast, when
the starting point is close to the target (x = -cos(2) ≈
-0.866), there are high changes of finding it in a short time, and
the exponential approximation (<ref>) fails.  It is
also worth noting that the exponential asymptotic decay of the
survival probability S_0(x,t) is determined by the smallest
eigenvalue of the diffusion operator, which is independent of the
starting point.  In this light, the approximation
(<ref>) can only be valid when τ(x) is close
to T = R^2/(Dν_0(ν_0+1)) and thus independent of x, requiring
the starting point to be far from the target and the target to be
small.  The situation is different for weakly reactive target with
κ R/D = 1 (Fig. <ref>b).  Even when the starting point
is close to the target, multiple failed attempts of binding to the
target allow for the particle to explore the whole domain until the
successful binding.  In other words, the second term in
Eq. (<ref>) is dominant, whereas the dependence on the
starting point x is irrelevant.  In summary, for a fixed x, the
approximation (<ref>) is getting more accurate when the
target size and/or reactivity decrease.  In particular, this
approximation can be employed to re-derive the asymptotic behavior of
the total flux in the small-target limit, as an alternative to the
rigorous analysis presented in Sec. <ref>.  At the
same time, the dependence on the starting point is not captured, and
the approximation fails when the starting point is close to a highly
reactive target, irrespectively of its size.  In this light, an even
simpler approximation S_0(x,t) ≈exp(-t/τ(0)) would be
more natural.






Figure <ref> shows the normalized concentration c(x,t)/C_0
for reversible reactions (with λ > 0 and μ = 0) as a
function of rescaled time Dt/R^2 for the same target and the same
starting points.  As a perfectly reactive target would yield
irreversible reactions, we do not consider the case κ = ∞
but compare two cases with high reactivity κ R/D = 10 and weak
reactivity κ R/D = 1.  To keep the same steady-state
concentration c_∞, we set λ R^2/D = 1 and λ
R^2/D = 0.1, respectively.  As discussed previously for irreversible
reactions, the approximate solution (<ref>) fails when
the starting point is close to the target, even for weakly reactive
target.  In turn, the approximation is getting more accurate when the
starting point is far from the target and the reactivity is low.














 §.§ Small-target asymptotic behavior



In many practical applications, the radius ρ of the target is
much smaller than the radius R of the sphere, i.e., ≪ 1 and
a is close -1.  In the conventional first-passage time framework,
this situation is known as the narrow escape problem
<cit.>.  Here we summarize the asymptotic behavior
of our solutions in the narrow escape limit → 0 (see details
in Appendix <ref>).

In the limit = 0, a point-like target cannot be accessed by
Brownian motion on the sphere and thus corresponds to no target
situation, in which the spectrum of the governing operator Ł is
simply ν_n = n for n = 0,1,2,…, whereas the eigenfunctions
P_ν(x) are the Legendre polynomials.  In particular, the
concentration c(x,t) remains constant, whereas the flux J(t) is
zero.  

When  is small but nonzero, finding such a small target is a long
search process.  In Appendix <ref>, we deduce the long-time
behavior of the total flux as


    J(t) ≃4π C_0 R^2/T  e^-t/T     (t →∞) ,

with a characteristic time T, which is determined by the smallest
eigenvalue of the diffusion operator:


    T = R^2/D ν_0(ν_0+1)≃R^2/Dν_0     (≪ 1).

Since the factor N_B = 4π C_0 R^2 ≃ 2π(1-a)C_0 R^2 is
simply the total initial number of diffusing particles B, the above
expression (<ref>) can be understood as an exponential
distribution of the first-passage times for N_B independent
particles.  The asymptotic behavior of the smallest degree ν_0 was
obtained in Appendix <ref> for four different scenarios:

(i) For a perfectly reactive target (κ = ∞), one gets


    ν_0 ≃1/2ln(2/)     (≪ 1) .

The associated time scale T = 2R^2 ln(2/)/D is close to the mean
FPT to a perfectly reactive target (the first term in
Eq. (<ref>)), see also <cit.>.

(ii) For irreversible reactions (λ = 0, κ < ∞),
the smallest degree ν_0 vanishes linearly with , 


    ν_0 ≃κ R/2D≈κρ/2D     (≪ 1).

The associated time scale T = 2R^2/(κρ) does not depend on
the diffusion coefficient and is close to the mean FPT to a weakly
reactive target (the second term in Eq. (<ref>)), see
also <cit.>.

(iii) For reversible reactions (λ > 0, μ = 0), the
smallest degree ν_0 is equal to 0 for any target size.  In this
case, the steady-state total flux is zero, J(∞) = 0, while the
steady-state concentration remains close to the initial concentration
C_0:


    c_∞≃ C_0 (1 + κρ/2λ)^-1 ,

as expected.  The approach to the steady state is controlled by the
next eigenvalue, ν_1(ν_1+1) = 2 + O(), which sets the time
scale R^2/(2D).  However, if the dissociation rate λ is
small, < 1, the approach is controlled by the solution in
Eq. (<ref>), which sets the longer time scale
1/λ.

(iv) For partly reversible reactions (λ > 0, μ > 0), we
obtain (see Appendix <ref>)


    ν_0 ≃μκρ/2D(μ+λ) ,

which is reduced to the second situation by setting λ = 0 and
to the third situation by setting μ = 0.  This value determines
the time scale T = 2R^2(1 + λ/μ)/(κρ).

In analogy to the exact solution (<ref>) for the mean
FPT, we propose an interpolation formula for the time scale T:


    T = 2R^2 ln(2R/ρ)/D + 2R^2(1 + λ/μ)/κρ .

In the limit κ = ∞, one retrieves the time scale for a
perfectly reactive target.  In turn, for moderate values of κ,
the second term provides the dominant contribution in the small-target
limit (ρ≪ R).


Figure <ref>(a) illustrates an excellent agreement between
numerically computed values of ν_0 and its asymptotic relations.
In addition, Fig. <ref>(b) shows the accuracy of the long-time
approximation (<ref>) for a target of angular size =
0.01.  One can see that our explicit formulas accurately describe the
long-time behavior of the total flux toward the target in different
regimes.










 §.§ Inclusion of constant flux from the bulk



Our description can be further extended e.g., by incorporating a
constant flux of particles from the bulk.  In this case, the
concentration is constantly fed, and the diffusion equation
(<ref>) can be modified by adding the constant flux density
J (in units mol/m^2/s) to the right-hand side.  In the Laplace
domain, this is equivalent to changing C_0 to C_0 + J/s in
Eq. (<ref>), as well as in the solutions (<ref>,
<ref>).  As a consequence, one gets


    c_ b(x,t) = c(x,t) + J/C_0∫_0^t dt'  c(x,t')

and 


    J_ b(t) = J(t) + J/C_0∫_0^t dt'  J(t'),

where c(x,t) and J(t) are the former solutions without bulk flux.
Substitution of the exact decompositions (<ref>,
<ref>) into Eqs. (<ref>, <ref>) yields


    c_ b(x,t) = c_ b(x,∞)     + C_0∑_n=0^∞[1 - J R^2/DC_0ν_n(ν_n+1)] 
       × c_n   P_ν_n(x)   e^-Dt ν_n(ν_n+1)/R^2

and


    J_ b(t) = J_ b(∞) + C_0 D∑_n=0^∞   [1 - J R^2/DC_0ν_n(ν_n+1)] 
       × J_n   e^-Dt ν_n(ν_n+1)/R^2 ,

where 


    c_ b(x,∞)     = J R^2/D∑_n=0^∞c_n/ν_n(ν_n+1)  P_ν_n(x) , 
    
    J_ b(∞)     = J R^2 ∑_n=0^∞J_n/ν_n(ν_n+1)

are the steady-state concentration and diffusive flux, respectively.
In contrast to the case without bulk flux, the concentration and the
diffusive flux do not vanish as t→∞ but reach a steady-state
regime.  Expectedly, this flux does not depend on the diffusion
coefficient D, except for the implicit dependence via ν_n.  Note
that the term n = 0 should be treated separately in the case of
fully reversible reactions (λ > 0, μ = 0), for which
ν_0 = 0.











§ CONCLUSIONS AND PERSPECTIVES



In this paper, we investigated diffusion of particles on the surface
of a sphere toward a partially reactive target with partly reversible
binding kinetics.  In this setting, a diffusing particle binds to the
target with a certain probability and then, after some residence time
on the target, it can either unbind to resume diffusion, or be
transformed into another particle.  These reaction mechanisms were
characterized by three parameters: the reactivity κ
(controlling binding to the target), the dissociation rate λ
(controlling unbinding from the target), and the transformation rate
μ (controlling the ultimate irreversible chemical transformation).
This is a rather general reaction model that englobes most of formerly
considered situations.  In particular, we distinguished four
scenarios: (i) an ideal sink (κ = ∞; λ does not
matter), (ii) a partially reactive target with irreversible binding
(κ < ∞ and λ = 0), (iii) a partially reactive
target with reversible binding and no chemical transformation (κ
< ∞, λ > 0, μ = 0), and (iv) a partially reactive
target with reversible binding and chemical transformation (κ <
∞, λ > 0, μ > 0).  Most former studies were
dedicated to the first and (less) to the second situation.

Employing the axial symmetry of the problem, we derived the exact
solution of the underlying diffusion-reaction equations in terms of
Legendre functions.  While the case of irreversible reactions on the
surface of a sphere was studied in the past <cit.>, the
coupling of differential equations for reversible reactions presented
a mathematical obstacle in getting an exact solution in time domain.
We solved this problem and derived a new equation determining the
eigenvalues of the governing diffusion operator.  While our solution
is valid for any target size, further simplifications were achieved in
the small-target limit.  Here, we determined the asymptotic behavior
of the eigenvalues of the governing operator.  In particular, the
smallest eigenvalue controls the long-time exponential decay of the
total flux.  Moreover, we proposed a simple interpolation formula
(<ref>) for the time scale T of this decay that includes
all the reaction parameters.

While surface diffusion on a cell membrane is an emblematic problem
that inspired our work, the proposed method is not limited to cellular
biology.  For instance, one can think of animal or human migration on
Earth, including the problems of epidemic spreading, in which
first-passage time phenomena are critically relevant, even though the
dynamics of such processes is usually more complex than ordinary
diffusion.  More generally, this geometric setting allows one to study
the effect of curvature onto diffusion-reaction processes, as compared
to a flat annulus geometry on the plane.  The latter case is simpler
and was much more thoroughly studied, especially for irreversible
reactions <cit.>.  Following the
same lines as in Sec. <ref>, one can derive the exact
solution describing partly reversible diffusion-influenced reactions
(see also <cit.>, in which such a solution was derived for
the case μ = 0).  A systematic analysis of the curvature effect
onto diffusion-reaction processes presents an interesting perspective
of this work.


Due to a finite surface area of the sphere and thus a finite
amount of diffusing particles, a steady-state solution of the
underlying diffusion equation in the presence of a sink is zero, so
that one has to look for a transient, time-dependent solution.
Moreover, the passage from the bounded surface on the sphere to the
whole plane does not resolve the difficulties <cit.>.  In
fact, the recurrent character of Brownian motion implies that there is
no steady-state solution in the plane with a finite concentration at
infinity.  This is in contrast to three-dimensional unbounded domains,
for which a nonzero steady-state regime can be reached.  In
particular, the reaction rate of particles diffusing in ^3 toward
a spherical sink was first computed a century ago by Smoluchowski
<cit.> and then extended to various target
configurations (see
<cit.> and references
therein).

Our exact solution for the concentration of particles revealed
limitations of some common practices.  In particular, the survival
probability of a particle diffusing in a bounded domain toward a
target is often approximated by a single exponential function (see
<cit.> and references therein).  This is a
long-time approximation, which is accurate when the starting point is
far from the target or when it is averaged over the whole domain.  In
both cases, the particle diffuses for long time until its binding to
the target and thus can largely explore the whole bounded domain,
resulting in such an exponential law.  However, if the particle starts
relatively close to the target, it can bind to the target very
rapidly, on a time scale much shorter than that of the diffusive
exploration of the domain <cit.>.
As a consequence, the exponential function fails to describe the
intricate short-time behavior of the survival probability.  This is
particularly relevant for reversible reactions because the dissociated
particle is released on the target.  For instance, we showed in
Sec. <ref> that the approximation by Prüstel and Tachiya
based on the exponential approximation may be accurate for distant
starting points but fails if the starting point is close to the
target.

The developed description of partly reversible diffusion-influenced
reactions is not limited to a spherical shape.  It can be easily
adapted to general manifolds or Euclidean domains.  While the
possibility of getting the exact solution strongly relied on the
symmetries of the spherical surface, small-target asymptotic results
may potentially be derived without knowledge of an exact solution, as
it was done for ideal sinks <cit.> or partially
absorbing traps on a cylindrical dendritic membrane
<cit.>.  Apart from considering other domains and
manifolds, the presented approach can be extended in different
directions.  First, we showed in Sec. <ref> how to
incorporate a constant flux density from the bulk and to get
stationary solutions.  In this way, surface diffusion toward the
target can be coupled to bulk diffusion (see also <cit.> and
references therein).  The bulk flux may account, e.g., for particles
coming from the exterior of a cell membrane or from the interior of
the cell.  The latter case can model “recycling” of proteins by the
cell.  Second, one can consider diffusion on a part of the sphere
delimited by two circles of latitude.  This is equivalent to solving
the diffusion-reaction equation (<ref>) for c̃(x,s) on an interval (a,b).  As previously, the endpoint a
corresponds to the target, whereas the circle of latitude determined
by cosθ = b can mimic a biological frontier that confines
proteins, e.g., in a lipid raft <cit.>.
Alternatively, if there are many targets that are more or less
uniformly distributed on the surface of a sphere, one can virtually
split this surface into “zones of influence” of each target.  Even
if the particle can cross the virtual border of each zone of
influence, such a crossing would simply mean changing the zone of
influence and can thus be modeled by reflecting boundary condition.  A
solution in such a geometric setting would involve Legendre functions
of two kinds, P_ν(x) and Q_ν(x).  The problem with multiple
small targets may also be addressed by extending the asymptotic
methods developed in <cit.>.  Third, one can investigate
first-passage quantities describing multiple independent particles
such as, e.g., the first moment when n among N particles are bound
to the target.  The reversible binding makes these statistics highly
nontrivial <cit.>.  Finally, an extension of the
proposed formalism beyond Brownian motion presents an important
perspective.  In fact, single-particle tracking experiments on cell
membranes witness considerable deviations from ordinary diffusion
<cit.>.  Some theoretical models of
anomalous diffusion such as continuous-time random walks or processes
with diffusing diffusivity can be implemented into our description, at
least in Laplace domain (see discussions in
<cit.> and references
therein).






The author acknowledges Prof. D. Calebiro for inspiring discussions.










§ SOME PROPERTIES OF LEGENDRE FUNCTIONS



The Legendre function P_ν(x) of the first kind is the solution of
the Legendre equation (<ref>), which is regular at x
= 1 and normalized such that P_ν(1) = 1.  For an integer ν,
P_ν(x) is the Legendre polynomial of degree ν.  In general,
the Legendre function is expressed through the hypergeometric
function:


    P_ν(x) =  _2F_1(-ν,ν+1; 1; (1-x)/2) ,

which is defined as


    _2F_1(a,b;c;z) = ∑_n=0^∞(a)_n   (b)_n/(c)_n  n!  z^n ,

with (a)_n = a(a+1)… (a+n-1) (and (a)_0 = 1).  For
non-integer ν, the Legendre function P_ν(x) has a logarithmic
divergence at x = - 1.  In turn, for any fixed x > -1, P_ν(x)
is an entire analytic function of ν because the hypergeometric
series in Eq. (<ref>) converges uniformly for any finite domain
of ν∈, see <cit.> (p.68) for a similar statement for
the hypergeometric function.

Many properties of the Legendre function are discussed in the classic
textbook by Hobson <cit.>, while their applications in electric
potential problems are summarized in <cit.> (see also
<cit.>).  Here we reproduce some properties that are
relevant for our work.

Multiplying Eq. (<ref>) by P_ν'(x), subtracting its
symmetrized version and integrating over x from a to 1, one
immediately deduces


    (ν-ν')(ν+ν'+1)/1-a^2∫_a^1 dx   P_ν(x)   P_ν'(x) 
           =  P_ν'(a) P'_ν(a) - P_ν(a) P'_ν'(a) ,

where P'_ν(a) = (∂_x P_ν(x))_x=a.  In turn, the
integral of Eq. (<ref>) from a to 1 yields


    ∫_a^1 dx   P_ν(x) = P_ν-1(a) - P_ν+1(a)/2ν+1 =
    1-a^2/ν(ν+1) P'_ν(a),

where we used the recurrence relation


    (ν+1) P_ν+1(x) = (2ν+1) x P_ν(x) - ν P_ν-1(x)

and the expression for the derivative with respect to x:


    P'_ν(x) = (ν+1) P_ν+1(x) - x P_ν(x)/x^2-1 .




The normalization constants b_n in Eq. (<ref>) can be
computed in a standard way.  On one hand, the Legendre equation
(<ref>) is multiplied by ∂_ν P_ν(x) and
integrated over x from a to 1.  On the other hand, the Legendre
equation is differentiated with respect to ν, multiplied by
P_ν(x) and then integrated.  Taking the difference of two
equations and integrating by parts, one gets


    ∫_a^1 dx  [P_ν(x)]^2 = 1/2ν+1
          ×.
    ((1-x^2) [∂_ν P_ν(x)   P'_ν(x) - P_ν(x)  ∂_ν P'_ν(x) ])|_x=a^x=1 .

Taking the derivative of Eq. (<ref>) with respect to
ν, one finds


    ∂_ν P'_ν(x) = P'_ν(x)/ν+1 + ν+1/x^2-1(∂_ν P_ν+1(x) - x ∂_ν P_ν(x))  ,

from which


    ∫_a^1 dx  [P_ν(x)]^2 = 1/2ν+1(P_ν(x)(P_ν+1(x)-xP_ν(x)) 
          + . (ν+1) [P_ν(x) ∂_ν P_ν+1(x) - P_ν+1(x) ∂_ν P_ν(x)
     ]) |_x=a^x=1 .

Since P_ν(1) = 1, the normalization coefficients b_n in
Eq. (<ref>) are obtained from


    ∫_a^1 dx  [P_ν(x)]^2 = -1/2ν+1(P_ν(a)(P_ν+1(a)-aP_ν(a)) 
       + (ν+1) [P_ν(a) ∂_ν P_ν+1(a) - P_ν+1(a) ∂_ν P_ν(a)  ])  .


For Dirichlet boundary condition P_ν(a) = 0, one recovers the
relation from <cit.>:


    ∫_a^1 dx [P_ν(x)]^2 = -1-a^2/2ν+1(∂_ν P_ν(x)   P'_ν(x))_x=a .

Some approximate formulas for computing ∂_ν P_ν(a) were
also provided in that reference.



Hobson derived the following representation of P_ν(x) for
negative x <cit.>:



    P_ν(x) = sin(πν)/π{[ln(1+x/2) + 2γ + ψ(ν+1) + ψ(-ν)] 
    P_ν(-x)
     + ∑_k=1^∞(-1)^k Γ(ν+k+1)/Γ(ν-k+1) (k!)^2 φ(ν,k) (1+x/2)^k } ,


where ψ(ν) = d logΓ(ν)/dν is the digamma
function, γ = 0.5772… is the Euler constant, and


    φ(ν,k)    =     - 2 (1/1 + … + 1/k) +
    (1/ν+1 + … + 1/ν+k)  
       +   (1/-ν + … + 1/-ν+k-1) .

Schelkunoff gave approximate expressions for Legendre functions of
nearly integral degree <cit.>, e.g.,

  

    P_ϵ(x)    ≃    1 + ϵln(1+x/2) , 
    
    P_1+ϵ(x)    ≃    x + ϵ[x-1 + x ln(1+x/2) ],


which are applicable for all |x|<1 and yield an error of the order
of ϵ^2.


To compute Legendre functions numerically at large ν, it is
convenient to use the Mehler-Dirichlet integral representation


    P_ν(cosθ) = √(2)/π∫_0^θ dαcos ((ν+1/2)α)/√(cosα - cosθ) ,

which is valid for θ∈ (0,π) and ν∈
<cit.>.

From this representation, one can also evaluate


    ∂_ν P_ν(cosθ) = - √(2)/π∫_0^θ dααsin ((ν+1/2)α)/√(cosα - cosθ) .

These two relations were used throughout this paper for numerical
computations of P_ν(x) and ∂_ν P_ν(x).

Note that the derivative of the Legendre function with respect to its
degree can also be written as <cit.>:


    ∂_ν P_ν(x) = π(πν) P_ν(x) - 1/π A_ν(x) ,

where


    A_ν(x)    =   sin(πν) ∑_k=0^∞((1-x)/2)^k Γ(k-ν) Γ(k+ν+1)/(k!)^2
       ×   (ψ(k+ν+1) - ψ(k-ν)) .

One also gets


    ∂_ν P'_ν(a) = P'_ν(a)/ν+1 + ν+1/a^2-1[∂_ν P_ν+1(a) - a ∂_ν P_ν(a)] .

However, these relations are less convenient than the integral
presentation (<ref>) for numerical computations.









§ SPECTRAL PROPERTIES




 §.§ Solutions ν_n^s of Eq. (<ref>)



If ν and ν' denote any two solutions of Eq. (<ref>)
(with a fixed q_s), the right-hand side of
Eq. (<ref>) is zero, ensuring the orthogonality of
the corresponding Legendre functions.  This property holds for Neumann
(q_s = 0), Robin (0 < q_s < ∞) and Dirichlet (q_s = ∞)
boundary conditions.

Moreover, MacDonald used Eq. (<ref>) to prove that
all zeros of P_ν(a) are real <cit.> (see also an
amendment to the proof in <cit.>).  His argument can be
directly applied to prove that all zeros of P'_ν(a) - q_s P_ν(a)
are also real for any fixed real q_s.  It is enough to rewrite the
right-hand side of Eq. (<ref>) as


    1/q_s( [q_s P_ν'(a) - P'_ν'(a)] P'_ν(a) - [q_s P_ν(a) - P'_ν(a)] P'_ν'(a) ).

If ν was a complex zero of P'_ν(a) - q_s P_ν(a), then its
complex conjugate, ν'= ν^*, would also be a zero (given that
P_ν^*(a) = P_ν^*(a)), implying that the integral in
Eq. (<ref>) would vanish.  This is not possible
because P_ν(x)   P_ν'(x) = |P_ν(x)|^2 for such a pair of
complex conjugate degrees ν and ν' (see <cit.> for
details).

As P_ν(x) is an entire analytic function of ν for any fixed x
> - 1 (see Appendix <ref>), whereas ν^s is a
solution of Eq. (<ref>) involving P_ν(a) and
P'_ν(a), we expect that P_ν^s(x) as a function of s has no
poles.  A rigorous proof of this claim is beyond the scope of this
paper.







 §.§ Solutions ν_n of Eq. (<ref>)



The above arguments are not directly applicable for the solutions
ν_n of Eq. (<ref>).  In fact, rewriting again the
right-hand side of Eq. (<ref>) in a form that is
consistent with Eq. (<ref>), one gets


    (ν-ν')(ν+ν'+1)/1-a^2∫_a^1 dx   P_ν(x)   P_ν'(x) 
           = P'_ν(a) ( P_ν'(a) - + - ν'(ν'+1)/q( - ν'(ν'+1)) P'_ν'(a) ) 
           - P'_ν'(a) ( P_ν(a) - + - ν(ν+1)/q( - ν(ν+1)) P'_ν(a) ) 
           -   (ν-ν')(ν+ν'+1)/q( - ν(ν+1))( - ν'(ν'+1)) P'_ν(a) P'_ν'(a) .

If ν and ν' are any two solutions of Eq. (<ref>),
the first two terms on the right-hand side vanish, but the last term
does not.  As a consequence, the functions P_ν(x) and
P_ν'(x), as well as all P_ν_n(x), determining c(x,t) and
J(t) via Eqs. (<ref>, <ref>), are not orthogonal to each
other.

At the same time, the above relation still allows us to prove that
all solutions of Eq. (<ref>) are real.  If there existed a
complex solution ν, then its complex conjugate, ν' = ν^*,
would also be a solution so that the above relation would become


    ∫_a^1 dx   |P_ν(x)|^2 = -   (1-a^2) |P'_ν(a)|^2/q | - ν(ν+1)|^2 .

However, this relation cannot hold because the left-hand side is the
integral of a positive function, whereas the right-hand side is
negative.








 §.§ Computation of the residues



The Laplace inversion formulas include the factor (∂_s
f_n)(s_n) which can be written as


    ∂_s f_n(s_n) = 1 + D/R^2 (2ν_n + 1) (∂_s ν_n)(s_n) ,

where ν_n is a shortcut notation for ν_n^s_n, i.e., the
n-th solution of Eq. (<ref>) evaluated at the pole s_n.
The last factor can be found by differentiating the boundary condition
(<ref>):


    ∂_s ν_n(s_n) = λ  P_ν_n(a)
     ( ∂_ν P'_ν(a) - q_s_n∂_ν P_ν(a))^-1_ν = ν_n/√(1-a^2) (s_n+λ+μ)^2 ,

so that one needs to evaluate the derivative with respect to the
degree ν.  Using Eqs. (<ref>,
<ref>) and the boundary condition (<ref>),
one gets


    (∂_ν P'_ν(a) - q_s_n∂_ν P_ν(a))_ν=ν_n  = P'_ν_n(a)/ν_n+1
        - 1/π[ν_n+1/a^2-1(A_ν_n+1(a) - a A_ν_n(a)) - q_s_n A_ν_n(a)] .

In practice, we used the integral representation (<ref>)
for evaluating Eq. (<ref>).





 §.§ Estimates for zeros ν_n^s



We discuss some properties of the solutions ν_n^s of
Eq. (<ref>) (with fixed a and q_s) that may facilitate
their numerical computation.  Even though most of the following
statements are expected to be classical, we could not find their
mathematical proofs in the literature.  As such proofs are beyond the
scope of the paper, most of the results of this section of the
Appendix remain conjectural from the mathematical point of view, even
so they have been checked numerically.  We emphasize that the
conjectural results of this section were not directly used in the main
text, except for guiding the author in the analysis.  As a
consequence, their conjectural character does not impair the quality
of main results.  In the following, we assume that q_s ≥ 0.

Since the operator Ł = ∂_x (1-x^2) ∂_x is
self-adjoint, the minimax principle should imply the standard
inequalities between Dirichlet, Robin, and Neumann eigenvalues (see
<cit.> for the case of the Laplace operator in
Euclidean domains), from which


    ν_n^N ≤ν_n^s ≤ν_n^D ,

where ν_n^D and ν_n^N are the zeros corresponding to the
limiting cases of Dirichlet (q_s = ∞) and Neumann (q_s = 0)
conditions:


    P_ν_n^D(a) = 0,       P'_ν_n^N(a) = 0 .

We also expect that ν_n^s are monotonous functions of a:


    n = ν_n^s(-1) ≤ν_n^s(a_1) ≤ν_n^s(a_2)       (a_1 < a_2) .



Using the asymptotic expansions of the Legendre function at large
ν <cit.>


    P_ν(cosθ)     = Γ(ν+1)/Γ(ν+3/2)(πsinθ/2)^-1/2
       ×cos((ν+1/2)θ - π/4) + O(ν^-1),

we get


    ν_n^D ≃ν_n^D, app = n + 3/4/1 - /π - 1/2     (n→∞),

where we used arccos(a) = π -.  Note that this asymptotic
relation is exact for a = 0 (or = π/2): ν_n^D = 2n+1.
We conjecture that for all n = 0,1,2,…,



    ν_n^D ≤ν_n^D, app     (a < 0) ,
          ν_n^D ≥ν_n^D, app     (a > 0) .


In other words, the asymptotic approximation ν_n^D, app is
either the upper bound (for a<0), or the lower bound (for a>0) for
the zeros ν_n^D.  These inequalities have been checked numerically
for several values of a and n.

Similarly, taking the derivative of P_ν(x) with respect to x,
one gets


    P'_ν(cosθ)    =   π/2Γ(ν+2)/Γ(ν+3/2)(πsinθ/2)^-3/2
       ×   sin((ν+1/2)θ - π/4) + O(ν^-1),

so that


    ν_n^N ≃ν_n^N, app = n + 1/4/1 - /π - 1/2     (n→∞).

This asymptotic relation becomes exact for a = 0 (or = π/2):
ν_n^N = 2n.  We conjecture that for all n = 0,1,2,…,



    ν_n^N ≥ν_n^N, app     (a < 0) ,
          ν_n^N ≤ν_n^N, app     (a > 0) .


In other words, the asymptotic approximation ν_n^N, app is
either the lower bound (for a<0), or the upper bound (for a>0) for
the zeros ν_n^N.  These inequalities have been checked
numerically for several values of a and n.

Combining these inequalities with Eq. (<ref>), we get the
following conjectural inequalities for a < 0:


    ν_n^N, app≤ν_n^N≤ν_n^s ≤ν_n^D ≤ν_n^D, app .

We note that the case a < 0 (small targets) is more relevant for
applications than a > 0 (large targets).  These inequalities allow
one to impose constraints on the range of ν, on which the n-th
zero ν_n^s is searched numerically.  At n = 0, the lower bound
is negative and should thus be replaced by 0.














§ EXAMPLE OF NUMERICAL IMPLEMENTATION



In this Appendix, we illustrate the computation of the total flux
J(t) shown in Fig. <ref>(b).  After truncation of the
spectral decomposition (<ref>) at = 20, this flux is
approximated by 20 exponential functions, determined by ν_n and
J_n, with n = 0,1,2,…,19.  For a given set of dimensionless
parameters (, κ̅, λ̅, and μ̅),
we first determined ν_n by solving Eq. (<ref>) by a
bisection method; then the coefficients J_n were computed from
Eq. (<ref>), as discussed in Sec. <ref>.

Figure <ref> shows the first 20 parameters ν_n and J_n
that were used for plotting Fig. <ref>(b) for a small target
= 0.01 and three configurations of reaction parameters.  For
convenience, ν_n are divided by n+1, to highlight their
asymptotic behavior ν_n ∝ n.  In all three cases, the
smallest solution ν_0 is small but strictly positive (its
dependence on  was shown in Fig. <ref>(a)).  As n
increases, ν_n/(n+1) approaches 1 from below, suggesting that
n+1 is the upper bound of ν_n.  Even though a rigorous proof of
this statement is beyond the scope of the paper (see also Appendix
<ref>), one can use the relation ν_n ∝ n to
estimate the validity range of the truncated decomposition.  In fact,
the omission of the terms with n ≥ is justified when Dt
ν_n(ν_n+1)/R^2 ≫ 1, i.e., Dt/R^2 ≫^-2.  This
inequality is satisfied for all figures in the paper.

For the considered three cases of reaction parameters, the
coefficients J_n are positive and decreasing with n.  The first
coefficient J_0 provides the dominant contribution to the total
flux, particularly for the case of a partial sink with irreversible
reaction (κ̅ = 1, λ̅ = μ̅ = 0).

















§ DERIVATION OF THE MFPT



The mean FPT to a partially reactive target on the surface of a sphere
was found in <cit.>.  Here we reproduce its computation for the
sake of completeness.  The MFPT τ(x) obeys the Poisson equation:
D Δτ = -1 that reads as


    ∂_x (1-x^2) ∂_x τ(x) = - R^2/D     (a < x < 1),

subject to the Robin boundary condition at x=a: D/Rsin
(∂_x τ)(a) = κτ(a).  A general solution is
obtained by integrating twice the above equation:


    τ(x) = R^2/2Dln(1-x^2) + c_1/2ln1+x/1-x + c_2 ,

with two constants c_1 and c_2.  To ensure the regularity of this
solution at x=1, we set c_1 = R^2/D so that


    τ(x) = R^2/Dln(1+x) + c_2 .

The constant c_2 is fixed by the Robin boundary condition, from
which Eq. (<ref>) follows.








§ ASYMPTOTIC BEHAVIOR FOR SMALL TARGETS



In this Appendix, we investigate separately the asymptotic behavior of
the solutions ν_n^s of Eq. (<ref>) and that of the
solutions ν_n of Eq. (<ref>) in the limit → 0.






 §.§ Solutions ν_n^s of Eq. (<ref>)


Here we analyze the asymptotic behavior of the solutions ν_n^s of
Eq. (<ref>) for a fixed q_s.  

For small , we have a ≃ -1 + ^2/2.  Using the Hobson's
representation (<ref>) and P_ν(-a) = P_ν(1-^2/2) = 1
+ O(^2), we get


    P_ν(a) ≃sin(πν)/π Ψ_(ν) + O(^2) ,

with


    Ψ_(ν) = ln(^2/4) + 2γ + 2ψ(ν+1) + π(πν),

where we used the reflection property of the digamma function:
ψ(1+z) - ψ(-z) = -π(πν).  Substituting this
expression into Eq. (<ref>) (which is equivalent to
Eq. (<ref>)), we get


    π(πν) + 2ψ(ν+1) ≃2/q_s ^2 + 2ln(2/) - 2γ ,

where we used that ψ(z+1) = ψ(z) + 1/z.  Since the right-hand
side of this relation is large as →0, ν should lie near the
poles of (πν).  We search thus an approximate solution of
Eq. (<ref>) as


    ν_n^s ≃ n + δ_n,     δ_n ≪ 1.

Note that ψ(ν+1) is regular at positive integers: ψ(n+1) =
H_n - γ, where H_n = 1/1 + … + 1/n are the harmonic
numbers (and H_0 = 0).  One gets thus


    1/δ_n≃2/q_s ^2 + 2ln(2/) - 2H_n .

This approximate solution is valid only when the right-hand side is
large (and δ_n is small), yielding a constraint of having not
too large n.

In what follows, we will need to evaluate P'_ν_n^s(a).  Using
again Eq. (<ref>), we find in the leading term


    P'_ν_n^s(a) ≃sin (πν_n)/π(1+a)≃2δ_n/1-a^2 .

Using Eqs. (<ref>, <ref>,
<ref>) and representations (<ref>,
<ref>), we also determine 


    (∂_s f_n)(s_n) ≃ 1       (→ 0).

Finally, given that b_n^2 = n+1/2 for = 0, this relation
should hold approximately for small .

In the following, we consider separately the cases of perfect (κ
= ∞) and imperfect (κ < ∞) target.








  §.§.§ Perfect sink


For a perfect sink (κ = ∞), one has q_s = ∞, and
the first term in Eq. (<ref>) vanishes, so that one
retrieves the asymptotic formula for the zeros of the Legendre
function from Ref. <cit.>:


    ν_n^s = ν_n ≃ n + 1/2ln(2/) - 2H_n .

For n > 0, the logarithmic term provides a small correction and can
thus be neglected so that ν_n ≃ n.  In contrast, when n =
0, the logarithmic term dominates.  One gets thus the asymptotic
behavior of c̃(x,s) and J̃(s) via
Eqs. (<ref>, <ref>).

Moreover, as q_s is infinite (and thus does not depend on s), the
degrees ν_n determine directly the poles via Eq. (<ref>).
In particular, one gets


    s_0 ≃ - D/2R^2ln(2/) .

Note that 1/|s_0| is close to the mean first passage time to the
target in this domain (see Appendix <ref>).
Substituting these results into Eq. (<ref>), we get the
long-time behavior (<ref>) of the total flux.




  §.§.§ Imperfect sink


For imperfect sink (κ < ∞), one can neglect the
logarithmic and constant terms in Eq. (<ref>), so that 


    ν_n^s ≃ n +  (s + μ)/2(s + λ + μ)  + o() .

This asymptotic relation determines c̃(x,s) and
J̃(s) via Eqs. (<ref>, <ref>).

Moreover, for irreversible reactions with λ = 0, q_s and
thus ν_n^s do not depend on s (and μ) so that one can use
the above relation to determine the asymptotic behavior of c(x,t)
and J(t) in Eqs. (<ref>, <ref>).  In fact,
Eq. (<ref>) yields


    s_n = - D/R^2(n(n+1) + (n+1/2)  + O(^2)) .

In a first approximation, one can neglect  in the higher-order
eigenmodes.  Keeping only the term with n = 0 with the smallest rate
|s_0|, we get the long-time behavior (<ref>) of the total
flux.





 §.§ Solutions ν_n of Eq. (<ref>)


For reversible reactions (with λ > 0), c(x,t) and J(t) are
determined by the solutions ν_n of Eq. (<ref>).  To
get their asymptotic behavior in the small-target limit, we rewrite
Eq. (<ref>) as


    +  - ν(ν+1) =  ( - ν(ν+1))  P_ν(a)/√(1-a^2)  P'_ν(a) .

The Hobson representation (<ref>) for P_ν(a) and its
derivative P'_ν(a) yields in the limit → 0:


    +  - ν(ν+1) =    ( - ν(ν+1)) [Ψ_(ν) + O(^2)]/2 + ^2 ν(ν+1)/2Ψ_(ν) + O(^2) ,

with Ψ_(ν) given by Eq. (<ref>).  As the function
Ψ(ν) diverges at any integer ν, one can search again for
solutions near integers, ν_n = n + δ_n, from which 


    Ψ_(ν_n) = 1/δ_n + ln(^2/4) + 2H_n + O(δ)

and thus


    δ_n ≃( - n(n+1))/2(+ -n(n+1))  + o().

For instance, one gets the following asymptotic relation for the
smallest degree:


    ν_0 = δ_0 ≈  /2( + )≈μκρ/2D(μ + λ) .

In the particular case when + = n'(n'+1) for some
integer n', Eq. (<ref>) does not hold for this n',
and one gets a slower approach to the limit:


    δ_n'≃( /2(2n'+1))^1/2 √() + o(√()).


Yet another solution of Eq. (<ref>) is possible when its
left-hand side is of the order of  and the corresponding ν is
not close to an integer.  In this case, Ψ_(ν) behaves
asymptotically as 2ln(/2) + O(1) and thus one gets the solution
in the leading order in :


    ν(ν+1) =  +  -  ln (1/) + O() .

As a consequence, the degree ν approaches the value √(1/4 +
+)-1/2 in the limit → 0.  We recall that this
value should not be integer (otherwise the relation
(<ref>) does not necessarily hold).  This solution
yields the eigenvalue, which is close to +.  We
emphasize that the above argument is not applicable when λ =
0, in which case q_s does not depend on s any more.


















89

Albert		B. Alberts, A. Johnson, J. Lewis, D. Morgan, M. Raff, K. Roberts, and P. Walter, 
				Molecular Biology of the Cell 
				(Garland Science, New York, NY, 2014).


Sungkaworn17		T. Sungkaworn, M.-L. Jobin, K. Burnecki, A. Weron, M. J. Lohse, and D. Calebiro,
				“Single-molecule imaging reveals receptor-G protein interactions at cell surface hot spots”,
				Nature 550, 543-547 (2017).

Bloomfield79		V. A. Bloomfield and S. Prager,
				“Diffusion-controlled reactions on spherical surfaces. Application to bacteriophage tail fiber attachment”,
				Biophys. J. 27, 447-453 (1979).

Chao81		N. Chao, S. H. Young, and M. Poo, 
				“Localization of Cell Membrane Components by Surface Diffusion into a Trap”, 
				Biophys. J. 36, 139-153 (1981).

Weaver83		D. L. Weaver,
				“Diffusion-mediated localization on membrane surfaces”,
				Biophys. J. 41, 81-86 (1983).

Linderman86		J. J. Linderman and D. A. Lauffenburger,
				“Analysis of intracellular receptor/ligand sorting. Calculation of mean surface and bulk diffusion times within a sphere”,
				Biophys. J. 50, 295-305 (1986).

Sano81		H. Sano and  M. Tachiya,
				“Theory of diffusion-controlled reactions on spherical surfaces and its application to reactions on micellar surfaces”,
				J. Chem. Phys. 75, 2870 (1981).


Berg85		O. G. Berg,
				“Orientation constraints in diffusion-limited macromolecular association. 
				The role of surface diffusion as a rate-enhancing mechanism”,
				Biophys. J. 47, 1 (1985).


Benichou10		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).

Benichou11		O. Bénichou, D. S. Grebenkov, P. Levitz, C. Loverdo, and R. Voituriez,
				“Mean First-Passage Time of Surface-Mediated Diffusion in Spherical Domains”,
				J. Stat. Phys. 142, 657-685 (2011).

Rojo11		F. Rojo and C. E. Budde,
				“Enhanced diffusion through surface excursion: A master-equation approach to the narrow-escape-time problem”,
				Phys. Rev. E 84, 021117 (2011).

Rupprecht12a		J.-F. Rupprecht, O. Bénichou, D. S. Grebenkov, and R. Voituriez,
				“Kinetics of Active Surface-Mediated Diffusion in Spherically Symmetric Domains”,
				J. Stat. Phys. 147, 891-918 (2012).

Rupprecht12b		J.-F. Rupprecht, O. Bénichou, D. S. Grebenkov, and R. Voituriez,
				“Exact mean exit time for surface-mediated diffusion”,
				Phys. Rev. E 86, 041135 (2012).







Berg77		H. C. Berg and E. M. Purcell,
				“Physics of chemoreception”,
				Biophys. J. 20, 193-239 (1977).


Zwanzig90		R. Zwanzig, 
				“Diffusion-controlled ligand binding to spheres partially covered by receptors: an effective medium treatment”,
				Proc. Natl. Acad. Sci. USA 87, 5856 (1990).

Zwanzig91		R. Zwanzig and A. Szabo,
				“Time dependent rate of diffusion-influenced ligand binding to receptors on cell surfaces”,
				Biophys. J. 60, 671 (1991).

Berezhkovskii04	A. Berezhkovskii, Y. Makhnovskii, M. Monine, V. Zitserman, and S. Shvartsman, 
				“Boundary homogenization for trapping by patchy surfaces”,
				J. Chem. Phys. 121, 11390 (2004).

Berezhkovskii06	A. M. Berezhkovskii, M. I. Monine, C. B. Muratov, and S. Y. Shvartsman, 
				“Homogenization of boundary conditions for surfaces with regular arrays of traps”,
				J. Chem. Phys. 124, 036103 (2006).

Muratov08		C. Muratov and S. Shvartsman, 
				“Boundary homogenization for periodic arrays of absorbers”,
				Multiscale Model. Simul. 7, 44-61 (2008).

Dagdug16		L. Dagdug, M. Vázquez, A. Berezhkovskii, and V. Zitserman, 
				“Boundary homogenization for a sphere with an absorbing cap of arbitrary size”,
				J. Chem. Phys. 145, 214101 (2016).

Lindsay17		A. E. Lindsay, A. J. Bernoff, and M. J. Ward, 
				“First Passage Statistics for the Capture of a Brownian Particle by a 
				Structured Spherical Target with Multiple Surface Traps”,
				Multiscale Model. Simul. 15, 74-109 (2017).

Bernoff18a		A. J. Bernoff and A. E. Lindsay, 
				“Numerical approximation of diffusive capture rates by planar and spherical surfaces with absorbing pores”,
				SIAM J. Appl. Math. 78, 266-290 (2018).

Bernoff18b		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).



Shoup81		D. Shoup, G. Lipari, and A. Szabo,
				“Diffusion-controlled bimolecular reaction rates. The effect of rotational diffusion and orientation constraints”,
				Biophys. J. 36, 697-714 (1981).

Grebenkov17		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).




Singer06c		A. Singer, Z. Schuss, and D. Holcman,
				“Narrow Escape, Part III Riemann surfaces and non-smooth domains”,
				J. Stat. Phys. 122, 491-509 (2006).



Coombs09		D. Coombs, R. Straube, and M. J. Ward,
				“Diffusion on a Sphere with Localized Traps: Mean First Passage Time, Eigenvalue Asymptotics, and Fekete Points”,
				SIAM J. Appl. Math. 70, 302 (2009).

Prustel13		T. Prüstel and M. Tachiya,
				“Reversible diffusion-influenced reactions of an isolated pair on some two dimensional surfaces”,
				J. Chem. Phys. 139, 194103 (2013).








Filoche99		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).

Grebenkov03		D. S. Grebenkov, M. Filoche, and B. Sapoval, 
				“Spectral Properties of the Brownian Self-Transport Operator”, 
				Eur. Phys. J. B 36, 221-231 (2003).

Grebenkov06		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). 

Grebenkov07a		D. S. Grebenkov,
				“Residence times and other functionals of reflected Brownian motion”,
				Phys. Rev. E 76, 041139 (2007). 

Singer08		A. Singer, Z. Schuss, Osipov, and D. Holcman,
				“Partially Reflected Diffusion”,
				SIAM J. Appl. Math. 68, 844 (2008).


Shoup82		D. Shoup and A. Szabo, 
				“Role of diffusion in ligand binding to macromolecules and cell-bound receptors”, 
				Biophys. J. 40, 33 (1982).

Lauffenburger		D. A. Lauffenburger and J. Linderman, 
				Receptors: Models for Binding, Trafficking, and Signaling
				(Oxford University Press, 1993).




Agmon90		N. Agmon and A. Szabo,
				“Theory of reversible diffusion-influenced reactions”,
				J. Chem. Phys. 92, 5270 (1990).




Tachiya80		M. Tachiya, 
				“Theory of diffusion-controlled dissociation and its applications to charge separation,” 
				in Extended Abstract of Annual Meeting on Photochemistry, Tsu, Japan, 1980 
				(Japan. Photochem. Assoc., 1980), pp. 256-257.


Agmon84		N. Agmon, 
				“Diffusion with back reaction”,
				J. Chem. Phys. 81, 2811 (1984).

Kim99			H. Kim and K. J. Shin,
				“Exact Solution of the Reversible Diffusion-Influenced Reaction for an Isolated Pair in Three Dimensions”,
				Phys. Rev. Lett. 82, 1578 (1999).







Goodrich54		F. C. Goodrich, 
				“Random walk with semi-adsorbing barrier”,
				J. Chem. Phys. 22, 588-594 (1954).



Collins49		F. C. Collins and G. E. Kimball,
				“Diffusion-controlled reaction rates”,
				J. Coll. Sci. 4, 425 (1949).

Sano79		H. Sano and M. Tachiya,
				“Partially diffusion-controlled recombination”,
				J. Chem. Phys. 71, 1276 (1979).

Sapoval94		B. Sapoval,
				“General Formulation of Laplacian Transfer Across Irregular Surfaces”,
				Phys. Rev. Lett. 73, 3314 (1994).


Grebenkov19c		D. S. Grebenkov,
				“Imperfect Diffusion-Controlled Reactions”, 
				in Chemical Kinetics: Beyond the Textbook, Eds. K. Lindenberg, R. Metzler, and G. Oshanin 
				(World Scientific, 2019).



Grebenkov19d		D. S. Grebenkov,
				“Spectral theory of imperfect diffusion-controlled
				reactions on heterogeneous catalytic surfaces”,
				J. Chem. Phys. 151, 104108 (2019).



Redner		S. Redner, 
				A Guide to First Passage Processes 
				(Cambridge: Cambridge University press, 2001).


Prustel12		T. Prüstel and M. Meier-Schellersheim,
				“Exact Green's function of the reversible diffusion-influenced reaction for an isolated pair in two dimensions”,
				J. Chem. Phys. 137, 054104 (2012).

Prustel13b		T. Prüstel and M. Meier-Schellersheim,
				“Theory of reversible diffusion-influenced reactions with non-Markoviandissociation in two space dimensions”,
				J. Chem. Phys. 138, 104112 (2013).





Metzler		R. Metzler, G. Oshanin, and S. Redner (Eds.) 
				First-Passage Phenomena and Their Applications
				(Singapore: World Scientific, 2014).

Holcman14		D. Holcman and Z. Schuss,
				“The Narrow Escape Problem”,
				SIAM Rev. 56, 213-257 (2014).






Torney83		D. C. Torney and H. M. McConnel,
				“Diffusion-Limited Reaction Rate Theory for Two-Dimensional Systems”,
				Proc. Royal Soc. London A 387, 147-170 (1983).



Carslaw		H. S. Carslaw and J. C. Jaeger,
				Conduction of Heat in Solids, 2nd Ed.
				(Oxford University Press, 1959).

Crank			J. Crank, 
				The Mathematics of Diffusion, 2nd Ed.
				(Oxford University Press, 1975).

Thambynayagam		R. K. M. Thambynayagam,
				The Diffusion Handbook: Applied Solutions for Engineers
				(New York: McGraw-Hill Education, 2011).




Smoluchowski17	M. Smoluchowski, 
				“Versuch einer Mathematischen Theorie der Koagulations Kinetic Kolloider Lösungen”,
				Z. Phys. Chem. 129, 129-168 (1917).

Galanti16		M. Galanti, D. Fanelli, S. D. Traytak, and F. Piazza, 
				“Theory of diffusion-influenced reactions in complex geometries”,
				Phys. Chem. Chem. Phys. 18, 15950 (2016).

Traytak18		S. D. Traytak and D. S. Grebenkov, 
				“Diffusion-influenced reaction rates for active “sphere-prolate spheroid” pairs and Janus dimers”, 
				J. Chem. Phys. 148, 024107 (2018).

Grebenkov18e		D. S. Grebenkov and D. Krapf, 
				“Steady-state reaction rate of diffusion-controlled reactions in sheets”, 
				J. Chem. Phys 149, 064117 (2018).

Grebenkov19a		D. S. Grebenkov and S. Traytak, 
				“Semi-analytical computation of Laplacian Green functions in three-dimensional 
				domains with disconnected spherical boundaries”, 
				J. Comput. Phys. 379, 91-117 (2019).




Benichou10a		O. Bénichou, C. Chevalier, J. Klafter, B. Meyer, and R. Voituriez,
				“Geometry-controlled kinetics”,
				Nature Chem. 2, 472-477 (2010).

Benichou14		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).




Godec16		A. Godec and R. Metzler,
				“Universal Proximity Effect in Target Search Kinetics in the Few-Encounter Limit”,
				Phys. Rev. X 6, 041037 (2016).

Grebenkov18a		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).

Grebenkov18b		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). 


Bressloff08		P. C. Bressloff, B. A. Earnshaw, and M. J. Ward,
				“Diffusion of protein receptors on a cylindrical dendritic membrane with partially absorbing traps”,
				SIAM J. Appl. Math. 68, 1223 (2008).

Novak07		I. L. Novak, F. Gao, Y.-S. Choi, D. Resasco, J. C. Schaff, and B. M. Slepchenko,
				“Diffusion on a Curved Surface Coupled to Diffusion in the Volume: Application to Cell Biology”,
				J. Comput. Phys. 226, 1271-1290 (2007).



Simons97		K. Simons and E. Ikonen,
				“Functional rafts in cell membranes”,
				Nature 387, 569-572 (1997).

Masson09		J.-B. Masson, D. Casanova, S. Türkcan, G. Voisinne, M. R. Popoff, M. Vergassola, and A. Alexandrou,
				“Inferring maps of forces inside cell membrane microdomains”,
				Phys. Rev. Lett. 102, 048103 (2009).

Lingwood10		D. Lingwood and K. Simons,
				“Lipid rafts as a membrane organizing principle”,
				Science 327, 46-50 (2010).



Grebenkov17c		D. S. Grebenkov, 
				“First passage times for multiple particles with reversible target-binding kinetics”, 
				J. Chem. Phys. 147, 134112 (2017). 

Lawley19		S. D. Lawley and J. B. Madrid,
				“First passage time distribution of multiple impatient particles with reversible binding”,
				J. Chem. Phys. 150, 214113 (2019).



Ritchie05		K. Ritchie, X.-Y. Shan, J. Kondo, K. Iwasawa, T. Fujiwara, and A. Kusumi,
				“Detection of Non-Brownian Diffusion in the Cell Membrane in Single Molecule Tracking”,
				Biophys. J. 88, 2266-2277 (2005).

Yamamoto14		E. Yamamoto, T. Akimoto, M. Yasui, and K. Yasuoka,
				Origin of subdiffusion of water molecules on cell membrane surfaces,
				Scient. Rep. 4, 4720 (2014).

He16			W. He, H. Song, Y. Su, L. Geng, B. J. Ackerson, H. B. Peng, and P. Tong,
				Dynamic heterogeneity and non-Gaussian statistics for acetylcholine receptors on live cell membrane,
				Nat. Comm. 7, 11701 (2016).
















Chechkin17		A. V. Chechkin, F. Seno, R. Metzler, and I. M. Sokolov,
				“Brownian yet Non-Gaussian Diffusion: From Superstatistics to Subordination of Diffusing Diffusivities”,
				Phys. Rev. X 7, 021002 (2017).

Lanoiselee18		Y. Lanoiselée, N. Moutal, and D. S. Grebenkov, 
				“Diffusion-limited reactions in dynamic heterogeneous media”, 
				Nat. Comm. 9, 4398 (2018).

Grebenkov19		D. S. Grebenkov, 
				“A unifying approach to first-passage time distributions in diffusing diffusivity and switching diffusion models,” 
				J. Phys. A 52, 174001 (2019).

Lanoiselee19		Y. Lanoiselée and D. S. Grebenkov, 
				“Non-Gaussian diffusion of mixed origins”, 
				J. Phys. A 52, 304001 (2019).







Erdelyi		A. Erdélyi (ed.) 
				Higher transcendental functions, Vol. 1 
				(McGraw-Hill Book Company , New York, 1953).



Hobson		E. W. Hobson, 
				The Theory of Spherical and Ellipsoidal Harmonics
				(Cambridge University Press, London, 1931).

Hall49		R. N. Hall, 
				“The application of non-integral Legendre functions to potential problems”,
				J. Appl. Phys. 20, 925-936 (1949).

Smythe		W. R. Smythe, 
				Static and Dynamic Electricity 
				(McGraw-Hill Book Company, Inc. New York, 1939)
				

Maier16		R. S. Maier,
				“Legendre functions of fractional degree: transformations and evaluations”,
				Proc. R. Soc. A 472, 20160097 (2016).


Schelkunoff41		S. A. Schelkunoff,
				“Theory of antennas of arbitrary size and shape”,
				Proc. I.R.E. 29, 493 (1941).




Abramowitz		M. Abramowitz and I. A. Stegun,
				Handbook of Mathematical Functions
				(Dover Publisher, New York, 1965).

MacDonald1899		H. M. MacDonald,
				“Zeroes of the Spherical Hormonic Pnm(mu) considerd as a Function of n”,
				Proc. Lond. Math. Soc. 31, 264-281 (1899).


DLMF			Digital Library of Mathematical Functions (NIST project),
				Legendre and Related Functions (<https://dlmf.nist.gov/14.11>)

Szmytkowski06		R. Szmytkowski,
				“On the derivative of the Legendre function of the first kind with respect to its degree”,
				J. Phys. A 39, 15147-15172 (2006).

Cohl11		H. S. Cohl,
				“On Parameter Differentiation for Integral Representations of Associated Legendre Functions”,
				SIGMA 7, 050 (2011).




Courant		R. Courant and D. Hilbert,
				Methods of Mathematical Physics 
				(Wiley-VCH, New York, 1989).

Grebenkov13		D. S. Grebenkov and B.-T. Nguyen, 
				“Geometrical structure of Laplacian eigenfunctions”, 
				SIAM Reviews 55, 601-667 (2013).








