The narrow escape problem in a circular domain with radial piecewise constant diffusivity

TL;DR

This paper investigates how spatially varying diffusivity in a circular domain affects the mean first passage time for particles to escape through a small boundary region, revealing optimal configurations and discontinuous dependencies.

Contribution

It introduces a model with piecewise constant diffusivity in a circular domain to analyze the narrow escape problem, highlighting the effects of inhomogeneous diffusion on MFPT.

Findings

MFPT can be minimized by tuning the outer shell width in annulus geometry.

In two-shell geometry, MFPT depends monotonically on model parameters.

The optimal starting point distance depends discontinuously on diffusivity ratio.

Abstract

The stochastic motion of particles in living cells is often spatially inhomogeneous with a higher effective diffusivity in a region close to the cell boundary due to active transport along actin filaments. As a first step to understand the consequence of the existence of two compartments with different diffusion constant for stochastic search problems we consider here a Brownian particle in a circular domain with different diffusion constants in the inner and the outer shell. We focus on the narrow escape problem and compute the mean first passage time (MFPT) for Brownian particles starting at some pre-defined position to find a small region on the outer reflecting boundary. For the annulus geometry we find that the MFPT can be minimized for a specific value of the width of the outer shell. In contrast for the two-shell geometry we show that the MFPT depends monotonously on all model parameters, in particular on the outer shell width. Moreover we find that the distance between the starting point and the narrow escape region which maximizes the MFPT depends discontinuously on the ratio between inner and outer diffusivity.

Equations

$$t(\mathbf{x})={\int }_{\Omega }d\mathbf{y}{\int }_0^{\infty }d\tau p(\mathbf{y},\tau |\mathbf{x}).$$

$$D_0{\nabla }^{2}t(\mathbf{x})=-1,\mathbf{x}\in \Omega$$

$$t(\mathbf{x})=0,\mathbf{x}\in \partial {\Omega }_{\epsilon }$$

$$\mathbf{n}\cdot \nabla t(\mathbf{x})=0,\mathbf{x}\in \partial \Omega \backslash \partial {\Omega }_{\epsilon }.$$

$$t(\mathbf{x})\simeq ⟨t⟩\simeq -\frac{|\Omega |}{\pi D_0}\ln \epsilon .$$

$$t(\mathbf{x})\simeq \frac{|\Omega |}{\pi D_0}\ln \frac{|\mathbf{x}-{\mathbf{x}}_0|}{\epsilon \mathcal{R}},$$

$${\nabla }^{2}G(\mathbf{x}|{\mathbf{x}}_0)=\frac{1}{|\Omega |},\mathbf{x}\in \Omega$$

$$\mathbf{n}\cdot \nabla G(\mathbf{x}|{\mathbf{x}}_0)=0,\mathbf{x}\in \partial \Omega$$

$$G(\mathbf{x}\to {\mathbf{x}}_0|{\mathbf{x}}_0)=-\frac{1}{\pi }\ln \frac{|\mathbf{x}-{\mathbf{x}}_0|}{\mathcal{R}}+R({\mathbf{x}}_0|{\mathbf{x}}_0),$$

$${\int }_{\Omega }d\mathbf{x}G(\mathbf{x}|{\mathbf{x}}_0)=0,$$

$$⟨t⟩=\frac{|\Omega |}{\pi D_0}\left[-\ln \frac{\epsilon }{4}+\pi R({\mathbf{x}}_0|{\mathbf{x}}_0)\right],$$

$$t(\mathbf{x})=⟨t⟩-\frac{|\Omega |}{D_0}G(\mathbf{x}|{\mathbf{x}}_0).$$

$$R(\mathbf{x}|{\mathbf{x}}_0)=G(\mathbf{x}|{\mathbf{x}}_0)+\frac{1}{\pi }\ln \frac{|\mathbf{x}-{\mathbf{x}}_0|}{\mathcal{R}}.$$

$$t(\mathbf{x})=\frac{|\Omega |}{\pi D_0}\left[\ln \frac{4|\mathbf{x}-{\mathbf{x}}_0|}{\epsilon \mathcal{R}}+\pi R({\mathbf{x}}_0|{\mathbf{x}}_0)-\pi R(\mathbf{x}|{\mathbf{x}}_0)\right].$$

$$⟨t⟩\underset{\epsilon \ll 1}{\simeq }\frac{{\mathcal{R}}^2}{D_0}\left[-\ln \frac{\epsilon }{4}+\frac{1}{8}\right],$$

$$t(\mathbf{x})\underset{\epsilon \ll 1}{\simeq }\frac{{\mathcal{R}}^2}{D_0}\left[\ln \frac{4|\mathbf{x}-{\mathbf{x}}_0|}{\epsilon \mathcal{R}}+\frac{{\mathcal{R}}^2-|\mathbf{x}{|}^2}{4{\mathcal{R}}^2}\right].$$

$$t(0)\underset{\epsilon \ll 1}{\simeq }\frac{{\mathcal{R}}^2}{D_0}\left[-\ln \frac{\epsilon }{4}+\frac{1}{4}\right],$$

$$t_{\max }\underset{\epsilon \ll 1}{\simeq }t(-{\mathbf{x}}_0)\underset{\epsilon \ll 1}{\simeq }\frac{{\mathcal{R}}^2}{D_0}\left[-\ln \frac{\epsilon }{4}+\ln 2\right],$$

$$\frac{D_{0}t(z)}{{\mathcal{R}}^2}=\frac{1-|z{|}^2}{4}+\ln \left|\frac{1-z+\sqrt{[1-z\exp (-i\epsilon /2)][1-z\exp (i\epsilon /2)]}}{2\sin (\epsilon /4)}\right|$$

$$D_{\Delta }{\nabla }^{2}t(\mathbf{x})=-1,\mathbf{x}\in \Omega$$

$$t(\mathbf{x})=0,\mathbf{x}\in \partial {\Omega }_{\epsilon }$$

$$\mathbf{n}\cdot \nabla t(\mathbf{x})=0,\mathbf{x}\in \partial \Omega \backslash \partial {\Omega }_{\epsilon }$$

$$\mathbf{n}\cdot \nabla t(\mathbf{x})=0,\mathbf{x}\in \partial {\Omega }_{\Delta }.$$

$$R({\mathbf{x}}_0|{\mathbf{x}}_0)=\frac{3}{8\pi }-\frac{1}{4\pi \delta (2-\delta )}-\frac{(1-\delta {)}^4\ln (1-\delta )}{2\pi {\delta }^2(2-\delta {)}^2}+\frac{2}{\pi }\sum _{n=1}^{\infty }\frac{(1-\delta {)}^{2n}}{n\left[1-(1-\delta {)}^{2n}\right]},$$

$$\frac{D_{\Delta }⟨t⟩}{{\mathcal{R}}^2}\underset{\epsilon \ll 1}{\simeq }\delta (2-\delta )\left[-\ln \frac{\epsilon }{4}+\frac{3}{8}+2\sum _{n=1}^{\infty }\frac{(1-\delta {)}^{2n}}{n\left[1-(1-\delta {)}^{2n}\right]}\right]-\frac{1}{4}-\frac{(1-\delta {)}^4\ln (1-\delta )}{2\delta (2-\delta )}.$$

$$S(X)=\sum _{n=1}^{\infty }\frac{X^{2n}}{n(1-X^{2n})}.$$

$$\frac{D_{\Delta }⟨t⟩}{{\mathcal{R}}^2}\simeq -2\delta \ln \frac{\epsilon }{4}+\frac{{\pi }^2}{3}+\mathcal{O}(\delta ,\epsilon ).$$

$$\frac{D_{\Delta }⟨t⟩}{{\mathcal{R}}^2}\underset{\epsilon =2\pi }{=}\frac{3\delta (2-\delta )}{8}-\frac{1}{4}-\frac{(1-\delta {)}^4\ln (1-\delta )}{2\delta (2-\delta )}.$$

$$\frac{D_{\Delta }}{{\mathcal{R}}^2}t(\partial \Omega \to \partial {\Omega }_{\epsilon })\simeq \delta (2-\delta )\left[-\ln \frac{\epsilon }{4}+2\sum _{n=1}^{\infty }\frac{(1-\delta {)}^{2n}}{n\left[1-(1-\delta {)}^{2n}\right]}\right].$$

