Diffusion-influenced reactions on non-spherical partially absorbing axisymmetric surfaces

TL;DR

This paper develops a perturbative analytical approach to estimate diffusion-controlled reaction rates on non-spherical, partially absorbing surfaces, validated through exact, approximate, and numerical solutions for spheroids.

Contribution

It introduces a simple first-order correction formula for reaction rates on non-spherical surfaces, extending the spherical case and validating it with numerical methods.

Findings

The correction formula accurately predicts reaction rates for non-spherical shapes.

Numerical solutions confirm the validity range of the perturbative approach.

The method quantifies how non-sphericity influences reaction rates.

Abstract

The calculation of the diffusion-controlled reaction rate for partially absorbing, non-spherical boundaries presents a formidable problem of broad relevance. In this paper we take the reference case of a spherical boundary and work out a perturbative approach to get a simple analytical formula for the first-order correction to the diffusive flux onto a non-spherical partially absorbing surface of revolution. To assess the range of validity of this formula, we derive exact and approximate expressions for the reaction rate in the case of partially absorbing prolate and oblate spheroids. We also present numerical solutions by a finite-element method that extend the validity analysis beyond spheroidal shapes. Our perturbative solution provides a handy way to quantify the effect of non-sphericity on the rate of capture in the general case of partial surface reactivity.

Equations

$$\text{ce}B+S->[k]S+P$$

$$k_S=4\pi DR{c}_{\infty }$$

$$x$$

$$y$$

$$z$$

$$r(\theta )=R+\delta f(\theta ),\delta \ll R$$

$${\nabla }^{2}c=0\vec{r}\in {\mathbb{R}}^3\setminus T$$

$$D\nabla c\cdot \hat{n}{|}_{\partial \Omega }=\kappa c{|}_{\partial \Omega }$$

$$\underset{||\vec{r}||\to \infty }{\lim }c=c_{\infty }$$

$$\hat{n}=\frac{{\vec{t}}_n}{||{\vec{t}}_n||}=\frac{{\vec{t}}_{\theta }\wedge {\vec{t}}_{\varphi }}{||{\vec{t}}_{\theta }\wedge {\vec{t}}_{\varphi }||}$$

$$k=-{\int }_{\partial \Omega }\vec{J}\cdot \hat{n}dS=2\pi D{\int }_0^{\pi }\nabla c\cdot \hat{n}d\theta$$

$${\nabla }^{2}u=0\vec{\rho }\in {\mathbb{R}}^3\setminus T^{\prime }$$

$$\displaystyle\left.\biggl{(}\nabla u\cdot\hat{n}+h(1-u)\biggr{)}\right|_{\partial\Omega^{\prime}}=0$$

$$\underset{\rho \to \infty }{\lim }u=0$$

$$u(\rho ,\theta )=u_0(\rho )+ϵu_1(\rho ,\theta )+\mathcal{O}({ϵ}^2)$$

$$u_0(\rho )=\frac{A_0}{\rho }$$

$$u_1(\rho ,\theta )=\sum _{\ell =0}^{\infty }{B}_{\ell }{\rho }^{-(\ell +1)}{P}_{\ell }(\cos \theta )$$

$${\vec{t}}_n$$

$$u{|}_{\partial {\Omega }^{\prime }}$$

$$\nabla u{|}_{\partial {\Omega }^{\prime }}$$

$$A_0=\frac{h}{1+h}$$

$$\sum _{\ell =0}^{\infty }(\ell +h+1)B_{\ell }{P}_{\ell }(\cos \theta )=\frac{h(2+h)}{1+h}f(\theta )$$

$${\int }_{-1}^1{P}_{\ell }(\mu )P_{{\ell }^{\prime }}(\mu )d\mu =\frac{2}{2\ell +1}{\delta }_{\ell {\ell }^{\prime }}$$

$$B_{\ell }=\frac{(2\ell +1)(2+h)h}{2(\ell +1+h)(1+h)}{\int }_0^{\pi }{P}_{\ell }(\cos \theta )f(\theta )\sin \theta d\theta$$

$$\frac{k}{k_S}$$

$$B_0^{\infty }=\frac{1}{2}{\int }_0^{\pi }f(\theta )\sin \theta d\theta$$

$$\frac{k}{k_S}=1+ϵB_0^{\infty }+\mathcal{O}({ϵ}^2)$$

$$\frac{1}{k}$$

$$\text{Da}=\frac{\kappa R}{D}=\frac{k_r}{k_S}=h$$

$$S$$

