Asymptotic Analysis of the Narrow Escape Problem in Dendritic Spine Shaped Domain: Three Dimension

TL;DR

This paper derives high-order asymptotic formulas for the mean first passage time of Brownian particles in a 3D dendritic spine-shaped domain, simplifying the complex boundary problem for better analytical understanding.

Contribution

It extends previous 2D analyses to 3D, providing a rigorous asymptotic expansion for the narrow escape problem in dendritic spines.

Findings

Derived high-order asymptotic expansion for mean first passage time

Simplified the boundary value problem from Dirichlet-Neumann to Robin-Neumann

Generalized previous 2D results to three-dimensional domains

Abstract

This paper deals with the three-dimensional narrow escape problem in dendritic spine shaped domain, which is composed of a relatively big head and a thin neck. The narrow escape problem is to compute the mean first passage time of Brownian particles traveling from inside the head to the end of the neck. The original model is to solve a mixed Dirichlet-Neumann boundary value problem for the Poisson equation in the composite domain, and is computationally challenging. In this paper we seek to transfer the original problem to a mixed Robin-Neumann boundary value problem by dropping the thin neck part, and rigorously derive the asymptotic expansion of the mean first passage time with high order terms. This study is a nontrivial generalization of the work in \cite{Li}, where a two-dimensional analogue domain is considered.

Tables

Table 1. Table 1 : Solutions and relative error for different values of ε 𝜀 \varepsilon .

$\epsilon$	$u$	$u_r$	$u_{\epsilon }$	$\left(u_{\epsilon }-u\right)/u$
0.10	145.01	145.37	144.48	-0.0036
0.09	177.53	177.45	177.01	-0.0029
0.08	222.81	223.28	222.31	-0.0022
0.07	288.57	289.11	288.11	-0.0016
0.06	389.47	390.12	389.07	-0.0010
0.05	556.12	556.91	555.80	-0.0006
0.04	861.64	862.63	861.46	-0.0002
0.03	1518.90	1520.30	1519.04	0.0001
0.02	3389.10	3391.10	3389.75	0.0002
0.01	13443.00	13456.00	13446.34	0.0002

Table 2. Table 2 : Solutions and relative error for different values of L 𝐿 L .

$L$	$u$	$u_r$	$u_{\epsilon }$	$\left(u_{\epsilon }-u\right)/u$
1.0	555.98	556.91	555.80	-0.0003
2.0	1090.80	1091.80	1090.63	-0.0002
3.0	1626.60	1627.60	1626.47	-0.0001
4.0	2163.40	2164.50	2163.30	-0.0000
5.0	2700.90	2702.30	2701.13	0.0001
6.0	3239.90	3241.20	3239.97	0.0000
7.0	3779.70	3781.10	3779.80	0.0000
8.0	4320.20	4321.90	4320.63	0.0001
9.0	4865.10	4863.80	4862.47	-0.0005
10.0	5408.30	5406.60	5405.30	-0.0006

Equations