$$\frac{D_{\Delta }{t}_{\max }}{{\mathcal{R}}^2}\simeq \frac{D_{\Delta }t(-{\mathbf{x}}_0)}{{\mathcal{R}}^2}\underset{\epsilon \ll 1}{\simeq }\delta (2-\delta )\left[-\ln \frac{\epsilon }{8}+4\sum _{k=0}^{\infty }\frac{(1-\delta {)}^{2(2k+1)}}{(2k+1)\left[1-(1-\delta {)}^{2(2k+1)}\right]}\right].$$

$$\frac{D_{\Delta }{t}_{\max }}{{\mathcal{R}}^2}\simeq -2\delta \ln \frac{\epsilon }{8}+\frac{{\pi }^2}{2}+\mathcal{O}(\delta ,\epsilon ).$$

$$\frac{D_{\Delta }{t}_{\max }}{{\mathcal{R}}^2}\underset{\epsilon =2\pi }{=}\frac{\delta (2-\delta )}{4}+\frac{(1-\delta )\ln (1-\delta )}{2}$$

$$\beta U(\mathbf{x})=-g\ln \frac{|\mathbf{x}|}{\mathcal{R}}\Theta (\mathcal{R}-\Delta -|\mathbf{x}|)$$

$$D_{\Delta }{\nabla }^{2}t(\mathbf{x})-D_{\Delta }\beta \nabla U(\mathbf{x})\cdot \nabla t(\mathbf{x})=-1,\mathbf{x}\in \Omega$$

Full text

I Introduction

The narrow escape problem is ubiquitous in biology and chemistry problems Schuss et al. (2007); Bressloff and Newby (2013); Chou and D’Orsogna (2014); Holcman and Schuss (2014); Iyer-Biswas and Zilman (2015) and consists in calculating the mean first passage time (MFPT) of a random walk or a Brownian motion from a starting point to a small absorbing window, the so-called escape or target region, on an otherwise reflecting boundary of a bounded domain. Mathematically the problem reduces to solving a Poisson equation with mixed Neumann-Dirichlet boundary conditions Redner (2001). These equations can be solved numerically with a partial differential equation (PDE) solver for any complex geometry boundaries. Moreover, new numerical techniques allow the efficient simulation of the stochastic process using Kinetic Monte Carlo algorithms Schwarz and Rieger (2013); Schwarz et al. (2016a, b).

Analytically, the narrow escape problem has been widely studied in the last decades Holcman and Schuss (2004); Kolokolnikov et al. (2005); Singer et al. (2006a, b, c); Schuss et al. (2007) by focusing on the limiting case of small escape regions on the domain boundary. In two dimensions, the leading order of the MFPT is known to be proportional to the logarithm of harmonic measure of the escape region, which is proportional to its perimeter in the case of a disk-like domain Holcman and Schuss (2004); Grebenkov (2016). Also the sub-leading behavior of the MFPT has been studied which plays an important role for the logarithmic dependence and depends generally on the pseudo Green’s function of the domain Condamin et al. (2007); Bénichou and Voituriez (2008); Cheviakov et al. (2010); Pillay et al. (2010); Chevalier et al. (2011); Cheviakov et al. (2012). An exact solution of the escape problem, via a general integral formula, was recently shown for any simply connected planar domain with an arbitrary diffusivity implying the conformal mapping of the unit disk onto this domain Grebenkov (2016). Moreover, some asymptotically exact formulas was found for non-spherical three-dimensional domains Gomez and Cheviakov (2015).

Recently the narrow escape problem was studied in spatially inhomogeneous environments Schwarz et al. (2016a). The type of spatial inhomogeneity was inspired by the spatial organization of the cytoskeleton of cells with a centrosome Schwarz et al. (2016b); Hafner and Rieger (2016, 2018) along which ballistic transport is possible in addition to simple diffusion: a circular or spherical domain is divided into two concentric shells, the inner shell allowing only radial ballistic transport (and diffusion) and the outer shell of width $\Delta$ allows multi-directional ballistic transport (and diffusion). In these studies it was shown that the MFPT for the narrow escape problem is optimizable with respect to the width of the outer shell. The physical reason for the optimization of the MFPT in the considered setup was argued to be related to the accelerated effective diffusion constant in the outer shell Schwarz et al. (2016a, b); Hafner and Rieger (2016, 2018), in reminiscence of the optimization of the reaction time in the presence of surface-mediated diffusion Bénichou et al. (2010, 2011); Calandre et al. (2012); Rupprecht et al. (2012a, b) with different models of absorption and desorption.

In view of this reminiscence the natural question arises whether the MFPT also shows a non-monotonous - and hence optimizable - dependence on the thickness of the outer shell if ballistic transport is neglected and only the diffusion constant in the outer shell is larger than the diffusion constant in the inner shell. Therefore in this paper we consider this two-shell geometry with only a diffusive motion with two different diffusivities in the center of the cell and in the cortex of width $\Delta$. Diffusive motion in such discontinuous media, presenting a discontinuity of the coefficient of diffusion throughout an interface, is little studied. The one-dimensional profile of concentration of diffusive particles in such domain, presenting a jump close to the interface, was investigated LaBolle et al. (2000); Alvarez-Ramirez et al. (2014a, b). Moreover, a simulation scheme for one-dimensional problems is presented in the Itō formalism taking into account the drift term proportional to the gradient of the diffusivity Lejay and Martinez (2006); Lejay and Pichot (2012); Lejay and Maire (2013) and compatible with the Fick’s semi-empirical law. In fact, different writings of Langevin equation are possible in heterogeneous media, where the diffusivity is spatially varying Vaccario et al. (2015).

In this article, we will be interested in the first passage problem for the disk-like geometries of radius ${\cal R}$ presenting one target of angle $\varepsilon$ located on the outer boundary, which represents an opening of length $\varepsilon{\cal R}$. The domain is separated into two shells: the inner shell of radius ${\cal R}-\Delta$ playing the role of the center of the cell, where the diffusive particle has a diffusivity $D_{0}$ and the outer shell of width $\Delta$ playing the role of the cortex, where the diffusive particle has a diffusivity $D_{\Delta}$. Since our main motivation was to investigate an idealized case in which the outer shell represents the actin cortex of a cell and therefore has a higher effective diffusion constant we mainly, but not exclusively, look at $D_{\Delta}>D_{0}$. We will look at two different expansions within this circular discontinuous geometry: the narrow escape problem where the escape region is small ($\varepsilon\ll 1$) and the thin outer region problem ($\Delta\ll{\cal R}$). We present in section II a well-known solution of the narrow escape problem of the disk (same diffusivity in the two shells), with a geometry represented in the Fig. 1a. In section III, we study the annulus geometry, shown in Fig. 1b, where the diffusive particle is excluded from the inner shell. Finally, we will study in section IV the discontinuous problem whose the geometry is shown in Fig. 1c. The section V concludes with a discussion and an outlook.

II Narrow escape problem for the disk-like geometries

We first look at the escape problem for the two-dimensional disk-like geometries. This problem is equivalent to the two-shell geometry with the same diffusivity in both shells such that $D_{\Delta}=D_{0}$ which is independent of the value of $\Delta$. The available space for the Brownian particles is a disk of radius ${\cal R}$, denoted as $\Omega$, and its area is equal to $|\Omega|=\pi{\cal R}^{2}$. The external boundary is located at the radius $|{\bf x}|={\cal R}$ and is denoted as $\partial\Omega$ whereas the absorbing part of this boundary located at $|\theta|\leq\varepsilon/2$ is denoted as $\partial\Omega_{\varepsilon}$. The geometry of this domain is sketched in Fig. 1a. The MFPT, denoted as $t({\bf x})$ for a Brownian particle starting at the point ${\bf x}$, can be related to the probability density function $p({\bf y},\tau|{\bf x})$ to be at the point ${\bf y}$ after a time $\tau$ without having reached the escape region via

[TABLE]

The equations for the function $t({\bf x})$ are then deduced from the backward equations for the probability density function Redner (2001)

[TABLE]

The physical parameters of this problem are: the diffusivity of the Brownian particle $D_{0}$, the radius of the disk ${\cal R}$ and the angle $\varepsilon$ of the escape region. Only one pertinent parameter remains to study the narrow escape (NE) problem: $\varepsilon$, using dimensionless variables after rescaling the length by the radius ${\cal R}$ and the time by ${\cal R}^{2}/D_{0}$. The numerical solution of Eqs. (2-4) for the MFPT obtained with the PDE solver FreeFem++ Hecht (2013) is shown in Fig. 2a for the escape angle $\varepsilon=0.2$.