$$f(\theta )=\{\begin{array}{ll}-{\sin }^2\theta & \text{prolate}\\ -{\cos }^2\theta & \text{oblate}\end{array}$$

$$f(\theta )=\{\begin{array}{ll}{\cos }^2\theta & \text{prolate}\\ {\sin }^2\theta & \text{oblate}\end{array}$$

$$\frac{k}{k_S}=\{\begin{array}{ll}1-\frac{2}{3}ϵ+\mathcal{O}({ϵ}^2)& \text{prolate}\\ 1-\frac{1}{3}ϵ+\mathcal{O}({ϵ}^2)& \text{oblate}\end{array}$$

$$\frac{k}{k_S}=\{\begin{array}{ll}1+\frac{1}{3}ϵ+\mathcal{O}({ϵ}^2)& \text{prolate}\\ 1+\frac{2}{3}ϵ+\mathcal{O}({ϵ}^2)& \text{oblate}\end{array}$$

$$\frac{k}{k_{S}}=\begin{cases}\displaystyle\frac{a_{E}}{b}\biggl{(}\frac{1}{2}\ln\biggl{(}\frac{1+a_{E}/b}{1-a_{E}/b}\biggr{)}+\frac{D}{\kappa b}\frac{\arcsin(a_{E}/b)}{(1-a_{E}^{2}/b^{2})^{1/2}}\biggr{)}^{-1}\hskip 11.38109pt\text{prolate}\\ \displaystyle\frac{a_{E}/b}{\arcsin(a_{E}/b)}\,\frac{1}{1+D/(\kappa b)}\hskip 93.89409pt\text{oblate}\\ \end{cases}$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: a* Centre de Biophysique Moléculaire (CBM) CNRS UPR4301 $\&$ Université d’Orléans, Orléans 45071, France; Tel: +33 238 255653; E-mail: Francesco.Piazza@cnrs-orleans.fr*††footnotetext: b* Laboratory of Condensed Matter Physics, CNRS – Ecole Polytechnique, IP Paris, F-91128 Palaiseau, France. *

1 Introduction

Reaction-diffusion processes govern the behavior of complex systems in many contexts ranging from chemistry and nanosciences 1 to bio-engineering 2, 3, 4, 5. They are also key in biology by controlling enzyme catalysis 6, 3, antigen-antibody encounter 7, 8, ligand-receptor binding 9, 10, 11, 12, 13, 14, fluorescence quenching 15, cellular nutrient uptake 16, 17, 18, 19, and oxygen uptake in lungs and placenta 20, 21, 22.

In many cases, diffusion relaxation times are fast with respect to typical time scales of interest. For example, the local concentration field of ligand molecules around a membrane receptor equilibrates more swiftly than lateral diffusion of the latter within the membrane and more rapidly than (or as fast as) large-scale conformational changes of the protein. Along the same lines, the concentration of a small substrate molecule equilibrates fast around a large core-shell nanoreactor 5 whose hydrogel shell carries a number of immobilized metal nanocatalysts. As a consequence, the problem can often be reduced to solving the stationary diffusion (i.e. Laplace) equation 2 in a non-equilibrium setting (i.e. corresponding to a steady non-zero flux of molecules from the bulk to the surface). In the simplest case, one is led to computing the total molecular flux to a reactive surface, modeling for example the uptake of nutrient molecules by a colony of algae or the binding of a ligand onto the receptor-covered cell surface.

In more mathematical terms, we consider the pseudo first-order contact reaction

[TABLE]

between point-like molecules $B$ that diffuse in solution with diffusion coefficient $D$ and static reactive surfaces $S$ that catalyze the transformation of $B$ into some inert product $P$. When the characteristic relaxation time for diffusive processes is small enough, the kinetics is essentially controlled by the steady-state rate $k$ 2. If $B$ and $S$ species are dilute enough and the concentration of $S$ is much smaller than the bulk concentration of $B$, $c_{\infty}$, the rate $k$ can be obtained by solving the Laplace equation with appropriate boundary conditions (BC) enforced on the surface $S$ and computing the total steady-state diffusive flux of $B$ molecules into an isolated $S$ 2, 23, 24, 25. When $S$ is a sphere of radius $R$, this leads to the well-known Smoluchowski rate 26, 27,

[TABLE]

which expresses the very blueprint of diffusion to capture: the number of particles absorbed at the sink surface scales linearly with the size of its surface.

In real-life applications, reactive boundaries tend to be more complex than perfectly absorbing spheres. On one hand, many interesting applications feature partially reactive surfaces (see a recent overview in 28). The partial reactivity can describe, e.g., the need for a molecule to overcome an energy activation barrier to react 29, or heterogeneous distribution of reactive patches on the otherwise inert boundary when small swiftly diffusing molecules can only react at specific parts of the boundary 30, 31, 32, 33, 34, 35, 36, 37. This is the case for example of a more realistic treatment of a cell surface 9, 10, 11, 12, 13, where receptor molecules only constitute small absorbing patches (or clusters) on an otherwise inert surface***We refer in this setting to the inability of the diffusing ligand to establish bonds at generic spots on the surface, as opposed to specific high-affinity bonds at locations occupied receptors. This does not exclude that the surface might be interacting non-specifically (and less strongly) with the ligand at generic locations.. Partial reactivity can also describe stochastic gating of exchange channels (like aquaporins) on the plasma membrane 38, 39, 40, 41, partial recombination 42, reactions on micelles 43, and transfer across semi-permeable membranes 44, 45. On the other hand, in many applications one has to deal with non-spherical boundaries 8. This in general makes the problem of determining the diffusive flux very hard already in the simpler case of a perfect sink 46, 47, the case of partially absorbing non-spherical surfaces being in general a rather formidable task. This is the case for example of binding to a protein, featuring a few specific binding pockets (local absorbing patches) on an overall reflecting surface that is most of the times non-spherical.

In this paper, we tackle the hard problem of computing the flux to non-spherical surfaces perturbatively. Our calculations lead to a handy, closed-form expression for the first-order correction to the Smoluchowski rate for arbitrary, partially absorbing axially symmetric perturbations of a spherical boundary. The paper is organized as follows. In section 2 we state the problem and work out the perturbative solution. In section 3 we report a detailed analysis of the validity regime of the first-order correction. For this purpose, the diffusive flux to partially reactive prolate and oblate spheroids is computed. We also present numerical computations for more sophisticated surfaces. Finally, we wrap our results and summarize in section 4.

2 Methodology

Let us consider the surface of revolution $\partial\Omega$, parameterized in spherical coordinates $(r,\theta,\phi)$ as (Fig. 1)

[TABLE]

with $\theta\in[0,\pi]$, $\phi\in[0,2\pi)$ and

[TABLE]

The idea is to consider an arbitrary small modulation of the sphere $S_{R}$ of radius $R$, as specified by a given function $f(\theta)$. We aim to compute the steady-state diffusive flux of point-like $B$ particles into the partially absorbing surface $\partial\Omega$. As put forward by Collins and Kimball 48, this situation can be accommodated for within the framework of a boundary value problem by endowing the reactive surface with a constant intrinsic reactivity $\kappa$ (in units m/s). Whenever a diffusing $B$ particle hits the surface, the reaction (1) occurs with a probability proportional to $\kappa$ 49, 50, 51. From a mathematical standpoint, one needs to compute the steady-state concentration field of $B$ particles, $c(\vec{r})$, satisfying the following boundary value problem:

[TABLE]

where $T$ is the interior of the surface $\partial\Omega$. The Robin boundary condition (5b) means that the diffusive flux density should equal the reactive flux density at each point of the surface. In this paper, we take the convention of considering the normal pointing away from the boundary into the bulk; in other words, we consider the inward normal vector, which can be computed as the cross product of the two local tangent vectors (see Fig. 1):

[TABLE]

In the limit $\kappa\to\infty$ the boundary becomes a sink (Dirichlet BC), i.e. $c|_{\partial\Omega}=0$, while in the opposite limit $\kappa\to 0$ (von Neumann BC), the surface becomes perfectly reflecting and no reaction occurs. The regularity condition (5c) means that the concentration field of $B$ molecules far from the boundary is constant and equal to the bulk concentration. This ensures that after a relaxation transient period, the system settles in a nonequilibrium steady state with a net diffusive flux onto the reactive boundary. Accordingly, the reaction rate can be computed as the total diffusive flux into the surface, i.e.

[TABLE]

where $\vec{J}=-D\nabla c$ is the flux density. It is expedient to introduce non-dimensional variables, $\vec{\rho}=\vec{r}/R$, $u(\vec{\rho})=1-c(\vec{\rho}R)/c_{\infty}$, $\epsilon=\delta/R\ll 1$, and $h=\kappa R/D$. The problem to be solved reads then

[TABLE]

where $\rho=||\vec{\rho}||$, and prime emphasizes that the surface $\partial\Omega$ and its interior $T$ are also rescaled.

The above problem can be solved perturbatively by looking for a solution in the form of a perturbative expansion in powers of $\epsilon$. Furthermore, the symmetry of the problem requires the solution to be axially symmetric. To the first order, we set then

[TABLE]

where $u_{0}$ and $u_{1}$ are both harmonic functions that can be expressed in general as

[TABLE]

Here $P_{\ell}(\cos\theta)$ denotes the Legendre polynomial of order $\ell$ and $\{A_{0},B_{0},B_{1},\dots\}$ is an infinite set of constants that are fixed by the boundary condition (8b).

In order to impose the BC (8b), we need to work out the gradient $\nabla u$ and the concentration field $u$ on the surface $\partial\Omega^{\prime}$, as well as the components of the inward normal vector $\vec{t}_{n}$. From the definition (6), it is straightforward to show that

[TABLE]

where $\hat{e}_{r}$ and $\hat{e}_{\theta}$ are the radial and polar tangent unit vectors relative to the sphere $S_{R}$ (see Fig. 1). With the help of Eqs. (11a) and (11c), the boundary conditions (8b) and (8c) give

[TABLE]

Multiplying Eq. (12b) by $P_{m}(\cos\theta)\sin\theta$, integrating and recalling the orthonormality of Legendre polynomials,

[TABLE]

one is immediately led to

[TABLE]

3 Results and discussion

The rate should be computed from Eq. (7). Taking into account Eqs. (11a), (11c) and (13), we get

[TABLE]

where

[TABLE]

Formula (15) is the main result of this paper, describing the first-order correction to the Smoluchowski rate for a non-spherical, axisymmetric partially reactive boundary. In the limit of a perfect sink $h\to\infty$, this becomes

[TABLE]

It is instructive to observe that the flux to the partially absorbing boundary and the flux to the sink are not related to the first order in $\epsilon$ through the usual additivity prescription for independent probabilities (per unit time):

[TABLE]

where $k_{r}=hk_{S}=4\pi R^{2}\kappa$ has dimensions of a rate constant and gauges the uniform chemical activity of the boundary. Eqs. (18) and (15) coincide to the first order in $\epsilon$ only for highly reactive surfaces ($h\gg 1$). It is worthwhile to mention that the information carried by the intrinsic reaction rate $k_{r}$ it is also known as Damköhler number Da,

[TABLE]

which quantifies to what extent the overall reaction rate is limited by diffusion.

3.1 Validity range of the perturbative calculation

In order to determine the validity range of our results, we consider the case of spheroidal perturbations (i.e. ellipsoids of rotation). As an example, Fig. 2 illustrates a prolate spheroidal perturbation of the form $(a,a,b)$ with $a<b$ that is obtained from a sphere of radius $R$ either by extension along the $z$-axis into $(R,R,R+\delta)$, or by compression at the equator into $(R-\delta,R-\delta,R)$. In the compression case, the surface area of the resulting spheroid is reduced as compared to that of the sphere $S_{R}$, whereas it increases in the extension case. Accordingly, the first-order correction to the Smoluchowski rate is expected to be negative in the compression case and positive in the extension case. A similar reasoning applies to the case of oblate spheroidal perturbations of the form $(b,b,a)$ with $a<b$.

3.2 Perfectly absorbing spheroids

In order to illustrate our calculations, we first treat the case of spheroidal perturbations of a spherical sink. Let us consider a perfectly absorbing sphere (i.e. $h\to\infty$) with radius $R=b$. With reference to Fig. 2 (right), the surface area of a prolate spheroid with axes $(a,a,b)$, with $a=R-\delta$, reads

[TABLE]

with $\epsilon=\delta/R=1-a/b$. An equivalent analysis can be performed in the case of oblate spheroids. Overall, it is not difficult to see (cf. also Appendix A) that the perturbation modulating function $f(\theta)$ for spheroidal perturbations of the compression type reads

[TABLE]

while for perturbation of the extension type, one has

[TABLE]

Taking into account expressions (21), the prescription (17) gives immediately the negative corrections for the compression case,

[TABLE]

where $k_{S}=4\pi Dbc_{\infty}$. The expressions (23) are in agreement with the known exact results (see 52, 46, 47 and Appendices). The corrections for the extension case can be computed easily from Eqs. (22),

[TABLE]

where the unperturbed rate is now $k_{S}=4\pi Dac_{\infty}$ (see also Fig. 1 left panel)†††These results can also be obtained by a simple rescaling argument. Setting $b=a/(1-\epsilon)$ into $k_{S}=4\pi Dbc_{\infty}$ in Eq. (23), one immediately recovers expressions (24) from $k/(4\pi Dac_{\infty})=k/[4\pi Dbc_{\infty}(1-\epsilon)]$..

These results have a very simple interpretation. Indeed, they correspond to the Smoluchowski rates into the respective equivalent spheres $S_{R_{e}}$, that is, spheres with the same surface (to the first order in $\epsilon$). From Eqs. (20) it can be seen immediately that $R_{e}=R\sqrt{1-4\epsilon/3}$ (prolate) and $R_{e}=R\sqrt{1-2\epsilon/3}$ (oblate). This suggests that in general the diffusive flux into a non-spherical boundary can be approximated by the Smoluchowski rate corresponding to its equivalent sphere to the first order in the difference of the relevant linear dimensions of the two surfaces. This agrees with the results of previous calculations by Berezhkovskii and Barzykin 52.

Our analysis shows that perturbing a sphere into a spheroid with given aspect ratio (to the first order in $\epsilon$) does not result in symmetric corrections to the Smoluchowski rate of the sphere, depending on whether the latter is compressed in the equatorial plane or elongated along the polar axis. More precisely, reducing a spherical surface to a prolate spheroid implies losing twice as much flux than reducing it to an oblate one. Conversely, increasing a spherical surface to a prolate spheroid implies gaining half as much flux than increasing it to an oblate one.

3.3 Partially absorbing spheroids

While the diffusive flux onto partially absorbing prolate and oblate spheroids can be obtained by solving infinite-dimensional systems of linear equations (see Appendices), its exact closed-form expressions are not available. However, it is possible to obtain approximate analytical expressions that turn out to be remarkably accurate for a wide range of aspect ratios and reactivities. By truncating the infinite set of linear equations to the zeroth order (see Appendices for the detailed calculations), we obtain in the case of compressed spheroids

[TABLE]

where $k_{S}=4\pi Dbc_{\infty}$ and $a_{E}=\sqrt{b^{2}-a^{2}}=b\sqrt{2\epsilon-\epsilon^{2}}$. These expressions are compared in Fig. 3 to the exact results (numerical solutions of the linear system to a fixed accuracy, cf Appendices) and to the perturbative expressions (23). One can see that the zeroth-order approximations (25) are barely distinguishable from the exact solutions over a very broad range of aspect ratios. More precisely, these expressions still yield accurate predictions of the total flux even for $a\ll b$, i.e. very elongated prolate spheroids and very flat oblate spheroids. Moreover, they are valid over the whole range of reactivities, from perfectly absorbing surfaces ($\kappa\to\infty$) to nearly reflecting boundaries ($\kappa\to 0$).

Figure 3 also illustrates clearly the range of validity of our perturbative calculations. The simple first-order perturbative corrections that one may compute easily from Eq. (15) provides a remarkably reliable estimate of the rate for values of $\epsilon=1-a/b$ as large as 0.5. In the case of oblate spheroids, Fig. 3 shows that the perturbative solution is an utterly reasonable approximation even for disk-like spheroids ($a=0$). Conversely, the flux onto thinner and thinner prolate spheroids become smaller and smaller, and consequently the perturbative solution fails in the limit $a\to 0$ (right-most upper panel in Fig. 3).

3.4 Beyond spheroidal shapes

In order to assess the accuracy of our perturbative approximation beyond spheroidal shapes, we solved the boundary value problem (8) through a finite-element method implemented in Matlab PDEtool (see Appendix C for details). We then compute numerically the total flux $k$ in the case of perfectly absorbing boundary conditions ($h=\kappa R/D\to\infty$), which is compared to our first-order approximation

[TABLE]

where $B_{0}^{\infty}$ is determined uniquely by the shape of the surface via Eq. (16). Even though this numerical scheme is applicable to any axisymmetric surface determined by a given function $r(\theta)$, we focus on two families of surfaces. The first family is determined by Legendre polynomials $P_{n}(z)$,

[TABLE]

Fig. 4 (top) illustrates the three surfaces corresponding to $n=2,3,4$ that we used in this study. Varying $\varepsilon$ from [math] to $1$, one progressively transforms the sphere of radius $R$ into three distinct shapes. A practical advantage of the choice of Eq. (27) is that the first-order correction to the total flux vanishes, as $B_{0}^{\infty}=0$ according to Eq. (16)) due to the orthogonality of Legendre polynomials. In other words, one has the somewhat surprising result that the total flux onto any surface of the form (27), with any $n>0$ and any $0\leq\varepsilon\leq 1$, is equal to $k_{S}$, up to the second-order term $\mathcal{O}(\varepsilon^{2})$, and thus $k_{\rm pert}=k_{S}$.

Fig. 4 (bottom) shows the difference between the total flux as computed numerically and its perturbative approximation $k_{\rm pert}$, for three surfaces of the kind (27) with $n=2,3,4$, as a function of $\varepsilon$. One can see that $(k-k_{\rm pert})/k_{S}$ vanishes as $\varepsilon^{2}$ in the limit $\varepsilon\to 0$, in agreement with our perturbative analysis. Remarkably, even for $\varepsilon=1$, which is far beyond the perturbation range, the perturbative prediction is within $30\%$ of the exact value.

In order to validate the perturbative approximation for shapes yielding nontrivial first-order corrections, we considered a second family of surfaces defined as

[TABLE]

For these surfaces, one has $B_{0}^{\infty}=1/(2n+1)$, so that the non-zero first-order correction reads

[TABLE]

Fig. 5 shows how the total flux to three surfaces of the kind (28) corresponding to $n=2,3,4$ deviates from our first-order perturbative approximation $k_{\rm pert}(n)$, eq. (29). Analogously to what reported in Fig. 4, this deviation grows as $\varepsilon^{2}$, confirming that the first-order term was evaluated correctly. Remarkably, the deviation at $\varepsilon=1$ is even smaller for the class of surfaces (28), as small as $\simeq 10\%$. This lower error (as compared to Fig. 4) reflects the fact that the second class of surfaces defined by Eq. (28) are somewhat rounder.

4 Conclusions and perspectives

In this paper we have tackled the problem of calculating the total stationary diffusive flux onto partially absorbing surfaces of revolution. We have outlined a general perturbative procedure that holds when the surface of revolution can be described as an axisymmetric perturbation of a sphere controlled by a single length scale (see Eq. (4)). We worked out an analytical expression for the reaction rate that agrees surprisingly well with the exact solutions for prolate and oblate spheroids, even in the case of rather non-spherical shapes with aspect ratio up to $0.5$. The exact and approximate formulas for the reaction rate for partially reactive prolate and oblate spheroids derived in this paper present their own interest for applications. Furthermore, we have tested our perturbative correction by comparing it to the exact values of the total flux to different families of non-spherical surfaces beyond spheroids. These comparisons show that our formula works surprisingly well even far from the perturbative limit, $\varepsilon=\mathcal{O}(1)$. More generally, our perturbative formula (15) for the case of a perfectly absorbing non-spherical boundary has a simple and general interpretation. The total flux to the first order in $\epsilon$ is the Smoluchowski rate into an equivalent sphere $S_{\rm eq}$, whose surface area equals that of the non-spherical surface of revolution $S$. Recalling the definitions (4) and (16), one has

[TABLE]

Comparing this result with the perturbative formula (17) for the total flux $k$ into the non-spherical sink $S$, one can conclude that

[TABLE]

where $R_{\rm eq}=\sqrt{S/4\pi}$. We note that the equivalence expressed by Eq. (31) proved in this paper had been conjectured in Ref. 52 based on intuitive arguments. However, our results prove that the equivalence (to $\mathcal{O}(\epsilon)$ at least) between a non-spherical surface and a sphere is only limited to perfectly absorbing boundary conditions (see again formula (15)). Furthermore, our method can in principle be employed to investigate whether formula (31) is valid beyond the first order in $\epsilon$. Moreover, the perturbative treatment presented in this paper can be also generalized to the case of non-uniform intrinsic reactivity $\kappa(\theta)$ (see 37), which constitutes an interesting extension for treating diffusion to non-spherical and non-uniform reactive surfaces.

As a final remark, we observe that the mathematics of diffusion in the presence of non-spherical reactive boundaries is important for the description of phoretic locomotion. Relevant examples include phoretic motion of non-spherical particles whose surface is only partially reactive (e.g. catalytic), the rest being chemically inactive 54, 53, or autophoresis purely induced by geometric asymmetries 55. Interestingly, in Ref.55 the authors perform a perturbative calculation of the autophoretic locomotion of non-spherical particles that includes computing the stationary concentration of a solute in the presence of partially absorbing boundary conditions. Their calculation is formally identical to our scheme, leading to an expression equivalent to our Eq. (14), even if their result is cast in the form of an expansion in series of Legendre polynomials. However, the focus of Ref.55 is not the reaction rate itself, which is the main subject of our analysis, rather it is the resulting phoretic velocity of the combined diffusion and hydrodynamic problem.

Appendix A Solution for partially absorbing prolate spheroids

Following Ref. 46, we write the concentration outside a prolate spheroid of the form $(a,a,b)$ (with semi-axes $a\leq b$) in prolate spheroidal coordinates $(\alpha,\theta,\phi)$ as

[TABLE]

with unknown coefficients $A_{n}$. Here $P_{n}(z)$ and $Q_{n}(z)$ are Legendre polynomials and functions of the first and second kind. The Robin boundary condition (8b) reads

[TABLE]

which should be satisfied by any $\theta\in(0,\pi)$. Here $\partial_{n}=\nabla\cdot\hat{n}$ is the normal derivative (directed inward), $\Lambda=D/\kappa$ is the reaction length 44, 20, 45, 56, 57 and $\alpha_{0}$ determines the surface of the prolate spheroid, $\cosh\alpha_{0}=b/a_{E}$, where $a_{E}=\sqrt{b^{2}-a^{2}}$. Multiplying this condition by $P_{m}(\cos\theta)\sin\theta$ and integrating over $\theta$ from [math] to $\pi$, one gets the infinite-dimensional system of linear equations on coefficients $A_{n}$

[TABLE]

where

[TABLE]

The total diffusive flux onto the prolate spheroid (i.e., the reaction rate) was computed in Ref. 46 as

[TABLE]

where $c_{\infty}$ is the concentration at infinity. For a perfectly reactive spheroid ($\Lambda=0$), the sum in Eq. (A) vanishes, and one gets a simple solution for the coefficients $A_{m}$:

[TABLE]

From Eqs. (36, 37), the exact form of the total flux onto the perfectly reactive prolate spheroid is recovered 58, 52

[TABLE]

where $k_{S}=4\pi Dc_{\infty}b$ is the Smoluchowski rate on the sphere of radius $b$.

In general, however, one needs to truncate the system (A) and then solve it numerically. Note that the same solution scheme is valid for Robin boundary condition (A) with a given function $g(\theta)$ on the left-hand side instead of $1$. The only difference is that $\delta_{m,0}$ on the left-hand side of Eq. (A) would be replaced by the integral of $g(\theta)$ with the Legendre polynomial $P_{m}(\cos\theta)$.

A.1 Numerical solution

For a numerical solution, one needs to truncate the system (A) to some order $n_{\rm max}$, to compute the matrix elements $F_{m,n}(\cosh\alpha_{0})$ for $m,n=0,1,\ldots,n_{\rm max}$, and then to invert the resulting matrix of size $(n_{\rm max}+1)\times(n_{\rm max}+1)$ numerically. The matrix elements $F_{m,n}(z)$ can be computed exactly by using the properties of Legendre polynomials. First, using the identity

[TABLE]

with

[TABLE]

and $a_{k}=\frac{\Gamma(k+1/2)}{\sqrt{\pi}\Gamma(k+1)}$ (with $a_{0}=1$), we get for even $m+n$

[TABLE]

where

[TABLE]

while $F_{m,n}(z)=0$ for odd $m+n$ due to the symmetry of Legendre polynomials. Second, the above integral can be computed via recursive relations:

[TABLE]

with

[TABLE]

As $F_{m,n}$ vanishes for odd $m+n$, the linear equations in the system (A) decouple into two subsystems, one for the coefficients $A_{2m}$ and the other for the coefficients $A_{2m+1}$. Moreover, as $\delta_{m,0}$ in the left-hand side of the system is not zero only for $m=0$, all the odd coefficients $A_{2m+1}$ are zero. In turn, the system for even coefficients $A_{2m}$ can be written as

[TABLE]

Truncating the sum to $n=n_{\rm max}$, computing $F_{2m,2n}(\cosh\alpha_{0})$ and solving numerically the truncated system with $n_{\rm max}+1$ equations onto the coefficients $A_{0},A_{2},\ldots,A_{2n_{\rm max}}$, one computes the reaction rate $k$ from Eq. (36).

A.2 Approximate solution

Solving the truncated system (A.1) for several values of $n_{\rm max}$, we checked that the coefficient $A_{0}$, determining the reaction rate, can be accurately computed from low-order truncations. In other words, the estimated $A_{0}$ converges very rapidly to its limit as $n_{\rm max}$ increases. As illustrated in Fig. 3, even the lowest-order truncation with $n_{\rm max}=0$ is accurate over the whole range of reactiony lengths $\Lambda$ even when the minor semi-axis $a$ is as small as $b/10$. This zeroth-order approximation yields

[TABLE]

where $k_{S}=4\pi Dc_{\infty}b$ is the Smoluchowski rate for a sphere of radius $b$. For a perfectly reactive spheroid, Eq. (45) at $\Lambda=0$ yields the exact result (38). As expected, one retrieves the Collins-Kimball result for a partially reactive sphere in the limit $a\to b$:

[TABLE]

The lowest-order approximation becomes less accurate in the limit $a\to 0$ when prolate spheroids are getting closer to the shape of a needle. In this case, one can consider the first-order truncation with $n_{\rm max}=1$, in which case the truncated $2\times 2$ system of equations can be again solved explicitly, yielding a cumbersome but remarkably accurate solution.

Figure 3 illustrates the behavior of the scaled reaction rate $k/k_{S}$ for a partially reactive prolate spheroid as a function of $D/(\kappa b)=\Lambda/b$. For all considered values of $\Lambda/b$ and $a/b$, the numerical solution of the truncated system of linear equations converges very rapidly with the truncated order $n_{\rm max}$. In particular, one can see that solutions with $n_{\rm max}=5$ (referred to as “exact solution”) and $n_{\rm max}=1$ are barely distinguishable over all $\Lambda/b$ even for a very narrow spheroid with $a/b=0.01$. This observation confirms that this explicit approximate solution is very accurate. A simpler zeroth-order approximate solution from Eq. (45) is also very accurate, in spite of small deviations at large reactivity (small $\Lambda/b$). Finally, the perturbative solution (55) is accurate for small and moderate $\epsilon=1-a/b$, and fails only when $\epsilon$ becomes close to $1$ (i.e., $a$ close to [math]).

A.3 Perturbative solution

An approximate solution of the above system of equations can be obtained in the limit when the prolate spheroid approaches a sphere, i.e., $a=b(1-\epsilon)$, so that $a_{E}^{2}=2b^{2}\epsilon+\mathcal{O}(\epsilon^{2})$, and thus $\cosh\alpha_{0}=1/\sqrt{2\epsilon}\to\infty$. In the limit $z\to\infty$, we have thus

[TABLE]

with

[TABLE]

(this integral can be evaluated explicitly). Denoting

[TABLE]

we rewrite the system (A.1) up to the first order in $\epsilon$ as

[TABLE]

where

[TABLE]

One gets thus from Eq. (50) with $m=0$:

[TABLE]

whereas $a_{2}=\mathcal{O}(\epsilon)$ so that the second term in the numerator can be neglected, yielding (with $f_{0,0}=1/3$)

[TABLE]

Given that $Q_{0}(z)=\frac{1}{2}\ln\frac{z+1}{z-1}$ and $Q^{\prime}_{0}(z)=\frac{1}{1-z^{2}}$, we obtain $e_{0}=-1/(1+h)$, where $h=b/\Lambda=\kappa b/D$. We get thus

[TABLE]

from which

[TABLE]

This expression agrees with our general formula (15) for $f(\theta)=-\sin^{2}\theta$ that approximately describes a prolate spheroid for small $\epsilon$.

Appendix B Oblate spheroid

As computations for an oblate spheroid of the form $(b,b,a)$ (with semi-axes $a\leq b$) are similar, we only sketch the main steps. In the oblate spheroidal coordinates $(\alpha,\theta,\phi)$, an axiosymmetric solution of the Laplace equation outside the spheroid reads

[TABLE]

where the coefficients $A_{n}$ are determined from the Robin boundary condition

[TABLE]

from which

[TABLE]

Note that the boundary of the oblate spheroid is defined again as $\cosh\alpha_{0}=b/a_{E}$, with $a_{E}=\sqrt{b^{2}-a^{2}}$. The total diffusive flux reads as (see Ref. 47 for technical details)‡‡‡ We note that Appendix C of Ref. 47 contains several technical misprints that do not alter the final results.

[TABLE]

As previously, only the even coefficients matter so that

[TABLE]

For a perfectly reactive oblate spheroid ($\Lambda=0$), the sum vanishes, so that $A_{0}=1/Q_{0}(i\sinh\alpha_{0})$, from which Eq. (58) helps to retrieve the exact formula for the reaction rate 58, 52

[TABLE]

where we used $Q_{0}(i\sinh\alpha_{0})=-i\arcsin(1/\cosh\alpha_{0})$.

For a partially reactive spheroid, the zeroth-order truncation of the system (B) yields the following approximation for the reaction rate

[TABLE]

In the limit $a\to b$, one retrieves the Collins-Kimball result (46) for a partially reactive sphere. In the opposite limit $a\to 0$, one finds the reaction rate for a partially absorbing disk of radius $b$:

[TABLE]

The accuracy of this approximation is illustrated in Fig. 3. An even higher accuracy can be achieved by using the truncation order $n_{\rm max}=1$, for which the resulting $2\times 2$ matrix can be solved explicitly. In comparison to Eq. (45) for the prolate spheroid, one can note that the geometric aspect is fully disentangled from the reactive aspect, which is accounted via the factor $1/(1+\Lambda/b)$, as for the sphere.

Skipping technical detials, we also provide the perturbative solution of Eqs. (B) up to the first order in $\epsilon$:

[TABLE]

This solution agrees with our general formula (15) for $f(\theta)=-\cos^{2}\theta$ that approximately describes an oblate spheroid for small $\epsilon$.

Appendix C Validation by a finite-element method

We seek a numerical solution of the boundary value problem (8) in the limit $h\to\infty$. For this purpose, we use a finite-element method implemented in Matlab PDEtool. As this method meshes the computational domain, the latter should be bounded by an artificially imposed outer boundary $\Gamma_{\alpha}$. For convenience, we choose $\Gamma_{\alpha}$ to be a sphere of radius $R_{\alpha}=R/\alpha$ centered at the origin: $\Gamma_{\alpha}=\{\vec{r}\in{\mathbb{R}}^{3}\leavevmode\nobreak\ :\leavevmode\nobreak\ \|\vec{r}\|=R_{\alpha}\}$. We aim thus to solve numerically the modified problem

[TABLE]

in the new bounded domain: $\Omega_{\alpha}=({\mathbb{R}}^{3}\setminus T^{\prime})\cap\{\vec{r}\in{\mathbb{R}}^{3}\leavevmode\nobreak\ :\leavevmode\nobreak\ \|\vec{r}\|<R_{\alpha}\}$. As $\alpha$ goes to zero, the outer boundary $\Gamma_{\alpha}$ moves away from the target $T^{\prime}$, $\Omega_{\alpha}$ blows up and approches $({\mathbb{R}}^{3}\setminus T^{\prime})$, and $u_{\alpha}$ converges to the solution $u$ of the original problem (8). Even though the boundary condition at the imposed outer boundary $\Gamma_{\alpha}$ does not matter in the limit $\alpha\to 0$, we choose the Robin condition (64c) to speed up the convergence. In fact, if the surface ${\partial\Omega}^{\prime}$ was a sphere of radius $1$, then the solution of Eqs. (8) would simply be $1/r$. Moreover, for an arbitrary surface ${\partial\Omega}^{\prime}$, the solution of Eqs. (8) is expected to asymptotically decay as $\gamma/r$ (with some constant $\gamma$), and the Robin condition (64c) is again consistent with this behavior. In summary, even so the condition (64c) could be replaced by more common Dirichlet or Neumann conditions, the choice (64c) is preferred.

The axial symmetry of the domain $\Omega_{\alpha}$ reduces the original three-dimensional problem to a two-dimensional one. The Laplace equation (64a) can be written in a standard matrix form in spherical coordinates $(r,\theta,\phi)$ as

[TABLE]

where $c$ is the diagonal $2\times 2$ matrix with the diagonal elements equal to $r^{2}\sin\theta$ and $\sin\theta$. The two-dimensional computational domain is then:

[TABLE]

i.e., it is a kind of rectangle of width $\pi$, whose bottom “edge” is curved and determined by a given function $r(\theta)$. The boundary conditions (64b, 64c) are set on the bottom and top “edges”, whereas the Neumann condition $\partial_{\theta}u_{\alpha}=0$ is imposed on the left and right edges.

Once the problem (64) is solved numerically, the total diffusive flux follows as

[TABLE]

Alternatively, the flux conservation allows one to compute the total flux by integrating over the outer boundary $\Gamma_{\alpha}$ and using Eq. (64c),

[TABLE]

which provides a much simpler way to compute the total flux. As $k_{\alpha}$ approaches $k$ of the original problem in the limit $\alpha\to 0$, this expression can be used as an approximation of $k$ for $\alpha$ small enough.

For each surface determined by a given function $r(\theta)$, we computed the total flux $k_{\alpha}$ by setting $R_{\alpha}=10R$ and the maximal mesh size $\eta_{\rm max}$ of $0.05$. To check the accuracy of numerical computations, $k_{\alpha}$ was also computed by using either $R_{\alpha}=20R$ (farther outer boundary), or $\eta_{\rm max}=0.025$ (finer mesh). For all considered surfaces, these changes did not almost affect $k_{\alpha}$ (the relative error being much smaller than $1\%$). This observation confirms the high accuracy of our numerical computation.

Acknowledgements

FP would like to thank Stefano Angioletti-Uberti for enlightening discussions.

Notes and references

• 1 G. Palazzo and D. Berti, in Colloidal Foundations of Nanoscience, 199–231, 2014.
• 2 S. Rice, Diffusion-Limited Reactions, Elsevier, Amsterdam, 1985.
• 3 N. Welsch, A. Wittemann, and M. Ballauff, J. Phys. Chem. B, 2009, 113, 16039–16045.
• 4 M. Galanti, D. Fanelli, S. Angioletti-Uberti, M. Ballauff, J. Dzubiella, and F. Piazza, Phys. Chem. Chem. Phys., 2016, 18, 20758–20767.
• 5 R. Roa, S. Angioletti-Uberti, Y. Lu, J. Dzubiella, F. Piazza, and M. Ballauff, Z. Phys. Chem., 2018, 232, 773–803.
• 6 I. V. Gopich and A. Szabo, Proc. Nat. Acad. Sci. USA, 2013, 110, 19784–19789.
• 7 L. Bongini, D. Fanelli, F. Piazza, P. De Los Rios, M. Sanner, and U. Skoglund, Phys. Biol., 2007, 4, 172–180.
• 8 F. Piazza, P. De Los Rios, D. Fanelli, L. Bongini, and U. Skoglund, Eur. Biophys. J., 2005, 34, 899–911.
• 9 H. C. Berg and E. M. Purcell, Biophys. J., 1977, 20, 193-239.
• 10 D. Shoup, G. Lipari, and A. Szabo, Biophys. J., 1981, 36, 697-714.
• 11 D. Shoup and A. Szabo, Biophys. J., 1982, 40, 33–39.
• 12 R. Zwanzig, Proc. Natl. Acad. Sci. USA, 1990, 87, 5856.
• 13 R. Zwanzig and A. Szabo, Biophys. J., 1991, 60, 671-678 (1991).
• 14 D. A. Lauffenburger and J. Linderman, Receptors: Models for Binding, Trafficking, and Signaling, Oxford University Press, 1993.
• 15 M. Tachiya and S. D. Traytak, J. Phys.: Conden. Matter, 2007, 19, 065111.
• 16 H. C. Berg, Random walks in biology, Princeton University Press, 1993.
• 17 L. Karp-Boss, E. Boss, P. A. Jumars, et al., Oceanogr. Mar. Biol., 1996, 34, 71–108.
• 18 T. Kiørboe, A mechanistic approach to plankton ecology, Princeton University Press, 2008.
• 19 A. Sozza, F. Piazza, M. Cencini, F. De Lillo, and G. Boffetta, Phys. Rev. E, 2018, 97, 023301.
• 20 B. Sapoval, M. Filoche, and E. Weibel, Proc. Nat. Ac. Sci. USA, 2002, 99, 10411-10416.
• 21 D. S. Grebenkov, M. Filoche, B. Sapoval, and M. Felici, Phys. Rev. Lett., 2005, 94, 050602.
• 22 A. S. Serov, C. Salafia, D. S. Grebenkov, and M. Filoche, J. Appl. Physiol., 2016, 120, 17-28.
• 23 F. Piazza, G. Foffi, and C. De Michele, J. Phys.: Conden. Matter, 2013, 25, 245101.
• 24 M. Galanti, D. Fanelli, S. D. Traytak, and F. Piazza, Phys. Chem. Chem. Phys., 2016, 18, 15950-15954.
• 25 D. S. Grebenkov and S. D. Traytak, J. Comput. Phys., 2019, 379, 91-117.
• 26 M. von Smoluchowski, Physik Z., 1916, 17, 557–571.
• 27 M. von Smoluchowski, Z. Phys. Chem., 1917, 92, 129–168.
• 28 D. S. Grebenkov, in Chemical Kinetics: Beyond the Textbook, Eds. K. Lindenberg, R. Metzler, and G. Oshanin World Scientific, 2019; DOI: 10.1142/q0209; available online as ArXiv: 1806.11471.
• 29 G. H. Weiss, J. Stat. Phys., 1986, 42, 3-36.
• 30 A. Berezhkovskii, Y. Makhnovskii, M. Monine, V. Zitserman, and S. Shvartsman, J. Chem. Phys., 2004, 121, 11390.
• 31 A. M. Berezhkovskii, M. I. Monine, C. B. Muratov, and S. Y. Shvartsman, J. Chem. Phys., 2006, 124, 036103.
• 32 C. Muratov and S. Shvartsman, Multiscale Model. Simul., 2008, 7, 44-61.
• 33 L. Dagdug, M. Vázquez, A. Berezhkovskii, and V. Zitserman, J. Chem. Phys., 2016, 145, 214101.
• 34 A. E. Lindsay, A. J. Bernoff, and M. J. Ward, Multiscale Model. Simul., 2017, 15, 74-109.
• 35 A. J. Bernoff and A. E. Lindsay, SIAM J. Appl. Math., 2018, 78, 266-290.
• 36 A. Bernoff, A. Lindsay, and D. Schmidt, Multiscale Model. Simul., 2018, 16, 1411-1447.
• 37 D. S. Grebenkov, J. Chem. Phys., 151, 104108.
• 38 O. Bénichou, M. Moreau, and G. Oshanin, Phys. Rev. E, 2000, 61, 3388-3406.
• 39 J. Reingruber and D. Holcman, Phys. Rev. Lett., 2009, 103, 148102.
• 40 S. D. Lawley and J. P. Keener, SIAM J. Appl. Dyn. Sys., 2015, 14, 1845-1867.
• 41 P. C. Bressloff, J. Phys. A., 2017, 50, 133001.
• 42 H. Sano and M. Tachiya, J. Chem. Phys., 1979, 71, 1276-1282.
• 43 H. Sano and M. Tachiya, J. Chem. Phys., 1981, 75, 2870-2878.
• 44 B. Sapoval, Phys. Rev. Lett., 1994, 73, 3314-3317.
• 45 D. S. Grebenkov, M. Filoche, and B. Sapoval, Phys. Rev. E, 2006, 73, 021103.
• 46 S. D. Traytak and D. S. Grebenkov, J. Chem. Phys., 2018, 148, 024107.
• 47 D. S. Grebenkov and D. Krapf, J. Chem. Phys., 2018, 149, 064117.
• 48 F. C. Collins and G. E. Kimball, J. Coll. Sci., 1949, 4, 425–437.
• 49 M. Filoche and B. Sapoval, Eur. Phys. J. B, 1999, 9, 755-763.
• 50 D. S. Grebenkov, M. Filoche, and B. Sapoval, Eur. Phys. J. B, 2003, 36, 221-231.
• 51 D. S. Grebenkov, in Focus on Probability Theory, Ed. L. R. Velle, pp. 135-169, Nova Science Publishers, 2006.
• 52 A. M. Berezhkovskii and A. V. Barzykin, J. Chem. Phys., 2007, 126, 106102.
• 53 M. N. Popescu, S. Dietrich, M. Tasinkevych, and J. Ralston, Eur. Phys. J. E, 2010, 31(4), 351-367.
• 54 W. E. Uspal, M. N. Popescu, M. Tasinkevych, and S. Dietrich, New J. Phys., 2018, 20(1), 015013.
• 55 S. Michelin and E. Lauga, Eur. Phys. J. E, 2015, 38(7), i2015-15007-6.
• 56 D. S. Grebenkov, Fractals, 2006, 14, 231-243.
• 57 D. S. Grebenkov, Phys. Rev. E, 2015, 91, 052108.
• 58 L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media, Pergamon, Oxford, 1984.