$$\{\begin{array}{ll}\Delta u=-1,& \text{mbox}in\Omega ,\\ \frac{\partial u}{\partial \nu }=0,& \text{mbox}on\partial {\Omega }_r,\\ u=0,& \text{mbox}on\partial {\Omega }_a,\end{array}$$

$$u\approx \frac{|{\Omega }_h|L}{\pi {\epsilon }^2},$$

$$\{\begin{array}{ll}\Delta u_{\epsilon }=-1,& \text{mbox}in{\Omega }_h,\\ \frac{\partial u_{\epsilon }}{\partial \nu }=0,& \text{mbox}on\partial {\Omega }_r,\\ \frac{\partial u_{\epsilon }}{\partial \nu }+\alpha u_{\epsilon }=\beta ,& \text{mbox}on{\Gamma }_{\epsilon }:=\partial {\Omega }_h\setminus \overline{\partial {\Omega }_r},\end{array}$$

$$u_{\epsilon }(x)\approx \frac{|{\Omega }_h|}{\pi \alpha {\epsilon }^2}+\frac{|{\Omega }_h|M}{{\pi }^2\epsilon }+\frac{\beta }{\alpha }-\frac{|{\Omega }_h|}{2\pi |x-x^{\ast }|}$$

$$\{\begin{array}{ll}{\Delta }_{x}N(x,z)=-{\delta }_z,& x\in \Omega ,\\ \frac{\partial N}{\partial {\nu }_x}=-\frac{1}{|\partial \Omega |},& x\in \partial \Omega ,\\ {\int }_{\partial \Omega }N(x,z)d\sigma (x)=0,& \end{array}$$

$$N(x,z)=\frac{1}{4\pi |x-z|}+R_{\Omega }(x,z),x\in \Omega ,$$

$$\begin{cases}-\Delta_{x}R_{\partial\Omega}(x,z)=0,&x\in\Omega,\\[4.30554pt] \displaystyle\frac{\partial R_{\Omega}}{\partial\nu_{x}}\Big{|}_{x\in\partial\Omega}=-\frac{1}{|\partial\Omega|}+\frac{1}{4\pi}\frac{\langle x-z,\nu_{x}\rangle}{|x-z|^{3}},&x\in\partial\Omega.\end{cases}$$

$$N_{\partial \Omega }(x,z)=\frac{1}{2\pi |x-z|}+R_{\partial \Omega }(x,z),x\in \Omega ,z\in \partial \Omega ,$$

$$\begin{cases}\Delta_{x}R_{\partial\Omega}(x,z)=0,&x\in\Omega,\\[4.30554pt] \displaystyle\frac{\partial R_{\partial\Omega}}{\partial\nu_{x}}\Big{|}_{x\in\partial\Omega}=-\frac{1}{|\partial\Omega|}+\frac{1}{2\pi}\frac{\langle x-z,\nu_{x}\rangle}{|x-z|^{3}},&x\in\partial\Omega,~{}z\in\partial\Omega.\end{cases}$$

$$R_{\partial \Omega }(x,z)=-\frac{1}{4\pi }H(z)\ln |x-z|+v_{\partial \Omega }(x,z),$$

$$u_{\epsilon }(x)=\frac{|{\Omega }_h|}{\pi \alpha {\epsilon }^2}+\frac{|{\Omega }_h|M}{{\pi }^2\epsilon }+\frac{\beta }{\alpha }-\frac{|{\Omega }_h|}{2\pi |x-x^{\ast }|}+\Phi (x)+O(\alpha ),$$

$$M={\int }_{{\Gamma }_1}{\int }_{{\Gamma }_1}\frac{1}{2\pi |x-z|}dxdz,$$

$${\int }_{{\Gamma }_{\epsilon }}\frac{\partial u_{\epsilon }}{\partial \nu }d\sigma =-|{\Omega }_h|.$$

$$g(x)={\int }_{{\Omega }_h}N(x,z)dz,x\in {\Omega }_h,$$

$$\{\begin{array}{ll}\Delta g=-1,& \text{mbox}in{\Omega }_h,\\ \frac{\partial g}{\partial \nu }=-\frac{|{\Omega }_h|}{|\partial {\Omega }_h|},& \text{mbox}on\partial {\Omega }_h,\\ {\int }_{\partial {\Omega }_h}gd\sigma =0.& \end{array}$$

$$u_{\epsilon }(x)=g(x)+{\int }_{{\Gamma }_{\epsilon }}{N}_{\partial {\Omega }_h}(x,z)\frac{\partial u_{\epsilon }(z)}{\partial {\nu }_z}d\sigma (z)+C_{\epsilon },$$

$$C_{\epsilon }=\frac{1}{|\partial {\Omega }_h|}{\int }_{\partial {\Omega }_h}{u}_{\epsilon }(z)d\sigma (z).$$

$$\frac{\beta }{\alpha }-\frac{1}{\alpha }{\varphi }_{\epsilon }(x)=g(x)+\frac{1}{2\pi }{\int }_{{\Gamma }_{\epsilon }}\frac{1}{|x-z|}{\varphi }_{\epsilon }(z)d\sigma (z)+{\int }_{{\Gamma }_{\epsilon }}{v}_{\partial {\Omega }_h}(x,z){\varphi }_{\epsilon }(z)d\sigma (z)+C_{\epsilon },x\in {\Gamma }_{\epsilon },$$

$$\begin{array}{r}\frac{\beta }{\alpha }-\frac{1}{\alpha \epsilon }{\tilde{\varphi }}_{\epsilon }(x)=g(\epsilon x)+\frac{1}{2\pi }{\int }_{{\Gamma }_1}\frac{1}{|x-z|}{\tilde{\varphi }}_{\epsilon }(z)d\sigma (z)+\epsilon {\int }_{{\Gamma }_1}{v}_{\partial {\Omega }_h}(\epsilon x,\epsilon z){\tilde{\varphi }}_{\epsilon }(z)d\sigma (z)+C_{\epsilon },\end{array}$$

$$L[\phi ](x)$$

$$\frac{\beta }{\alpha }-\frac{1}{\alpha \epsilon }{\tilde{\varphi }}_{\epsilon }(x)=g(\epsilon x)+(L+\epsilon L_1){\tilde{\varphi }}_{\epsilon }(x)+C_{\epsilon }.$$

$$\big{[}I+\alpha\varepsilon(L+\varepsilon L_{1})\big{]}\tilde{\phi}_{\varepsilon}(x)=\varepsilon\alpha\left(\frac{\beta}{\alpha}-g(\varepsilon x)-C_{\varepsilon}\right).$$

$${\left(I+\alpha \epsilon (L+\epsilon L_1)\right)}^{-1}=I-\alpha \epsilon (L+\epsilon L_1)+O({\alpha }^2{\epsilon }^2).$$

$${\tilde{\varphi }}_{\epsilon }(x)=\epsilon \alpha \left[I-\alpha \epsilon (L+\epsilon L_1)+O({\alpha }^2{\epsilon }^2)\right]({\tilde{C}}_{\epsilon }+O(\epsilon )).$$

$${\tilde{C}}_{\epsilon }:=\frac{\beta }{\alpha }-g(x^{\ast })-C_{\epsilon }.$$

$${\tilde{\varphi }}_{\epsilon }(x)=\epsilon \alpha {\tilde{C}}_{\epsilon }-(\epsilon \alpha {)}^2{\tilde{C}}_{\epsilon }(L[1]+\epsilon L_1[1])+O(\alpha {\epsilon }^2).$$

$$\pi \alpha {\epsilon }^2{\tilde{C}}_{\epsilon }-{\alpha }^2{\epsilon }^3{\tilde{C}}_{\epsilon }{\int }_{{\Gamma }_1}(L[1]+\epsilon L_1[1])(x)=-|{\Omega }_h|+O(\alpha {\epsilon }^3),$$

$$\begin{array}{rl}{\tilde{C}}_{\epsilon }& =\left(I+\frac{\alpha \epsilon }{\pi }{\int }_{{\Gamma }_1}(L[1]+\epsilon L_1[1])(x)+O({\alpha }^2{\epsilon }^2)\right)\left(-\frac{|{\Omega }_h|}{\pi \alpha {\epsilon }^2}+O(\epsilon )\right)\\ & =-\frac{|{\Omega }_h|}{\pi \alpha {\epsilon }^2}-\frac{|{\Omega }_h|}{{\pi }^2\epsilon }M-|{\Omega }_h|v_{\partial \Omega }(x^{\ast },x^{\ast })+O(\alpha ),\end{array}$$

$$C_{\epsilon }=\frac{|{\Omega }_h|}{\pi \alpha {\epsilon }^2}+\frac{|{\Omega }_h|}{{\pi }^2\epsilon }M+|{\Omega }_h|v_{\partial \Omega }(x^{\ast },x^{\ast })+\frac{\beta }{\alpha }-g(x^{\ast })+O(\alpha ).$$

$${\tilde{\varphi }}_{\epsilon }(x)=-\frac{|{\Omega }_h|}{\pi \epsilon }-\frac{\alpha |{\Omega }_h|}{{\pi }^2}(M-\pi L[1])+O({\alpha }^2\epsilon ),$$

Full text

1 Introduction

The narrow escape problem (NEP) in diffusion theory, which goes back to Lord Rayleigh [18], is to calculate the mean first passage time (MFPT) of a Brownian particle to a small absorbing window on the otherwise reflecting boundary of a bounded domain. NEP has recently attracted significant attention from the point of view of mathematical and numerical modeling due to its relevance in molecular biology and biophysics. The small absorbing window often represents a small target on a cellular membrane, such as a protein channel, which is a target for ions [8], a receptor for neurotransmitter molecules in a neuronal synapse [5], a narrow neck in the neuronal spine, which is a target for calcium ions [13], and so on. A main concern for NEP is to derive an asymptotic expansion of the MFPT when the size of the small absorbing window tends to zero. There have been several significant works deriving the leading-order and higher-order terms of the asymptotic expansions of MFPT for regular and singular domains in two and three dimensions [1, 2, 3, 6, 9, 10, 11, 19, 20].

In this paper we focus on the MFPT of calcium ion in the three-dimensional dendritic spine shape domain. For many neurons in the mammalian brain, the postsynaptic terminal of an excitatory synapse is found in a specialized structure protruding from the dendritic draft, known as a dendritic spine (Figure 1(a)). Dendritic spines function as biochemical compartments that regulate the duration and spread of postsynaptic calcium fluxes produced by glutamatergic neurotransmission at the synapses [16]. Calcium ions are ubiquitous signaling molecules that accumulate in the cytoplasm in response to diverse classes of stimuli and, in turn, regulate many aspects of cell function [7]. The rate of diffusional escape from spines through their narrow neck, which we call MFPT, is one factor that regulates the retention time of calcium ions in dendritic spines.

Most spines have a bulbous head, and a thin neck that connects the head of the spine to the shaft of the dendrite (Figure 1). Spines are separated from their parent dendrites by this thin neck and compartmentalize calcium during synaptic stimulation. It has been given that many organelles inside the spine head do not affect the nature of the random motion of ions, mainly due to their large size relative to that of ions [12]. So we assume the environment inside the spine is homogenous, which means the movement of the calcium ions is pure diffusion.

The NEP can be mathematically formulated in the following way. Let $\Omega_{h}$ be a bounded simply connected domain in $\mathbb{R}^{3}$. Let $\Omega_{n}$ be a cylinder with length $L$ and radius $r=O(\varepsilon)$ which is much smaller than the length. Connecting this cylinder with $\Omega_{h}$, we have the geometry for the spine (Figure 1(b)). The connection part between $\Omega_{h}$ and $\Omega_{n}$ is a small interface which is denoted by $\Gamma_{\varepsilon}$. Let $\Omega=\Omega_{h}\cup\Omega_{n}$ denote the domain for the whole spine. Suppose that the boundary $\partial\Omega$ is decomposed into the reflecting part $\partial\Omega_{r}$ and the absorbing part $\partial\Omega_{a}$, where $\partial\Omega_{a}$ is the end of the thin cylinder neck. We assume that the area of $\partial\Omega_{a}$, $|\partial\Omega_{a}|=O(\varepsilon^{2})$ is much smaller than the area of the whole boundary. The NEP is to calculate the MFPT $u$ which is the unique solution to the following boundary value problem,

[TABLE]

where $\nu$ is the outer unit normal to $\partial\Omega$. The asymptotic analysis for NEP arises in deriving the asymptotic expansion of $u$ as $\varepsilon\rightarrow 0$, from which one can estimate the escape time of the calcium ions. In [15, 4], the first order of the asymptotic expansion has been obtained numerically as

[TABLE]

where $|\Omega_{h}|$ denotes the volume of the spine head. In this study, we derive higher order asymptotic solution to (1) by means of the Neumann-Robin model which is proposed in [14] to deal with the narrow escape time in a two-dimensional analogue domain.

In the Robin-Neumann model the solution to the original boundary value problem (1) in the singular domain $\Omega$ is approximated by the solution to the following boundary value problem in the smooth domain $\Omega_{h}$:

[TABLE]

where $\Omega_{h}$ is the spine head of $\Omega$ mentioned in Figure 1, and $\Gamma_{\varepsilon}:=\partial\Omega_{h}\setminus\overline{\partial\Omega_{r}}$ is the connection part between the spine head and the spine neck. Here $\alpha>0$ and $\beta>0$ are constants to be determined. We assume $\alpha<\alpha_{0}$ for some a priori constant $\alpha_{0}>0$ and $\varepsilon$ is sufficiently small so that $\alpha\varepsilon\ll 1$. We shall apply the layer potential technique to derive the asymptotic solution $u_{\varepsilon}$ to (3) as

[TABLE]

for $x\in\Omega_{h}$ and away from $\Gamma_{\varepsilon}$, where $M$ is a computable constant, $x^{*}$ is a fixed point in $\Gamma_{\varepsilon}$.

This study is organized as follows. In section 2, we review the Neumann function for the Laplacian in $\mathbb{R}^{3}$, which is a major tool for our study. In section 3, we derive the asymptotic solution for Robin-Neumann model. In section 4, we apply the Robin-Neumann boundary model to approximate the MFPT of calcium ion in dendritic spine. Numerical experiments are also given in this section to confirm the theoretical results. This study ends with a short conclusion in section 5.

2 Neumann function in $\mathbb{R}^{3}$

Let $\Omega$ be a bounded domain in $\mathbb{R}^{3}$ with $C^{2}$ smooth boundary $\partial\Omega$, and let $N(x,z)$ be the Neumann function for $-\Delta$ in $\Omega$ with a given $z\in\Omega$. That is, $N(x,z)$ is the solution to the boundary value problem

[TABLE]

where $\nu$ is the outer unit normal to the boundary $\partial\Omega$.

If $z\in\Omega$, then $N(x,z)$ can be written in the form

[TABLE]

where $R_{\Omega}(x,z)$ has weaker singularity than $1/|x-z|$ and solves the boundary value problem

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the inner product in $\mathbb{R}^{3}$.

If $z\in\partial\Omega$, then Neumann function on the boundary is denoted by $N_{\partial\Omega}$ and can be written as

[TABLE]

where $R_{\partial\Omega}(x,z)$ has weaker singularity than $1/|x-z|$ and solves the boundary value problem

[TABLE]

The structure of $R_{\partial\Omega}$ is given in [17] as

[TABLE]

where $z\in\partial\Omega$, $x\in\Omega\cup\partial\Omega$, where $H(z)$ denotes the mean curvature of $\partial\Omega$ at $z$, and $v_{\partial\Omega}$ is a bounded function.

3 Derivation of the asymptotic expansion

The goal in this section is to derive the asymptotic expansion of $u_{\varepsilon}$ to (3) as $\varepsilon\rightarrow 0$. For simplicity we assume the connection part $\Gamma_{\varepsilon}$ lies in a plane. The general case where $\Gamma_{\varepsilon}$ is curved can be handled with minor modifications.

Theorem 3.1**.**

The solution $u_{\varepsilon}$ to the boundary value problem (3) has the following asymptotic expansion,

[TABLE]

where $M$ is a constant given by

[TABLE]

$x^{*}$* is a fixed point in $\Gamma_{\varepsilon}$, and $\Phi(x)$ is a bounded function depending only on $\Omega_{h}$. The remainder $O(\alpha)$ is uniform in $x\in\Omega_{h}$ satisfying dist$(x,\Gamma_{\varepsilon})\geq c$ for some constant $c>0$.*

Proof. By integrating the first equation in (3) over $\Omega_{h}$ using the divergence theorem we get the compatibility condition

[TABLE]

Let us define $g(x)$ by

[TABLE]

which is seen to solve the boundary value problem

[TABLE]

Applying the Green’s formula and using (3) and (8), we obtain

[TABLE]

where

[TABLE]

Let $x\in\Gamma_{\varepsilon}$. Substitute the Robin boundary condition and the structure of Neumann function (5) into (9), we obtain

[TABLE]

where $\phi_{\varepsilon}(x)=\partial u_{\varepsilon}(x)/\partial\nu_{x}$. Note that the mean curvature $H(z)=0$ for $z\in\Gamma_{\varepsilon}$ since $\Gamma_{\varepsilon}$ is assumed to be flat.

By a simple change of variables, the above equation can be written as

[TABLE]

where $\Gamma_{1}=\{x/\varepsilon:x\in\Gamma_{\varepsilon}\}$, and $\tilde{\phi}_{\varepsilon}(x)=\varepsilon\phi_{\varepsilon}(\varepsilon x)$, $x\in\Gamma_{1}$.

Define two integral operators $L,~{}L_{1}:L^{\infty}(\Gamma_{1})\rightarrow L^{\infty}(\Gamma_{1})$ by

[TABLE]

Since $v_{\partial\Omega}$ is bounded, one can easily see that $L_{1}$ is bounded independently of $\varepsilon$. The integral operator $L$ is a also bounded (see the proof in Appendix A).

So we can write (10) as

[TABLE]

Collecting $\tilde{\phi}_{\varepsilon}$ terms, we have

[TABLE]

Assume here $\alpha<\alpha_{0}$ and $\alpha\varepsilon\ll 1$. It is easy to see that

[TABLE]

Noting that $g$ is $C^{1}$ in $\bar{\Omega}_{h}$, and $g(\varepsilon x)=g(x^{*})(1+O(\varepsilon))$ on $\Gamma_{1}$, we have from (11) that

[TABLE]

where

[TABLE]

By the compatibility condition (7), we can see that $\tilde{C}_{\varepsilon}=O((\alpha\varepsilon^{2})^{-1})$. Then collecting terms we have

[TABLE]

Plug (13) into the compatibility condition (7), we obtain

[TABLE]

which implies

[TABLE]

where $M=\int_{\Gamma_{1}}L[1]dx$, and $\int_{\Gamma_{1}}L_{1}[1](x)=\pi^{2}v_{\partial\Omega}(x^{*},x^{*})+O(\varepsilon)$. Hence from (12), we have

[TABLE]

Substitute $\tilde{C_{\varepsilon}}$ into (13), we have

[TABLE]

and hence

[TABLE]

In order to obtain the solution to (3), it remains to calculate the second term in (9). Combining (4), (5) and (15) yields

[TABLE]

provided that ${\rm dist}(x,\Gamma_{\varepsilon})\geq c$ for some constant $c>0$.

Finally, combining(9), (14) and (16), we obtain

[TABLE]

where

[TABLE]

is a bounded function depending only on $\Omega_{h}$.

4 Application to the narrow escape problem in dendritic spine

In this section we use the asymptotic solution to the Robin-Neumann model to approximate the calcium ion diffusion time in a dendritic spine domain. Without loss of generality, let $\Omega_{n}$ be placed along the $x_{1}$-axis with $\Gamma_{\varepsilon}$ at $x_{1}=0$ and $\partial\Omega_{a}$ at $x_{1}=L$. Since $\varepsilon\ll 1$, we assume $u(x)$ is constant in each cross section of $\Omega_{n}$ and thus is a function of $x_{1}$ only. The three-dimensional problem (1) restricted in $\Omega_{n}$ is then approximated by the one-dimensional problem

[TABLE]

By direct calculation we obtain the solution to the above problem as

[TABLE]

where $C$ is a constant. Solution $u$ satisfies the Robin boundary condition at $x=0$. Evaluating $u$ and $du/dx$ at $x_{1}=0$ yields the Robin condition

[TABLE]

By the continuity of $u$ and $\partial u/\partial\nu$ on $\Gamma_{\varepsilon}$, we obtain the Robin-Neumann boundary value problem in $\Omega_{h}$:

[TABLE]

where $\alpha=1/L$ and $\beta=L/2$. Applying Theorem 3.1 we obtain the asymptotic expansion of $u_{\varepsilon}$ as

[TABLE]

Note that the leading term coincides with (2), which is obtained in [15, 4] using numerical simulation.

4.1 Numerical experiments

In the rest of this section we shall conduct numerical experiments to verify the asymptotic expansion (17). We shall compare the asymptotic solution (17) with the solution $u$ to the original problem (1) obtained numerically with the finite element method. We shall also confirm the coefficients in the first two terms in (17).

For simplicity we confine ourself to the case when the connection part $\Gamma_{\varepsilon}$ is a disk of radius $\varepsilon$. In this case the constant $M$ in $\eqref{eq:2}$ has an explicit and elegant value $M=8/3$ (see Appendix B). The spine neck $\Omega_{n}$ is chosen to be a cylinder with radius $\varepsilon$ and length $L$, and whose axis is parallel to the normal of $\Gamma_{\varepsilon}$. To nondimensionalize our problem, the numerical results are regarded as using consistent units throughout this section.

For the first experiment, we select the spine head as an unit ball and set $\varepsilon=0.1,L=1.0$. In Figure 2(a) we plot the numerical solution of $u$ on $\partial\Omega_{h}$, i.e. the whole surface of the spine. Observe that the MFPT $u$ is relatively large in the spine head and decreases monotonically to zero towards the end of the spine neck. This is consistent with our intuition about the underlying physical process. In Figure 2(b) we plot the asymptotic solution $u_{\varepsilon}$ according to (17) on $\partial\Omega_{h}$, i.e. the surface of the spine head. Note that $u$ is relatively constant but is smaller for points closer to the connection part, which is also consistent with the physical intuition. In Figure 2(c) we plot the relative error between $u$ and $u_{\varepsilon}$, computed as $(u_{\varepsilon}-u)/u$, on the surface of the spine head. The maximal relative error is seen to be about $0.2\%$.

For the next experiment, we fix the neck length $L=1.0$ and let the neck radius $\varepsilon$ decreases from $0.1$ to $0.01$ in a step size of $0.01$. For each value of $\varepsilon$, we compare the value of the numerical solution $u$ to the original problem (1), the numerical solution of the Robin-Neumann model (3), which is denoted by $u_{r}$, and the asymptotic solution $u_{\varepsilon}$ given in (17). In Table 1 we list the value of $u,u_{r},u_{\varepsilon}$, as well as the relative error $(u_{\varepsilon}-u)/u$, at the center of the spine head. Clearly we have a good match of the solutions and small relative error for all values of $\varepsilon$.

The asymptotic solution $u_{\varepsilon}$ in (17) at a fixed $x$ can be considered as a quadratic polynomial of $1/\varepsilon$ with leading coefficients $|\Omega_{h}|L/\pi=4/3$ and $|\Omega_{h}|M/\pi^{2}\approx 1.13$. We now confirm these coefficients by fitting the value of $u$ by a quadratic polynomial of $1/\varepsilon$ in Table 1. The result is plotted in Figure 3. Clearly the coefficients of the fitting polynomial matches well with those of the asymptotic solution.

Another important parameter in the asymptotic solution (17) is $L$, the length of the spine neck. For the next experiment, we fix the neck radius $\varepsilon=0.05$ and let the neck radius $L$ increases from $1.0$ to $10.0$ in a step size of $1.0$. In Table 2 we list the value of $u,u_{r},u_{\varepsilon}$, as well as the relative error $(u_{\varepsilon}-u)/u$, at the center of the spine head. Clearly we have a good match of the solutions and small relative error for each value of $L$.

Finally we conduct numerical experiments on three different shapes of the spine head, which may correspond to different types of spine. The numerical solution $u$, the asymptotic solution $u_{\varepsilon}$ and the relative error between them are shown in Figure 4. The first row shows the results when the spine head is a relatively flat three dimensional domain. The second one is a relatively thin domain and the third one is a non-convex domain. From the relative error $(u_{\varepsilon}-u)/u$, which is described in the third column, we can see that the relative error is small enough to show that asymptotic formula $u_{\varepsilon}$ is a good approximation to the MFPT.

5 Conclusion

In this study, we used the Robin-Neumann model to compute the NEP in a three-dimensional domain with a long neck, which is often referred as a dendritic spine domain. This is a following up paper of [14], where the Robin-Neumann model is presented to solve the NEP in a two-dimensional analogue of a dendritic spine domain. In this paper, we derived asymptotic expansion formula for three-dimensional Robin-Neumann model. Our results demonstrate that this asymptotic expansion formula could approximate the MFPT up to at least second leading order using this model, which has not been reported previously.

The work of Hyundae Lee was supported by the National Research Foundation of Korea grant (NRF-2015R1D1A1A01059357). The work of Yuliang Wang was supported by the Hong Kong RGC grant (No. 12328516) and the NSF of China (No. 11601459).

6 Appendix

6.1 Appendix A.

Lemma 6.1**.**

(Single layer potential on an surface) Let $\Gamma_{1}\subset\mathbb{R}^{2}$ be a bounded domain with $\sup_{x,z\in\Gamma_{1}}|x-z|\leq 2$. The integral operator $L$: $L^{\infty}(\Gamma_{1})\mapsto L^{\infty}(\Gamma_{1})$ defined by

[TABLE]

is bounded, i.e.

[TABLE]

where $C$ is a constant independent of $x$.

Proof.

Let $B_{x,2}\subset\mathbb{R}^{2}$ denote the disk centered at $x$ with radius $2$. For given $x\in\Gamma_{1}$ we have

[TABLE]

Hence

[TABLE]

Therefore $L$ is bounded.

∎

6.2 Appendix B.

The value of

[TABLE]

is $\frac{16}{3}\pi$, where $\Gamma_{1}$ is unit disk.

Proof.

First, fix point $y$, draw a circle centered at $y$ with radius $r$, where the distance between $x$ and $y$ is $r$. Let $s$ be the distance between [math] and $y$.

Then we have

[TABLE]

When $r\leq 1-s$, by simple calculation, we can calculate the first term of (18) as

[TABLE]

Changing the order of $ds$ and $dr$, the second term of (18) becomes

[TABLE]

By calculations, we get the value of the first term of (19) as $\frac{4}{3}+\frac{\pi}{6}-\frac{3}{4}\sqrt{3}$, and second term $-\frac{1\pi}{3}+\frac{3}{4}\sqrt{3}$.

Finally, we get

[TABLE]

∎