In this section, we will reproduce the main results of the literature Kolokolnikov et al. (2005); Singer et al. (2006b); Schuss et al. (2007); Pillay et al. (2010); Chevalier et al. (2011); Cheviakov et al. (2012); Caginalp and Chen (2012); Rupprecht et al. (2015) for the MFPT of the disk geometry, which will constitute the starting point of our analysis. The leading order of the MFPT in the NE limit ($\varepsilon\ll 1$) is given by the constant solution equal to the spatial average value of the MFPT Holcman and Schuss (2004)

[TABLE]

Applying the Benichou’s method Condamin et al. (2007); Bénichou and Voituriez (2008), one obtains straightforwardly the space-dependent solution in leading order

[TABLE]

where ${\bf x_{0}}=({\cal R},0)$ is the center of the escape region. However, the space average of the MFPT remains equal to Eq. (5) since the average value of $\ln(|{\bf x}-{\bf x_{0}}|/{\cal R})$ is zero over the disk. To obtain the order ${\cal O}(\varepsilon^{0})$ result for the MFPT, which is important even if $\varepsilon$ is as small as $10^{-2}$ due to the logarithmic behavior, we introduce the pseudo Green’s function $G({\bf x}|{\bf x_{0}})$, following Refs. Pillay et al. (2010); Chevalier et al. (2011), defined by the equations

[TABLE]

where $R({\bf x_{0}}|{\bf x_{0}})$ is the (unknown) regular part of the Green’s function at the center of the escape region ${\bf x_{0}}$. The spatial average of the MFPT, hereinafter called Global Mean First Passage Time (GMFPT), is given by

[TABLE]

and the expression of the MFPT $t({\bf x})$ is

[TABLE]

A general proof of these equations is shown in the appendix A. The equations (7-10) can be solved for the disk-like geometries we study here, for the regular part of the Green’s function defined by

[TABLE]

The MFPT function $t({\bf x})$ satisfies thus the expression

[TABLE]

For the disk geometry represented in Fig. 1, the regular part of the Green’s function satisfies the Eq. (75). Its value close to the escape region is equal to $R({\bf x_{0}}|{\bf x_{0}})=1/8\pi$. The GMFPT expression is

[TABLE]

and the expression of the MFPT for a Brownian particle starting at ${\bf x}$ is

[TABLE]

At the leading order, this solution is comparable to Benichou’s solution (6) close to the escape region. The expression of the MFPT for a Brownian particle starting at the center of the disk ${\bf x}={\bf 0}$, denoted hereinafter as CMFPT, is then

[TABLE]

and the expression of the maximal value of the MFPT, denoted hereinafter as MMFPT, is

[TABLE]

corresponding to the value of the MFPT for a particle starting at the position $-{\bf x_{0}}=(-{\cal R},0)$, at the maximal distance to the escape region. The GMFPT, CMFPT and MMFPT values are shown in Fig. 2b as a function of $\varepsilon$. The narrow escape expressions given by Eqs. (15), (17) and (18) are compatible with the numerical solutions of the PDEs (2-4) for $\varepsilon\lesssim 2\pi/3$.

These numerical solutions are obtained with the finite element method using the software package FreeFem++ Hecht (2013), with a relative interpolation error of order $10^{-4}$ for a mesh grid built with around $10^{4}$ vertices. For our study of the GMFPT, we can compare the numerically interpolated value of $\pi R({\bf x_{0}}|{\bf x_{0}})$ which is approximatively equal to $0.124728$ with the analytical solution $1/8$ which represents a relative error of $0.217\%$. In the following sections, we considerer the data generated by FreeFem++ as the exact numerical solution.

Furthermore, Singer et al. Singer et al. (2006b) solved directly the partial differential equations for the MFPT (2-4) and gave an exact expression of the MFPT for all values of $\varepsilon$, as an infinite series with implicit coefficients written as an integral. The explicit form of these coefficients was given by Rupprecht et al. Rupprecht et al. (2015) in terms of Legendre polynomials. Finally, Caginalp and Chen Caginalp and Chen (2012) derived a fully explicit closed formula for the mean first-passage time which writes

[TABLE]

with respect to the complex variable $z=r\exp(i\theta)$ such that $r=|{\bf x}|/{\cal R}$ and $\theta$ represents the angle between ${\bf x}$ and the center of the escape region ${\bf x_{0}}$. This exact expression is asymptotically compatible with Eq. (16) in the $\varepsilon\ll 1$ limit. We also mention here that the result of Rupprecht et al. Rupprecht et al. (2015) provides the full distribution of the first passage times.

To solve the annulus and two-shell geometry problems, we will use the same formalism to solve the equations for the regular part of the pseudo Green’s function $R({\bf x}|{\bf x_{0}})$. Then we apply Eqs. (11) and (14) to obtain the expressions of the GMFPT and the MFPT for a particle starting at the position ${\bf x}$, which are correct in the first sub-leading order of the small escape angle $\varepsilon$ expansion. In addition, we can mention that the recent exact solution of the escape problem for any simply connected planar domain with an arbitrary diffusivity Grebenkov (2016). This exact solution is compatible with the Eq. (19) for the disk geometry.

III The narrow escape problem for the annulus geometry

In this section, we look at the NE problem for a Brownian particle inside an annulus of width $\Delta$, ${\cal R}-\Delta\leq|{\bf x}|\leq{\cal R}$, denoted as $\Omega$, with an area $|\Omega|=\pi\Delta(2{\cal R}-\Delta)$. The outer boundary $|{\bf x}|={\cal R}$ is denoted as $\partial\Omega$ while the inner boundary $|{\bf x}|={\cal R}-\Delta$ is denoted as $\partial\Omega_{\Delta}$. The escape region is part of the outer boundary defined by $|\theta|\leq\varepsilon/2$ and is denoted as $\partial\Omega_{\varepsilon}$. The opposite case with an escape region on the inner boundary was studied by Holcman et al. Schuss et al. (2007). This problem is equivalent to the two-shell geometry in the limit case $D_{0}=0$ by which the diffusive particle is excluded from the inner shell ($|{\bf x}|\leq{\cal R}-\Delta$). The equations for the MFPT for a Brownian particle starting at the position ${\bf x}$, denoted as $t({\bf x})$, are deduced from the Eq. (1) and the backward equations for the probability density function Redner (2001)

[TABLE]

This problem can be solved with only two dimensionless parameters: the dimensionless annulus width $\delta=\Delta/{\cal R}$ with respect to the radius ${\cal R}$ and the angle $\varepsilon$ of the escape region, after rescaling the length by the radius ${\cal R}$ and the time by ${\cal R}^{2}/D_{\Delta}$. Numerical solutions of these PDEs are shown in the Fig. 3 for several values of the annulus width $\Delta\in[0,{\cal R}]$ and the escape angle $\varepsilon=0.2$. When $\Delta\rightarrow{\cal R}$, the solutions of the section II are valid, replacing $D_{0}$ by $D_{\Delta}$ in analytical expressions, comparing the Fig. 3f with the Fig. 2a.

Fig. 3 demonstrates several interesting features of the MFPT: by comparing (a) with (f) one sees that particles starting close to the escape region find the target faster in the thin annulus (a), i. e. small $\Delta$, than in the thick annulus (f), i. e. large $\Delta$. The inaccessible central region represents an obstacle impeding the particles starting close to the escape region from diffusing away when it is sufficiently large. On the other hand particles starting far away from the escape region, for instance on the opposite side, reach the target faster in the thick annulus (f), i. e. large $\Delta$, than in the thin annulus (a), i. e. small $\Delta$. Those particle are impeded by the central obstacle to move quickly towards the target, they remain jammed for a longer time when the obstacle is large. Thus the MFPT for particles starting close to the escape region increases with increasing $\Delta$, whereas the MFPT of particles starting far away from the escape region decreases with decreasing $\Delta$. After averaging over the initial positions these two opposing tendencies could in principle lead to a minimum of the GMFPT as function of $\Delta$. This is actually what we will find in the next subsection.

III.1 The narrow escape limit ($\varepsilon\ll 1$)

The derivation of the sub-leading order of the MFPT is presented in appendix A, where it is shown that the MFPT satisfies Eq. (14) with the volume replaced by $|\Omega|=\pi\Delta(2{\cal R}-\Delta)$ and $D_{0}$ by $D_{\Delta}$. The regular part of the Green’s function $R({\bf x}|{\bf x_{0}})$ is calculated in the appendix B.2 and its expression is given by Eq. (90) in terms of (dimensionless) polar coordinates. Its value at the center of the escape region is then equal to

[TABLE]

which implies the GMFPT expression from Eq. (11)

[TABLE]

In Fig. 4, the GMFPT is plotted as a function of the escape angle $\varepsilon$ for several annulus widths $\Delta\in[0,{\cal R}]$. In the limit $\varepsilon\ll 1$, we observe that the analytical prediction given by Eq. (25) is verified in comparison with the numerical solutions obtained with FreeFem++. Moreover, we remark in the inset of the Fig. 4 that the value of the regular part of the Green’s function at ${\bf x_{0}}$ plotted with respect to the dimensionless annulus width $\delta$, extract from numerical solutions, is compatible with the Eq. (24).

The value of $R({\bf x_{0}}|{\bf x_{0}})$ depends on the function $S(X)$ defined by the series expression (for $|X|<1$)

[TABLE]

When $X\rightarrow 0$, this function behaves like $S(X)\simeq-\ln(1-X^{2})$. Then, in the limit $\delta\rightarrow 1$, the regular part of the pseudo Green’s function behaves like $R({\bf x_{0}}|{\bf x_{0}})\sim 1/8\pi$, corresponding to the result of the section II which implies that the GMFPT satisfies the Eq. (15). When $X\rightarrow 1$, the leading order of the series is $S(X)\simeq\pi^{2}(1-X)^{-1}/12$. Thus, in the limit $\Delta\rightarrow 0$, the regular part of the pseudo Green’s function behaves like $R({\bf x_{0}}|{\bf x_{0}})\sim\pi{\cal R}/6\Delta$ which implies the leading order of the GMFPT in the $\varepsilon\rightarrow 0$ and $\delta\rightarrow 0$ limits

[TABLE]

When $\varepsilon=2\pi$, the mixed Neumann-Dirichlet boundary condition becomes a fully absorbing (Dirichlet) boundary condition. The MFPT for a particle starting at ${\bf x}$ is then radially symmetric and the GMFPT expression is

[TABLE]

The expression (25) of the GMFPT in the narrow escape limit can be decomposed as the sum of the GMFPT to reach the boundary $\partial\Omega$ given by Eq. (28) and the GMFPT starting at the boundary $\partial\Omega$ to reach the escape region which writes in the NE regime as

[TABLE]

Moreover, the maximal value of the MFPT is located in the NE regime at $-{\bf x_{0}}=(-{\cal R},0)$. Its expression is then obtained from Eqs. (14) and (90) with ${\bf x}=-{\bf x_{0}}$

[TABLE]

In the $\varepsilon\rightarrow 0$ and $\delta\rightarrow 0$ limits, this expression of the MMFPT goes to the value

[TABLE]

In the opposite regime ($\varepsilon\simeq 2\pi$), the maximal value of the MFPT is located on the circle of radius $|{\bf x}|={\cal R}-\Delta$ and its expression is then given by

[TABLE]

for a fully absorbing boundary.

The GMFPT calculated by applying the PDE solver FreeFem++ to Eqs. (20-23) is shown in Fig. 5a as a function of the width of the annulus $\Delta$ for a fixed value of $\varepsilon$. One observes that the GMFPT approaches a constant in the thin annulus limit. Furthermore, the GMFPT is minimized for a specific value of the width $\Delta$ of the annulus. For an escape angle $\varepsilon\lesssim 1$, the minimum is pronounced and clearly visible in Fig. 5a, whereas for smaller values of $\varepsilon$ the minimum is very shallow and nearly invisible but indicated by the black dots in Fig. 6.

In connection with Fig. 3 we have already realized that the MFPT for particles starting close to the escape region increases with increasing $\Delta$, whereas the MFPT of particles starting far away from the escape region decreases with decreasing $\Delta$ such that these opposing tendencies could lead to a minimum of the GMFPT as a function of $\Delta$. This picture is confirmed by the behavior of the GMFPT presented in Fig. 6a, where the black dots indicate the value of $\Delta$ where for fixed $\varepsilon$ the minimum occurs. The starting point with the maximum MFPT is always $-{\bf x_{0}}=(-{\cal R},0)$ and Fig. 6b shows clearly that for $\varepsilon<0.15$ the MFPT of a particle starting opposite to the escape region first decreases with increasing $\Delta$, reaches a minimum at a specific width and then increases with $\Delta$.

It is an interesting question whether the existence of the minimum discussed here persists when the impenetrable central obstacle, $|{\bf x}|<{\cal R}-\Delta$, is replaced by a soft repulsive barrier, which may occur because of geometric constraints in a cell. This situation is modeled by a Bessel process Ryabov et al. (2015) which is an overdamped Brownian motion in a logarithmic potential $U({\bf x})$ for the inverse temperature $\beta$

[TABLE]

repulsive when $g>0$ and presenting a divergence at the center of the disk. $\Theta$ denotes here the Heaviside function. The MFPT starting at ${\bf x}$ to reach the escape region of angle $\varepsilon$ satisfies then the equations

[TABLE]

The boundary condition on $|{\bf x}|={\cal R}-\Delta$ is given by the conservation of the backward flux: $-D_{\Delta}{\bf n}\cdot\nabla t({\bf x})$ which does not depend on the value of the potential. The gradient of potential is

[TABLE]

When $g=0$, the MFPT of the disk geometry is recovered (data not shown). In Fig. 7a, the numerical solution MFPT obtained with the PDE solver FreeFem++ is shown for the width $\Delta=0.25{\cal R}$, the escape angle $\varepsilon=0.2$ and the potential strength $g=100$. The solution in the inner shell $|{\bf x}|<{\cal R}-\Delta$ is only $\theta$-dependent since the Brownian particle moves rapidly to the outer shell following the radial potential, equivalent to the boundary condition ${\bf n}\cdot\nabla t=0$ on the inner boundary. This solution is compatible with the annulus geometry in the outer shell (see Fig. 3b). A difference appears only when the potential is too soft, for example $g<10$, since MFPT is considerably changed in the limit $\Delta\ll{\cal R}$.

In Fig. 7b, the GMFPT is shown as a function of the width $\Delta$ for $\varepsilon=0.3$ and several values of $g$ and compared to the GMFPT obtained for the annulus geometry. For $g=1$, the GMFPT of the annulus geometry is compatible only for $\Delta>0.9{\cal R}$ and the minimum is reached for $\Delta=0$ instead of $\Delta\simeq 0.4{\cal R}$. A minimum exists for a width $\Delta\neq 0$ for $g\gtrsim 10$ and for $g>50$ the $\Delta$-dependance of the GMFPT is similar to the one of the annulus geometry. The minimum of the GMFPT is then not present for all soft potentials and a sufficiently strong repulsion is needed for it to appear. The lower limit of $g$, denoted as $g_{l}$, for the existence of a minimum depends on the escape angle $\varepsilon$ (data not shown) and decreases with $\varepsilon$.. In the narrow escape regime ($\varepsilon\ll 1$), this limit value is then large $g_{l}\gg 1$, which implies that the minimum appears only for a strong potential, like in the annulus geometry.

III.2 The pseudo one-dimensional limit ($\Delta\ll{\cal R}$)

In the $\Delta\ll{\cal R}$ limit, we can consider a pseudo one-dimensional diffusion on the surface $|{\bf x}|={\cal R}$ and Eqs. (20-21) become for the function $t(\theta)$ dependent only on the orthoradial coordinate

[TABLE]

while the Eqs. (22-23) are automatically satisfied. The unique solution of this Cauchy system is

[TABLE]

The GMFPT expression is then given by the spatial average over $\theta\in[0,2\pi]$ of the last expression

[TABLE]

This expression reproduces the asymptotic value $\pi^{2}/3$ of the GMFPT as predicted by Eq. (27) in the limit $\varepsilon\rightarrow 0$ and $\delta\rightarrow 0$ simultaneous under the condition $-\delta\ln\varepsilon\ll 1$. The limit $\varepsilon\rightarrow 2\pi$ gives a vanishing GMPFT which is plausible since the whole domain is absorbing the Brownian particle. The Eq. (41) is shown in Fig. 5b and compared with the numerical solution for $\delta=10^{-3}$ and several values for the escape angle $\varepsilon\geq 10^{-2}$. We conclude that the prediction of the pseudo-one-dimensional approximation is valid in the limit $\delta\ll 1$ for a fixed value of escape angle $\varepsilon$.

The NE expression of the GMFPT given by Eq. (25) is then valid for small values of $\varepsilon$ such that $-\delta\ln\varepsilon>1$ whereas the pseudo-1d expression of the GMFPT given by Eq. (41) is valid for small values of $\delta$ such that $-\delta\ln\varepsilon\ll 1$. These two validity domains are shown in Fig. 6 and a comparison between analytical solutions and numerical solution in both domains is proposed in Fig. 5a. We can then see that the minimum value of the GMFPT is reached in the narrow escape regime.

In the $\Delta\ll{\cal R}$ limit, the maximal MFPT is reached at $\theta=\pi$ and its expression is deduced from Eq. (40)

[TABLE]

A numerical verification of this expression is also shown in Fig. 5b for $\delta=10^{-3}$ and several values of escape angle $\varepsilon\geq 10^{-2}$. In the $\varepsilon\rightarrow 0$ and $\delta\rightarrow 0$ limits, such that $-\delta\ln\varepsilon\leq 1$, this expression is equivalent to Eq. (31). The NE expression of the MMFPT given by Eq. (30) is then valid for small values of $\varepsilon$ such that $-\delta\ln\varepsilon>1$ whereas the pseudo-1d expression of the GMFPT given by Eq. (42) is valid for small values of $\delta$ such that $-\delta\ln\varepsilon\ll 1$. These validity domains are shown in Fig. 5b.

IV Narrow escape problem for the two-shell geometry

In this section, we investigate the NE problem for the two-shell geometry where the search domain is separated into two concentric regions: the outer shell of width $\Delta$ for a radius in ${\cal R}-\Delta\leq|{\bf x}|\leq{\cal R}$ denoted as $\Omega_{\Delta}$ and the inner shell for a radius $|{\bf x}|\leq{\cal R}-\Delta$. The inner boundary located at $|{\bf x}|={\cal R}-\Delta$ is denoted as $\partial\Omega_{\Delta}$ whereas the outer boundary located at $|{\bf x}|={\cal R}$ is denoted as $\partial\Omega$. The escape region is located on the outer boundary for $|\theta|\leq\varepsilon/2$ and is denoted as $\partial\Omega_{\varepsilon}$. The diffusivity of the Brownian particle is different in the two shells, denoted as $D_{0}$ and $D_{\Delta}$ in the inner shell and the outer shell respectively (see Fig.1c). The MFPT for a diffusive particle starting in the outer shell and in the inner shell is denoted as $t_{1}({\bf x})$ and $t_{2}({\bf x})$ respectively. The equations satisfied by the functions $t_{i}({\bf x})$ are deduced from the Eq. (1) and the backward equations for the probability density function Redner (2001)

[TABLE]

The boundary conditions (45) on the inner boundary guarantee the continuity of the concentration and the flux of the Brownian particles across the interface, implying the continuity of the probability density function $p({\bf y},\tau|{\bf x})$ and the normal component of the current defined by the Fick’s law as $-D({\bf x})\nabla p({\bf y},\tau|{\bf x})$. Note that the model studied here and the resulting equations (43-46) are different from the diffusion equation with a spatially varying diffusion constant $D({\bf x})\nabla^{2}t({\bf x})=-1$, which has been studied and solved exactly in Ref. Grebenkov (2016).

Numerical solutions of these PDEs are shown in Fig. 8 for $\Delta=0.25{\cal R}$ and several values of the ratio $D_{\Delta}/D_{0}$. This problem has only three dimensionless parameters: the ratio of diffusivity between outer and inner shells $D_{\Delta}/D_{0}$, the width of the outer shell $\delta=\Delta/{\cal R}$ and the angle $\varepsilon$ of the escape region, after rescaling the lengths by the radius ${\cal R}$ and the time by ${\cal R}^{2}/D_{\Delta}$.

IV.1 The narrow escape limit ($\varepsilon\ll 1$)

In the limit $\varepsilon\rightarrow 0$ the MFPT satisfies in sub-leading order Eq. (14) where the volume of the space is $|\Omega|=\pi{\cal R}^{2}$ and $D_{0}$ is replaced by the diffusivity of the region close to the hole: $D_{\Delta}$. The regular part of the Green’s function $R({\bf x}|{\bf x_{0}})$ is calculated in the appendix B.3 and its expression is given by Eqs. (114-115) in terms of (dimensionless) polar coordinates. Its value at the center of the escape region is then equal to

[TABLE]

We observe that the value $R({\bf x_{0}}|{\bf x_{0}})$ goes to $1/8\pi$ when $D_{\Delta}\rightarrow D_{0}$ or when $\Delta\rightarrow{\cal R}$. In fact, in these two limits, the physical problem is identical to the disk geometry discussed in Sec. II. In the limit $D_{\Delta}\gg D_{0}$, we note that $R({\bf x_{0}}|{\bf x_{0}})$ diverges, implying a divergence of the GMFPT. This is due to the fact that the particles starting at $|{\bf x}|<{\cal R}-\Delta$ are slowed down by the small diffusivity in the central core. Moreover, we note that in the limit $\Delta\rightarrow 0$, the value of $R({\bf x_{0}}|{\bf x_{0}})$ diverges as $\ln(\Delta/{\cal R})$. The expression of the GMFPT is then given by Eq. (11)

[TABLE]

From the Eqs. (14) and (115), the expression of the CMFPT for a diffusive particle starting at the center of the disk is in the narrow escape limit

[TABLE]

The result for the GMFPT and CMFPT, Eqs. (48) and (49), valid in the narrow escape limit $\varepsilon\rightarrow 0$, are shown in Fig. 9a as a function of $\varepsilon$ for $D_{\Delta}=5D_{0}$ and $\delta=0.25$. The analytical solutions are compared with numerical solutions of the PDEs (43-46) obtained with FreeFem++ and with the results of Brownian particle simulations. In figure 9b, the analytical predictions for $R({\bf x_{0}}|{\bf x_{0}})$ given by Eq. (47) are compared with the numerical fits of the constant sub-leading order of the GMFPT obtained with FreeFem++.

For the numerical simulation of Brownian motion, we apply the algorithm presented in the appendix C to the radial coordinate of the Brownian particle. When the particle is on the inner boundary defined by the identity $|{\bf x}|={\cal R}-\Delta$, the probability to go to the inner or the outer shell is given respectively by Eq. (126) and Eq. (127). The movement is then decomposed into radial and orthoradial coordinates assuming that - for very small displacement - the inner boundary is approximatively a wall with a zero curvature. This algorithm is valid for any values of $\delta$ and $\varepsilon$. The relative error of these stochastic simulations is inversely proportional to the square-root of their number of realizations $N=10^{6}$. The relative error is then of order $10^{-3}$.

Furthermore, when $\varepsilon=2\pi$, the mixed Neumann-Dirichlet boundary condition becomes a fully absorbing (Dirichlet) boundary condition implying that the MFPT function is radial. The GMFPT and CMFPT expressions are respectively

[TABLE]

The GMFPT and the CMFPT expressions in the narrow escape limit given respectively by Eqs. (48) and (49) can be decomposed into two terms: respectively the GMFPT and the CMFPT to reach the surface given by Eqs. (50) and (51) and the GMFPT to reach the escape region starting on the surface, which satisfies the expression

[TABLE]

Moreover, the starting point that has the maximum MFPT is not always located at $-{\bf x_{0}}=(-{\cal R},0)$ [see Fig. 8]. In fact, if the diffusivity of the inner shell is much bigger than the one of the outer shell ($D_{\Delta}>D_{0}$), the time spent in the inner shell can be bigger than the one spent in the outer shell which implies that the MFPT to reach the escape region - in the outer shell - is maximized in the inner shell. Because of the symmetries of the problem (concentric and symmetric with respect to the horizontal axis), the maximum must be located on the horizontal axis, since it depends on the distance to the escape region, and more precisely on the segment between $-{\bf x_{0}}$ and ${\bf 0}$. We denote the distance between the starting point and the center of the disk that gives the maximum MFPT (MMFPT) with $r_{\rm max}$ and with $D_{\Delta}^{c}>D_{0}$ the lower limit of the diffusivity $D_{\Delta}$ for which the starting point with the MMFPT is located in the inner shell. If $D_{\Delta}<D_{\Delta}^{c}$, the MMFPT starting point is located in the outer shell ($|{\bf x}|>{\cal R}-\Delta$) by definition of $D_{\Delta}^{c}$. The value of $r_{\rm max}$ is then restricted between ${\cal R}$ and ${\cal R}-\Delta$. Since the MFPT is an increasing function of the distance to the escape region, when the diffusivity is constant, the MMFPT is then necessarily situated at $-{\bf x_{0}}$ and $r_{\rm max}={\cal R}$ [see Fig. 8c]. If $D_{\Delta}>D_{\Delta}^{c}$, the MMFPT starting point is located in the inner shell ($|{\bf x}|<{\cal R}-\Delta$) by definition of $D_{\Delta}^{c}$. The value of $r_{\rm max}$ is then lesser than ${\cal R}-\Delta$. For a diffusivity $D_{\Delta}\gtrsim D_{\Delta}^{c}$, the MMFPT starting point is located in the maximum distance to the escape region within the inner shell implying that $r_{\rm max}\lesssim{\cal R}-\Delta$ [see Fig. 8d].

When $D_{\Delta}\gg D_{\Delta}^{c}$, the MFPT can be decomposed into the sum of the MFPT to reach the inner boundary ${\bf x}={\cal R}-\Delta$ and the MFPT to reach the escape region from this last point. The MMFPT starting point is then located at the center of the disk which implies that $r_{\rm max}=0$ [see Fig. 8f]. Since the MMFPT starting point stays for $D_{\Delta}>D_{\Delta}^{c}$ in a region of constant diffusivity, $r_{\rm max}$ decreases continuously from ${\cal R}-\Delta$ to [math].

To summarize, the location of the MMFPT moves discontinuously from the position $-{\bf x_{0}}$ ($r_{\rm max}={\cal R}$) for $D_{\Delta}<D_{\Delta}^{c}$ to a position inside the inner shell ($r_{\rm max}<{\cal R}-\Delta$) for $D_{\Delta}>D_{\Delta}^{c}$ and approaches continuously the center of the disk ($r_{\rm max}=0$) in the limit $D_{\Delta}\gg D_{\Delta}^{c}$. This behaviour is demonstrated in Fig. 10a for $\Delta=0.25{\cal R}$ and $\varepsilon=0.05$. The limit value of the diffusivity at the discontinuity is equal to $D_{\Delta}^{c}\simeq 4.5D_{0}$. In Fig. 10b, we show the position of the maximum $r_{\rm max}$ with respect to the parameters $\Delta/{\cal R}$ and $D_{\Delta}$ for the escape angle $\varepsilon=0.05$. The function $D_{\Delta}^{c}(\Delta)$ can be extracted from the discontinuity of the location of the MMFPT and appears to be continuous. The discontinuity is less pronounced for $\Delta\ll{\cal R}$ since the width of the outer shell is small while the discontinuity is more pronounced for $\Delta\simeq{\cal R}$ but occurring for a larger diffusivity - for a geometry close to the disk.

IV.2 The thin outer shell limit ($\Delta\ll{\cal R}$)

Next we calculated the GMFPT numerically with FreeFEM++ and show the result in Fig. 11a as a function of the width $\Delta$ for fixed values of the escape angle $\varepsilon$ and for $D_{\Delta}=5D_{0}$. We compare it with the narrow escape expression (48) and note that it is valid only for not too small values of $\Delta$ since the analytical solution diverges when $\Delta\ll{\cal R}$ while the numerical solution approaches a constant. This behavior shows that the diffusive motion is enhanced by the excursions in the thin outer shell ($\Delta\ll{\cal R}$), which implies a finite GMPFT to reach a small escape region on the outer boundary. This kind of behavior was already observed for surface-mediated diffusion problems Bénichou et al. (2010) where the leading order of the narrow escape problem is constant and no more logarithmic divergent with the escape angle $\varepsilon$. This behavior only happens for $D_{\Delta}>D_{0}$ as shown in Fig. 11b. In fact, for $D_{\Delta}=D_{0}$ the two-shell geometry reduces to the disk geometry studied in section II presenting the classical logarithmic behavior. And for $D_{\Delta}<D_{0}$, the behavior is super-logarithmic since the GMFPT diverges much faster than the logarithmic behavior of the $D_{\Delta}=D_{0}$ case: the Brownian motion in the outer shell is slowed down due to the diffusion constant value, smaller than in the inner shell, and the MFPT to reach the escape region located in it, is strongly increased.

Moreover, Fig. 12 shows that the GMFPT is a monotonic function of $\Delta/{\cal R}$ and $D_{\Delta}/D_{0}$ for a fixed value of $\varepsilon$. The same observation holds for the CMFPT (the data is not shown here). The GMFPT as a function of $\Delta$ is a strictly decreasing function if and only if $D_{\Delta}>D_{0}$. When the value of $\Delta$ is increased, the GMFPT is decreased if $D_{\Delta}>D_{0}$ since the area with the biggest diffusivity (the outer shell) increases. Strictly the opposite behavior happens when $D_{\Delta}<D_{0}$ since the area with the lowest diffusivity increases. The GMFPT as a function of $D_{\Delta}/D_{0}$ is a decreasing function since the MFPT is inversely proportional to the diffusivity. Then, there is no minimum of the MFPT for the two-shell geometry, i. e. for $\Delta\neq 0$ and $\Delta\neq{\cal R}$.

The limiting case $D_{0}\ll D_{\Delta}$ is not equivalent to the annulus geometry studied in the section III since the probability to find a particle in the inner shell is not zero. In fact, the diffusive particle can enter in the inner shell and the stationary probability (without the hole) to find it in the inner shell is the same as in the outer shell and equal to $1/|\Omega|$, at the contrary of the annulus geometry for which the Brownian particle is excluded from the inner shell. The MFPT starting far from the escape region in the outer shell will then be impacted by this difference: the particle will have a high probability, equal to $(1-\delta)^{2}$ from ergodicity, to enter in the inner shell which plays the role of a trap where the particle stays a long time since $D_{0}\ll D_{\Delta}$. Since the inner boundary does not play the role of a barrier but just separates two regions with different spreading of particles, the MFPT obtain for the annulus geometry will be considerably modified and then explain the differences between both solutions.

V Discussion

We have studied the narrow escape problem for two different disk-like geometries with an inner and an outer region, the latter of width $\Delta$. The annulus geometry for which the Brownian particle is excluded from the center shell and the two-shell geometry where the particle has two different diffusion constants. The narrow escape problem consists to look at the asymptotic behavior of the MFPT when the escape region is reduced to a small escape angle $\varepsilon\ll 1$ when the width $\Delta$ is fixed. Furthermore, we analyzed the behavior of the MFPT in the thin outer shell limit $\Delta\ll{\cal R}$ when the size of the escape region $\varepsilon$ is fixed.

For the annulus geometry, both limits are compatible such that the expression of the global mean first passage time (GMFPT) is given at the leading order by

[TABLE]

for the dimensionless width $\delta=\Delta/{\cal R}$. Whereas the GMFPT is an increasing function of $\varepsilon$, it can be optimized with respect to the width $\Delta$ of the annulus. The minimum of GMFPT can be explained by the reduction of the exploration area for Brownian particles starting close to the escape region while the particles starting far from it remain jammed for a longer time when the obstacle is large.

For the two-shell geometry, both limits are not interchangeable since the GMFPT is equal to $\langle t\rangle\propto-\ln(\varepsilon/4)+C(\Delta)$ in the narrow escape limit with a sub-leading order term that diverges in the thin outer shell limit: $C(\Delta)\underset{\Delta\rightarrow 0}{\propto}\ln(\Delta/{\cal R})$, whereas the limit of the GMFPT is constant only when $D_{\Delta}>D_{0}$. A similar behavior was recently derived Grebenkov et al. (2017) for the MFPT to reach a small target of height $\varepsilon$ on a cylinder of radius $\Delta$, exhibiting similar asymptotically logarithmic features. Moreover, we found that the distance of the starting point with the maximum MFPT jumps discontinuously from ${\cal R}$, corresponding to a point opposite to the escape region, to a value close to ${\cal R}-\Delta$, corresponding to a point within the inner shell, when the outer diffusion constant $D_{\Delta}$ is increased from $D_{0}$ to a sufficiently high value.

It turns out that the GMFPT is a strictly monotonous function for all three main parameters : $\varepsilon$, $\Delta/{\cal R}$ and $D_{\Delta}/D_{0}$. This is in contrast to the optimizability of the MFPT reported in Refs. Schwarz et al. (2016a, b). In a model for spatially inhomogeneous intermittent search that alternates stochastically between ballistic and diffusive motion where ballistic motion is only radial in the inner shell and multi-directional in the outer shell, it was found that the MFPT for the narrow escape problem can be minimized for a thin outer shell, provided the outward radial transport in the inner shell is enhanced. Thus in the model considered in Refs. Schwarz et al. (2016a, b), it is the enforced re-injection of the searcher into the outer shell that is crucial for the optimization effect. It leads to an increased probability to find the particle in the outer shell, which is where the target is: the narrow escape region. On the other hand, the optimizability of the MFPT found for surface mediated diffusion Bénichou et al. (2010) is due to a low desorption and a high absorption probabilities to the surface also enforcing the particle to stay on the surface, where again the target is. Analogously, we mention the result of Ref. Grebenkov and Oshanin (2017) concerning the optimization of the MFPT with respect to the range of attractive and repulsive potential on the surface which creates a competition between bulk and boundary events.

This analysis shows that it would be interesting to introduce for the model considered in this paper a similar mechanism leading to an increased probability to find the target in the outer shell. A straightforward possibility would be to make the boundary between the inner and outer shell semi-permeable, allowing only transitions from the inner to the outer shell, to favor the finding of the target in the outer shell. A less constrained possibility would be to introduce appropriate asymmetric reflection and transmission probabilities (or energy barriers) for this boundary. Alternatively, one could consider stochastic resetting into the cortex region when leaving the outer shell Majumdar et al. (2015). Furthermore, these analytical solutions of the narrow escape problem for a two-shell geometry let suppose that the spatially inhomogeneous intermittent search Schwarz et al. (2016a, b) in this two-shell geometry could be solvable in the narrow escape limit using the same method. These extensions of the problem analyzed in the present paper will be considered in a forthcoming publication.

VI Acknowledgment

This work was performed with financial support from the Deutsche Forschungsgemeinschaft DFG within the Collaborative Research Center SFB 1027. M. M. thanks T. Guérin for interesting discussions.

Appendix A Narrow escape expression of the MFPT

In this appendix, we present the derivation of the sub-leading term of the MFPT mainly inspired by the method of Ward and Keller Ward and Keller (1993) which was improved by Pillay et al. Pillay et al. (2010) and Chevalier et al. Chevalier et al. (2011). The escape region with an angle $\varepsilon\ll 1$ is considered as a perturbation on the boundary. The solution far from the hole (outer solution) takes on the form

[TABLE]

where $\tau({\bf x})$ represents the space-dependent sub-leading term. The multiplicative constant $4$ is just taken for convenience. It is of order ${\cal O}(\varepsilon^{0})$ usually included in the unknown $\tau$ function. From Eqs. (2-4), the function $\tau({\bf x})$ satisfies the equations

[TABLE]

The solution close to the hole (inner solution) is written in terms of the inner variable ${\bf\tilde{x}}=2({\bf x}-{\bf x_{0}})/\varepsilon{\cal R}$, such that the hole corresponds to the segment (independent of the value of $\varepsilon$) defined by $\tilde{x}=0$ and $|\tilde{y}|\leq 1$, and the reflecting boundaries corresponds to the half-lines defined by $\tilde{x}=0$ and $|\tilde{y}|\geq 1$, considering $\varepsilon$ enough small to have a flat boundary. In these coordinates, one writes $t({\bf x})\equiv v({\bf\tilde{x}})$. The equations for $v({\bf\tilde{x}})$ are then for $\varepsilon\rightarrow 0$

[TABLE]

These equations can be solved in the elliptic coordinates ($\mu$, $\nu$) defined by $\tilde{x}=\sinh\mu\sin\nu$ and $\tilde{y}=\cosh\mu\cos\nu$. The Laplace’s equation writes in this system of coordinates: $\partial_{\mu\mu}v+\partial_{\nu\nu}v=0$. The escape region corresponds to $\mu=0$ and $\nu\in[0,2\pi]$, implying the boundary condition $v(\mu=0,\nu)=0$. The reflecting boundary corresponds to $\nu=0$ and $\mu>0$, involving the condition $\partial_{\nu}v(\mu,\nu=0)=0$. The $\nu$-independent solution of this system is $v(\mu,\nu)=A\mu$. The outer and the inner solutions are matched in an intermediate region such that ${\bf x}\rightarrow{\bf x_{0}}$ and $|{\bf\tilde{x}}|\rightarrow\infty$ (i. e. $\mu\rightarrow\infty$ for elliptic coordinates). The inner solution behaves thus like

[TABLE]

Comparing with the outer solution (54), this implies the value of $A=|\Omega|/(\pi D_{0})$ and the origin of the factor $4$ is seen. The behavior of the function $\tau({\bf x})$ close to the hole is thus given by

[TABLE]

Following Ref. Pillay et al. (2010), the pseudo Green’s function $G({\bf x}|{\bf x_{0}})$ is introduced via the equations

[TABLE]

where $\delta({\bf x})$ represents the Dirac’s distribution which can be omitted here since ${\bf x_{0}}$ belongs to the boundary of the domain and $R({\bf x_{0}}|{\bf x_{0}})$ is the (unknown) regular part of the Green’s function at the center of the escape region. The function $\tau({\bf x})$ is then given by the relation obtained from Eqs. (55-56) and (61)

[TABLE]

where the spatial average $\langle\tau\rangle$ satisfies

[TABLE]

The spatial average of the MFPT is then given by

[TABLE]

reproduced inside the main text via Eq. (11) and the expression of the MFPT $t({\bf x})$ is

[TABLE]

corresponding to the given Eq. (12) in the main text. These both equations depends on the pseudo Green’s function defined by Eqs. (62-65) reproduced in the main text via Eqs. (7-10) where the Dirac distribution is removed.

Appendix B Computations of the regular part of the Green’s function $R({\bf x}|{\bf x_{0}})$

In this appendix, we will compute the expression of the regular part of the Green’s function $R({\bf x}|{\bf x_{0}})$, defined by Eq. (13), for the three different geometries shown in Fig. 1. We will start to a rewriting of Eqs. (7-10) for this regular part. The Eq. (9) is satisfied by the definition of $R({\bf x}|{\bf x_{0}})$. Since the function $\ln|{\bf x}-{\bf x_{0}}|$ is a solution of the Laplace’s equation, for all planar domains where ${\bf x}\neq{\bf x_{0}}$, the equation (7) becomes

[TABLE]

Moreover, the spatial average of $\ln(|{\bf x}-{\bf x_{0}}|/{\cal R})$ vanishes over the disk - via the integral over the orthoradial part - the equation (10) gives

[TABLE]

The condition on the boundary of the domain (8) is changed to

[TABLE]

We use the (dimensionless) polar coordinates $(r,\theta)$ defined as $r=|{\bf x}|/{\cal R}$ and ${\bf x}\cdot{\bf x_{0}}/{\cal R}^{2}=r\cos\theta$. This definition allows us to rename the regular part of the Green’s function as $R(r,\theta)$ without any loss of generality. Since the normal vector ${\bf n}$ is always radial for all studied geometries, the relation for the radial derivative of $\ln|{\bf x}-{\bf x_{0}}|$ at the point $(r,\theta)$ is then

[TABLE]

B.1 Known solution: the disk geometry

The equation on the boundary of the domain (72) rewrites

[TABLE]

which is $\theta$-independent. We can thus look for a radial solution, given by

[TABLE]

involving the value at the center of the escape region $R({\bf x_{0}}|{\bf x_{0}})\equiv R(1,0)=1/8\pi$.

B.2 The annulus geometry

The equations of the Green’s function for the annulus geometry can be deduced from the Eqs. (20-23) of the MFPT. Its regular part satisfies then the bulk equation (70) and the normalization condition (71). From Eq. (72), the boundary conditions respectively on the outer and the inner boundaries become

[TABLE]

where we used the relation (73). Since this last equation is $\theta$-dependent, the solution $R(r,\theta)$ is not radially symmetric for this geometry. We thus look for a solution with the global form

[TABLE]

due to the periodicity of the polar coordinate $\theta$ and the symmetry of the problem under the $\theta\rightarrow-\theta$ transformation. The functions $f_{n}(r)$ should satisfy the equations resulting from the bulk equation (70)

[TABLE]

The global solution is thus given by $f_{n}(r)=a_{n}r^{n}+b_{n}r^{-n}$ (with $n\geq 1$) and $f_{0}(r)=r^{2}/(4|\Omega|)+a_{0}+b_{0}\ln r$, which just rewrites as

[TABLE]

We will now look at the expression for $a_{n}$ and $b_{n}$ coefficients. The first condition we use is given by the equation (71). Since the average value of $\cos(n\theta)$ vanishes for $n\geq 1$, the remaining equation is

[TABLE]

From the equation (76), we get the relation

[TABLE]

which implies immediately from the orthogonality of $\cos(n\theta)$ that $a_{n}=b_{n}$ for $n\geq 1$ and

[TABLE]

This last relation implies the value of $a_{0}$ from the expression (82)

[TABLE]

The boundary condition at $r=1-\delta$ given by Eq. (77) gives the relation, using the equality $a_{n}=b_{n}$,

[TABLE]

We again use the orthogonality of $\cos(n\theta)$. For $n=0$, we get

[TABLE]

This implies again that the Eq. (84) is satisfied which ensures that a solution of the system exists. For $n\neq 0$, we get the relation

[TABLE]

which implies the expression of $a_{n}$ as

[TABLE]

After determining all these constants, we get the general expression of the regular part of the Green’s function $R(r,\theta)$

[TABLE]

Its value at the center of the escape region ${\bf x_{0}}$ - $R({\bf x_{0}}|{\bf x_{0}})\equiv R(1,0)$ - is then

[TABLE]

which is reproduced in Eq. (24) of the main text.

In the limit $\delta\rightarrow 1$, the expression (90) tends to the disk formula (75)

[TABLE]

for a radius which is let constant such that $r>1-\delta$. Moreover, the derivative of the function on the inner boundary (77) transforms as

[TABLE]

and is not applicable at the origin $r=0$, where $\theta$ is undefined. In the limit $\delta\rightarrow 1$, the regular part of the Green’s function of the annulus geometry tends to the one of the disk geometry except at the origin point. This will not change the MFPT expression since it is a continuous function all over the domain.

B.3 The two-shell geometry

The equations of the Green’s function for the two-shell geometry can be deduced from the Eqs. (43-46) of the MFPT. We define $R_{1}(r,\theta)$ and $R_{2}(r,\theta)$ the value of its regular part respectively on the inner shell and the outer shell in the (dimensionless) polar coordinates. The bulk equations verified by $R_{i}(r,\theta)$ are thus

[TABLE]

The boundary conditions at $r=1-\delta$ are given by

[TABLE]

where we used the relation (73), and the reflective boundary condition at $r=1$ becomes

[TABLE]

Finally, the normalization condition gives the last relation, fixing the additive constant of $R$

[TABLE]

The general solutions of the bulk equations (94) and (95) can be written as

[TABLE]

From the non singularity of the problem at $r=0$, we should impose $d_{n}=0$ for $n\geq 0$ to have a non divergent solution $R_{2}(r,\theta)$. The boundary condition at $r=1$ given by Eq. (98) writes

[TABLE]

The orthogonality of $\cos(n\theta)$ imposes thus $b_{0}=0$ and $a_{n}=b_{n}$ for $n\geq 1$. The continuity at $r=1-\delta$, expressed by the Eq. (96) becomes

[TABLE]

The orthogonality of $\cos(n\theta)$ imposes thus

[TABLE]

From the relation $\langle\cos(n\theta)\rangle=0$, the condition (98) becomes

[TABLE]

which simplifies as

[TABLE]

This equation gives thus

[TABLE]

Finally, to determine the $a_{n}$ and $c_{n}$ coefficients for $n\geq 1$, we solve the condition (97) at $r=1-\delta$

[TABLE]

The orthogonality of $\cos(n\theta)$ imposes thus

[TABLE]

which implies the expressions of $c_{n}$ and $a_{n}$ coefficients

[TABLE]

We have then determine the general expression of the regular part of the Green’s function $R(r,\theta)$ with

[TABLE]

Its value at the center of the escape region ${\bf x_{0}}$ - $R({\bf x_{0}}|{\bf x_{0}})\equiv R_{1}(1,0)$ - is then

[TABLE]

which is reproduced in Eq. (47) of the main text.

In the limit $D_{\Delta}\rightarrow D_{0}$, we can remark that the expressions (114-115) become

[TABLE]

corresponding to the disk formula (75). Moreover, in the limit $\delta\rightarrow 1$, the function $R_{1}(r,\theta)$ goes also to the disk formula and the function $R_{2}(r,\theta)$, valid only at $r=0$, is equal to $-1/8\pi$ corresponding to the value of the disk formula at $r=0$. Then, in both limit $D_{\Delta}\rightarrow D_{0}$ and $\Delta\rightarrow{\cal R}$, the regular part of the Green’s function tends to the disk formula (75), which implies that the MFPT starting at any position ${\bf x}$ of the domain correctly satisfies the well-known formula (16).

Appendix C Simulation of the Langevin equation

In this appendix we consider the one-dimensional problem on the unit segment $x\in[0,1]$, corresponding to the two-shell geometry with a fully absorbing external boundary ($\varepsilon=2\pi$), where $x$ plays the role of the radial component $|{\bf x}|/{\cal R}$. The extremity at $x=0$ is reflecting whereas the one at $x=1$ is absorbing. The diffusivity of the particle is denoted as $D_{0}$ for $x<1-\delta$ and $D_{\Delta}$ for $x>1-\delta$. The one dimensional MFPT for a Brownian particle starting at the position $x<1-\delta$ (respectively at $x>1-\delta$) is denoted as $t_{1}(x)$ (respectively $t_{2}(x)$). The equations in the bulk are then Redner (2001)

[TABLE]

The boundary conditions are $t_{1}^{\prime}(0)=0$ and $t_{2}(1)=0$, and the matching conditions at $x=1-\delta$ are $t_{1}(1-\delta)=t_{2}(1-\delta)$ and $D_{0}t_{1}^{\prime}(1-\delta)=D_{\Delta}t_{2}^{\prime}(1-\delta)$ corresponding respectively to the continuity of the concentration and the flux of diffusive particles across the boundary. The general solution is

[TABLE]

The CMFPT expression is given here by

[TABLE]

and the GMFPT is then

[TABLE]

We want now to compare these exact one-dimensional solutions to the numerical simulations’ ones. We recall here that the general one dimensional Fokker-Planck equation for a diffusive particle with a inhomogeneous diffusivity $D(x)$ writes

[TABLE]

with a current $J(x,t)$ satisfying the Fick’s semi-empirical law. The Langevin equation written with the Itō convention, corresponding to a full explicit discretization, is

[TABLE]

where $dB_{t}$ represents the Wiener process at time $t$ with a variance $dt$. The convective term of the Itō formalism is then given by $D^{\prime}(x)$. For our two-shell geometry, with piecewise constant diffusivities, this convective term is thus proportional to the Dirac’s distribution.

To get a realizable numerical simulation for this problem, we follow the algorithm shown by Lejay and Pichot (2012); Lejay and Maire (2013), satisfying the flux continuity on the inner boundary. When the diffusive particle is exactly at $x=1-\delta$, the probability to go to the left (through a region of diffusivity $D_{0}$) is

[TABLE]

whereas the probability to go to the right (through a region of diffusivity $D_{\Delta}$) is

[TABLE]

If $D_{0}=D_{\Delta}$, we get back to a random walker with $p_{\rm left}=p_{\rm right}=1/2$ equivalent to a diffusive motion with a constant diffusivity. In figure 13, we show that this algorithm is compatible with exact expressions of CMFPT and GMFPT given respectively by Eqs. (122) and (123) which justify the use of this algorithm for the 2d problem and the numerical data shown in Fig. 9.

Remarkably, we can mention that the fully implicit discretization (isothermal convention) writes Vaccario et al. (2015)

[TABLE]

instead of Eq. (125), with a zero convective term but difficult to integrate numerically.