# Global Optimization of the Mean First Passage Time for Narrow Capture
Problems in Elliptic Domains

Jason Gilbert Electronic mail<EMAIL_ADDRESS>Department of Mathematics
and Statistics, University of Saskatchewan  Alexei Cheviakov Corresponding
Author. Alternative English spelling: Alexey Shevyakov. Electronic mail:
<EMAIL_ADDRESS>Department of Mathematics and Statistics, University
of Saskatchewan

###### Abstract

Narrow escape and narrow capture problems which describe the average times
required to stop the motion of a randomly travelling particle within a domain
have applications in various areas of science. While for general domains, it
is known how the escape time decreases with the increase of the trap sizes,
for some specific 2D and 3D domains, higher-order asymptotic formulas have
been established, providing the dependence of the escape time on the sizes and
locations of the traps. Such results allow the use of global optimization to
seek trap arrangements that minimize average escape times. In a recent paper,
the escape time expansion for a 2D elliptic domain was derived, providing the
dependence of the average MFPT on sizes and locations of small internal traps.
The goal of this work is to systematically seek global minima of MFPT for
$1\leq N\leq 50$ traps, and compare the corresponding putative optimal trap
arrangements for different values of the domain eccentricity.

##  1 Introduction

The narrow capture problem, as described here, concerns the average time
required for a particle undergoing Brownian motion to encounter a region
within the domain, referred to as a trap, which causes its motion to cease. It
is mathematically defined as a Neumann-Dirichlet Poisson problem

$\begin{split}&\Delta u=-\dfrac{1}{D}\ ,\quad x\in\bar{\Omega}\
;\qquad\bar{\Omega}=\Omega\setminus\mathop{\cup}_{j=1}^{N}\Omega_{\varepsilon_{j}}\
;\\\\[5.0pt] &\partial_{n}u=0\ ,\quad x\in\partial\Omega\ ;\qquad u=0\ ,\quad
x\in\partial\Omega_{\varepsilon_{j}}\ ,\quad j=1,\ldots,N\ ;\end{split}$ (1.1)

where $\Omega\subset\mathbb{R}^{n}$, $n=2,3$, denotes the physical domain of
the problem; $\\{\Omega_{\varepsilon_{j}}\\}_{j=1}^{N}$ are small volume trap
domains within $\Omega$; $\bar{\Omega}$ is the domain except the traps, where
the motion of particles takes place, and $\partial_{n}$ denotes the normal
derivative on the boundary of the domain. The diffusivity $D$ of the medium
filling $\bar{\Omega}$ is assumed constant. Problem (1.1) describes the
distribution of the mean capture time, the time $u(x)$ needed for a particle
to be captured by any trap, averaged over a large number of launches from the
same point $x\in\bar{\Omega}$. An illustration of the problem is provided in
Figure 1.

Figure 1:  (a) A two-dimensional narrow capture problem in the unit disk
containing internal traps with absorbing boundaries
$\\{\partial\Omega_{\epsilon_{j}}\\}$. (b) A three-dimensional narrow capture
problem, a sample Brownian particle trajectory, leading to a capture in a trap
(lowermost) denoted by purple color (color online).

The boundary conditions on $\partial\Omega$ are reflective: $\partial_{n}u=0$,
whereas the traps $\Omega_{\varepsilon_{j}}$ are characterized by immediate
absorption of a boundary particle, which is manifested by a Dirichlet boundary
conditions $u=0$ on all $\partial\Omega_{\varepsilon_{j}}$.

The above generic formulation affords a variety of physical interpretations,
ranging from biological to electrostatic (see, e.g., Ref. [1, 2] for an
overview of applications). In this work, it will be strictly considered in
terms of a particle undergoing Brownian motion [3]. In this case, the problem
regards the stopping time [1, 4] of a Brownian particle, confined to a domain
with a reflective boundary which contains a number of absorbing “traps”. When
the path of the particle intersects the boundary of one of these traps, the
particle is captured, meaning that the process of Brownian motion is stopped.
This stopping time may be interpreted as the time required for the particle to
leave the confining domain, thus it is often referred to as the first passage
time [5, 6, 7]. As Brownian motion is an inherently random process, there is
no single first passage time which accurately characterizes the process;
instead, the mean first passage time (MFPT), representing the average time
taken for a particle to encounter a trap based on its first observed position,
is used. Interpreting the problem (1.1) accordingly, $u$ is the MFPT; $D$ is
the diffusion coefficient, representing the mobility of the Brownian particle;
$\Omega_{\varepsilon_{j}}$ is the portion of the domain occupied by the
$j_{th}$ trap.

Given the physical domain and the number and sizes of the traps, it is of
interest to ask whether there is an _optimal_ arrangement of traps within the
domain which minimizes the MFPT, or in other words, maximizes the rate at
which Brownian particles encounter the absorbing traps. Related work dedicated
to similar optimization, in the case that the traps are asymptotically small
relative to the size of the domain, for various kinds of confining domains
with interior or boundary traps, can be found, for example, in Refs [8, 9, 10,
11, 12] and references therein. Both (putative) globally optimal and multiple
locally optimal arrangements of boundary and volume traps have been computed
in various settings.

In this contribution, we consider the narrow capture problem for the case of a
2D elliptical domain, elaborating on previous work on the subject for the case
of a circular domain [8], and the more general case of the elliptical domain
[9]. The paper is organized as follows. In Section 2 we briefly summarise the
previous results which this work is based on. This includes a review of the
approximations used in the case that the traps are small relative to the
domain size, and are well-separated; as well as an explanation of the choice
of merit function, and relevant constraints, for each case.

In Section 3 the optimization method chosen in our study is described. This
includes an explanation of our selection of algorithms, as well the details of
their use, and some decisions made to facilitate comparison to previous
studies. Specifically, we seek the optimal configuration of traps for numbers
of traps $N\leq 50$, and elliptic domain eccentricities of $0$, $1/8$, $1/4$,
and $1/2$.

In Section 4 the results of the study are presented. These results include the
optimized values of the merit functions (related to the domain-average
Brownian particle capture time) for each number of traps, and each domain
eccentricity; in the case of the unit disk, our results are compared to those
of a previous study. Additionally, the distributions of traps for some
illustrative cases are shown, and bulk measures of trap distribution are
calculated for each of the optimized configurations.

In Section 5, the results presented in Section 4 are discussed, and some
remarks are provided regarding on the distribution of traps, the
interpretation of the bulk quantities used, and some related problems.

##  2 Problem Statement

The goal of this contribution is to further explore the optimal configuration
of absorbing traps in terms of the asymptotic expansions for the circular and
elliptical domains for which asymptotic MFPT formulae depending on trap
positions are now available [8, 9]. To this end, numerical methods will be
used to approximate the optimum configurations of larger numbers of traps than
have previously been considered. In the case of the unit disk, the parameter
space used in this study is general and does not assume any simplifying
restrictions that have been imposed in previous works to reduce computation
complexity.

In essence, the problem at hand is to find the trap positions which minimize
the spatial average of the MFPT $u(x)$, denoted by $\overline{u}_{0}$. This
section will summarize the equations used to minimize this quantity, as
derived in [8, 9]. Here the unit disk will be considered a special case of the
elliptical domain with zero eccentricity. In either case, when the domain
contains $N$ identical circular traps, each of radius $\varepsilon$, the
average MFPT satisfies [9]

$\overline{u}_{0}\ =\ \dfrac{|\Omega|}{2\pi D\nu N}\ +\
\dfrac{2\pi}{N}\textbf{e}^{T}\mathcal{G}\mathcal{A}\ ,$ (2.1)

where the column vector $\mathcal{A}$ is the solution of the linear system

$\left[I+2\pi\nu\left(I-E\right)\mathcal{G}\right]\mathcal{A}\ =\
\dfrac{|\Omega|}{2\pi DN}\textbf{e}\ .$

Here $E\equiv\textbf{e}\textbf{e}^{T}/N$, $\textbf{e}=(1,\ldots,1)^{T}$,
$\nu\equiv-1/\log\varepsilon$, and the Green’s matrix $\mathcal{G}$ depends on
the trap locations $\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{N}\\}$, such that

$\mathcal{G}_{ij}\ =\
\left\\{\begin{array}[]{ll}G(\mathbf{x}_{i};\mathbf{x}_{j})\quad\text{if}\
i\neq j\ ,\\\\[5.0pt] R(\mathbf{x}_{i};\mathbf{x}_{j})\quad\text{if}\ i=j\
;\end{array}\quad\text{where}\quad
R(\mathbf{x}_{i};\mathbf{x}_{j})=\lim_{\mathbf{x}_{i}\to\mathbf{x}_{j}}G(\mathbf{x}_{i};\mathbf{x}_{j})\right.$
(2.2)

where $G$ is the Green’s function of the domain.

Examining the equation (2.1), it can be seen that the first term depends only
on the combined trap size but does not depend on the trap locations. The
second term explicitly depends on the trap locations is the quantity and
therefore can be optimized by adjusting trap positions. The merit function
subject to optimization can be chosen as

$q\left(\mathbf{x}\right)\ =\ \textbf{e}^{T}\mathcal{G}\mathcal{A}\ .$ (2.3)

For a value of $q$ to be a permissible solution, it is required that
$\overline{u}_{0}\geq 0$, as it is a measure of time; all traps must be within
the domain; no trap may be in contact with any other trap (or the two must
instead be represented as a single non-circular trap).

While the preceding statements are true for both the circular and the
elliptical domain, the elements of the matrix $\mathcal{G}$ are populated by
evaluating the Green’s function of the domain for each pair of traps, and the
form of this function differs for the two cases considered here.

In the case of the circular domain, the elements of the Green’s matrix
$\mathcal{G}$ are computed using the Neumann Green’s function [8]

$G(\mathbf{x};\mathbf{x}_{0})\ =\
\dfrac{1}{2\pi}\left(-\log|\mathbf{x}-\mathbf{x}_{0}|\ -\
\log\left|\mathbf{x}|\mathbf{x}_{0}|-\dfrac{\mathbf{x}_{0}}{|\mathbf{x}_{0}|}\right|\
+\ \dfrac{1}{2}\left(|\mathbf{x}|^{2}+|\mathbf{x}_{0}|^{2}\right)\ -\
\dfrac{3}{4}\right)\ ,$ (2.4a) and its regular part $R(\mathbf{x})\ =\
\dfrac{1}{2\pi}\left(-\log\left|\mathbf{x}|\mathbf{x}|-\dfrac{\mathbf{x}}{|\mathbf{x}|}\right|\
+\ |\mathbf{x}|^{2}\ -\ \dfrac{3}{4}\right)\ .$ (2.4b)

The Green’s function for the elliptical domain, in the form of a quickly
convergent series, has been derived in [9] using elliptical coordinates
$(\xi,\eta)$. It has the form

$\begin{split}G(\mathbf{x},\mathbf{x}_{0})\ =\
\dfrac{1}{4|\Omega|}\left(|\mathbf{x}|^{2}+|\mathbf{x}_{0}|^{2}\right)\ -\
\dfrac{3}{16|\Omega|}(a^{2}+b^{2})\ -\ \dfrac{1}{4\pi}\log\beta\ -\
\dfrac{1}{2\pi}\xi_{>}\\\\[5.0pt] \ -\
\dfrac{1}{2\pi}\sum^{\infty}_{n=0}\log\left(\prod_{j=1}^{8}|1-\beta^{2n}z_{j}|\right)\
,\qquad\mathbf{x}\neq\mathbf{x}_{0}\ ,\end{split}$ (2.5)

where $a$ and $b$ are the major and minor axis of the domain, respectively;
the area of the domain is $|\Omega|=\pi ab$, the parameter
$\beta=(a-b)/(a+b)$, and the values $\xi_{>}=\mathop{\hbox{\rm
max}}(\xi,\xi_{0})$, $z_{1},\ldots,z_{8}$ are defined in terms of the
elliptical coordinates $\xi$ and $\eta$ as follows.

The Cartesian coordinates $(x,y)$ and the elliptical coordinates $(\xi,\eta)$
are related by the transformation

$x=f\cosh\xi\cos\eta\ ,\quad y=f\sinh\xi\sin\eta\ ,\quad f=\sqrt{a^{2}-b^{2}}\
.$ (2.6a) For the major and minor axis of the elliptical domain, one has
$a=f\cosh\xi_{b}\ ,\quad b=f\sinh\xi_{b}\
,\quad\xi_{b}=\tanh^{-1}\left(\dfrac{b}{a}\right)=-\dfrac{1}{2}\log\beta\ .$
(2.6b)

For the backward transformation, to determine $\xi(x,y)$, equation (2.6a) is
solved to give

$\xi=\dfrac{1}{2}\log\left(1-2s+2\sqrt{s^{2}-s}\right)\ ,\quad
s\equiv\dfrac{-\mu-\sqrt{\mu^{2}+4f^{2}y^{2}}}{2f^{2}}\ ,\quad\mu\equiv
x^{2}+y^{2}-f^{2}\ .$ (2.7a) In a similar fashion, $\eta(x,y)$ is found in
terms of $\eta_{\star}\equiv\sin^{-1}(\sqrt{p})$ to be $\eta\ =\
\left\\{\begin{array}[]{ll}\eta_{\star}\ ,&\text{if}\ x\geq 0,\ y\geq
0\\\\[5.0pt] \pi-\eta_{\star}\ ,&\text{if}\ x<0,\ y\geq 0\\\\[5.0pt]
\pi+\eta_{\star}\ ,&\text{if}\ x\leq 0,\ y<0\\\\[5.0pt] 2\pi-\eta_{\star}\
,&\text{if}\ x>0,\ y<0\\\\[5.0pt] \end{array}\right.\ ,\quad\text{where}\quad
p\equiv\dfrac{-\mu+\sqrt{\mu^{2}+4f^{2}y^{2}}}{2f^{2}}\ .$ (2.7b)

Finally, the $z_{j}$-terms appearing in the infinite sum of equation (2.5) are
defined via $\xi$, $\eta$, and $\xi_{b}$ as

$\begin{array}[]{lll}z_{1}\equiv e^{-|\xi-\xi_{0}|+i(\eta-\eta_{0})}\
,&z_{2}\equiv e^{|\xi-\xi_{0}|-4\xi_{b}+i(\eta-\eta_{0})}\ ,&z_{3}\equiv
e^{-(\xi+\xi_{0})-2\xi_{b}+i(\eta-\eta_{0})}\ ,\\\\[4.0pt] z_{4}\equiv
e^{(\xi+\xi_{0})-2\xi_{b}+i(\eta-\eta_{0})}\ ,&z_{5}\equiv
e^{(\xi+\xi_{0})-4\xi_{b}+i(\eta+\eta_{0})}\ ,&z_{6}\equiv
e^{-(\xi+\xi_{0})+i(\eta+\eta_{0})}\ ,\\\\[4.0pt] z_{7}\equiv
e^{|\xi-\xi_{0}|-2\xi_{b}+i(\eta+\eta_{0})}\ ,&z_{8}\equiv
e^{-|\xi-\xi_{0}|-2\xi_{b}+i(\eta+\eta_{0})}\ .&\\\ \end{array}$ (2.8)

The regular part of the Neumann Green’s function, $R$, can be expressed in
similar terms as $G$ in equation (2.5) but requires a restatement of the
$z_{j}$-terms given in (2.8). It is given by

$\begin{split}R(\mathbf{x}_{0})\ =\ &\dfrac{|\mathbf{x}_{0}|^{2}}{2|\Omega|}\
-\ \dfrac{3}{16|\Omega|}(a^{2}+b^{2})\ +\ \dfrac{1}{2\pi}\log(a+b)\ -\
\dfrac{\xi_{0}}{2\pi}\ +\
\dfrac{1}{4\pi}\log\left(\cosh^{2}\xi_{0}-\cos^{2}\eta_{0}\right)\\\\[5.0pt]
&\ -\ \dfrac{1}{2\pi}\sum_{n=1}^{\infty}\log\left(1-\beta^{2n}\right)\ -\
\dfrac{1}{2\pi}\sum_{n=0}^{\infty}\log\left(\prod_{j=2}^{8}\left|1-\beta^{2n}z^{0}_{j}\right|\right).\end{split}\
$ (2.9a) Here $z^{0}_{j}$ denotes the limiting value of $z_{j}$, as defined in
equation (2.8), as $(\xi,\eta)\to(\xi_{0},\eta_{0})$, given by
$\begin{array}[]{lll}&z^{0}_{2}=\beta^{2}\ ,&z^{0}_{3}=\beta e^{-2\xi_{0}}\
,\\\\[5.0pt] z^{0}_{4}=\beta e^{2\xi_{0}}\
,&z^{0}_{5}=\beta^{2}e^{2(\xi_{0}+i\eta_{0})}\
,&z^{0}_{6}=e^{2(-\xi_{0}+i\eta_{0})}\ ,\\\\[5.0pt] z^{0}_{7}=\beta
e^{2i\eta_{0}}\ ,&z^{0}_{8}=\beta e^{2i\eta_{0}}\ .&\\\ \end{array}$ (2.9b)

##  3 Global optimization

In this section, the methods used to find the optimum trap configurations
minimizing the average MFPT (2.1) are discussed. We include a description of
the general strategy for optimization, the algorithms used, and their specific
implementations.

For both problems the same general approach was taken to finding the optimum.
This was to search iteratively, switching between global and local searches
after each iteration. The global search used the particle swarm algorithm
PSwarm [13], as implemented in the freely available software package OPTI
[14]. The default values for the social and cognitive parameters were chosen,
meaning the local optimum known by each particle tended to be as attractive as
the known global optimum. These values were chosen with the intent that the
parameter space would be explored as broadly as possible. For the local
search, the Nelder-Mead algorithm [15], as implemented in MATLAB R2020, was
used. The algorithms were chosen based on their generality and simplicity to
adapt them to the current study.

In addition, for the case of the elliptical domain, special care needed to be
taken to ensure that the traps were well-separated. If the traps were to come
into contact, or overlap with one another, the asymptotic equation (2.1) can
yield non-physical values $\overline{u}_{0}<0$, (which is a common feature of
asymptotic formulas that replace finite traps with “point traps” in various
narrow escape and narrow capture setups). In the MFPT optimization, the traps
are effectively repelled from one another, as well as from their “reflections”
in the domain boundary; this is mathematically manifested in the fact that the
Green’s function (2.4a) for the disk domain grows logarithmically as
$\mathbf{x}\to\mathbf{x}_{0}$, as well as when $|\mathbf{x}|\to 1$. In
particular, $q$ increases as distance between traps decreases, as traps begin
to overlap, $q$ decreases extremely rapidly, appearing to the optimization
algorithm to be a favourable configuration. Though this problem can be
addressed by artificially assigning $q$ a very high value when an unacceptable
configuration is encountered, this approach was found to significantly reduce
the effectiveness of the global search, as many evaluations of $q$ would be
wasted on these configurations. Instead, an optimum was first found using the
iterative approach taking $\varepsilon=0$, following which a local search was
carried out using these coordinates as an initial guess for a local search,
and taking $\varepsilon=0.05$ in order to facilitate comparison with previous
studies [9].

Defining the eccentricity of the elliptic domain according to the usual
formula

$\kappa\ =\ \sqrt{1-\left(\dfrac{b}{a}\right)^{2}}\ ,$ (3.1)

where $a$ is the major axis and $b$ the minor, optimum configurations were
computed for $N\leq 50$ for domain eccentricities of $\kappa=0$, $1/8$, $1/4$,
and $1/2$. For each eccentricity value, the axes of the ellipse were chosen so
that the area of the domain was fixed at $|\Omega|=\pi$, to allow for natural
comparisons.

##  4 Optimization Results

A comparison of $\overline{u}_{0}$ for the optimal configurations found for
the unit disk in previous work [8] to those presented here, found using the
more accurate approximation [9], shows that the newly computed optimums for
the unit disk are consistently better than previous results. Though it is
convention to discuss optimal configurations in terms of their interaction
energy, the quantity used as the merit function, here the $\overline{u}_{0}$
is used for comparison. This serves the purpose of normalizing the results
obtained using the expression for $\overline{u}_{0}$ used in previous work
[8],

$\overline{u}_{0}\ =\ \dfrac{|\Omega|}{2\pi D\nu
N}\left(1+\dfrac{2\pi\nu}{N}\textbf{e}^{T}\mathcal{G}\textbf{e}\right)\ ,$
(4.1)

to equation (2.1), the newly derived expression [9],

$\overline{u}_{0}\ =\ \dfrac{|\Omega|}{2\pi D\nu
N}\left(1+\dfrac{4\pi^{2}D\nu}{|\Omega|}\textbf{e}^{T}\mathcal{G}\mathcal{A}\right)\
.$ (4.2)

Though it would be sufficient to compute only the second term for comparison,
$\overline{u}_{0}$ serves as a more consistent standard for comparison. Table
1 compares $\overline{u}_{0}$ for each $N$ reported in the previous study,
computed using the results found in Table 2 of Ref. [8], while Figure 2
depicts the relative difference between the two.

The computed optimal values of the merit function (2.3) that correspond to
putative globally optimal minima of the average MFPT (2.1), for the domain
eccentricities $\kappa=0$ (circular disk), $1/8$, $1/4$, and $1/2$, are
presented in Table 2 below, and are graphically shown in Figure 3. While the
first three plots are nearly identical, the plot (d) for the largest
eccentricity value differs significantly for small $N$, but becomes similar to
the other plots for larger $N$.

Plots comparing the optimal configurations of select $N$ for each of the
eccentricities considered in this study are shown in Figures 4–7. Each plot
shows the position of the trap within the domain, along with a visualization
of a Delaunay triangulation [16] calculated using the traps as vertices, to
illustrate the distribution of, and relative distance between, traps. In
addition, it was of interest to see how the ring-like distribution of traps
would change with the eccentricity of the domain. To visualize this change, a
scaling factor was calculated for each trap such that each trap would lie on a
scaled copy of the elliptic domain boundary. These scaling factors are shown
in the lower subplots in Figures 4–7. For the case of $N=5$, Figure 4 can be
compared to the optimal configurations presented in Ref [9], through which it
can be seen that the two are qualitatively similar and exhibit the same
relationship between trap distribution and domain eccentricity.

In order to examine the distribution of traps in terms of their mutual
distance, a Delaunay triangulation was computed to identify approximate
nearest neighbours to each trap (see upper plots in Figures 4–7). In general,
for a configuration of $N$ traps distributed, in some sense, “uniformly” over
the elliptic domain of area $|\Omega|=\pi$, the average “area per trap” is
given by $A(N)=|\Omega|/N=\pi/N$. Likening an optimal arrangement of $N$ traps
to a collection of circles packed into an enclosed space, the (average)
distance $\langle d\rangle$ between two neighbouring traps would be the
distance between the centers of two identical circles representing the area
occupied by each trap; it would be related to the area per trap as
$A(N)=\pi\langle d\rangle^{2}/4$. One consequently finds that the average
distance between neighbouring traps, equivalent to the diameter of one of the
circles, is given by

$\langle d\rangle\ =\ \sqrt{\dfrac{4|\Omega|}{\pi N}}\ =\ \dfrac{2}{\sqrt{N}}\
,$ (4.3)

Extending this comparison to the traps nearest the boundary, the smallest
distance between a trap and the boundary was taken to be the radius of a
circle surrounding the trap, and the diameter of this circle was compared to
the smallest distance between two traps. This essentially provides a measure
of the distance between a near-the-boundary trap and its “reflection” in the
Neumann boundary.

In Figure 8, for each of the four considered eccentricities of the elliptic
domain, the mean pairwise distance between neighbouring traps is plotted as a
function of $N$, along with minimum pairwise distance between traps, and
$2\times$ minimal distance to the boundary. These are compared with the
average distance formula (4.3) coming from the “area per trap” argument. It
can be observed that the simple formula (4.3) may be used as a reasonable
estimate of common pairwise distances between traps in an optimal
configuration. The case of $\kappa=0.5$ serves as somewhat of an exception.
This could either be due to a significant change in the distribution of the
traps as the eccentricity increases, or it could be an artefact introduced by
the calculation of the Delaunay triangulation.

As mentioned in Section 3, the optimal configuration for the case
$\varepsilon=0$ was used as an initial guess when used to search for the
optimum in the case $\varepsilon=0.05$. Using this approach it was found that
the optimal configuration for the two cases were not identical, implying that
the optimal configuration significantly depends on the size of the traps.

$N$ | $\overline{u}_{0}$ | $\overline{u}_{0}^{\prime}$
---|---|---
6 | 0.11648 | 0.11648
7 | 0.09299 | 0.09297
8 | 0.07660 | 0.07660
9 | 0.06518 | 0.06512
10 | 0.05653 | 0.05624
11 | 0.04920 | 0.04900
12 | 0.04291 | 0.04278
13 | 0.03805 | 0.03796
14 | 0.03380 | 0.03375
15 | 0.03042 | 0.03038
16 | 0.02747 | 0.02745
17 | 0.02502 | 0.02499
18 | 0.02286 | 0.02280
19 | 0.02078 | 0.02076
20 | 0.01909 | 0.01907
21 | 0.01756 | 0.01755
22 | 0.01626 | 0.01624
23 | 0.01512 | 0.01510
24 | 0.01411 | 0.01403
25 | 0.01314 | 0.01307
Table 1: Average MFPT in the unit disk for previously computed optimal
configurations $\overline{u}_{0}$ ( Ref. [8], Table 2 ), compared to that of
the newly computed configurations for the unit disk
$\overline{u}_{0}^{\prime}$. Here it can be seen that the new values are
consistently smaller, at most differing in the third significant figure. A
plot of the difference between these two, relative to the previous results,
can be found in Figure 2.  Figure 2: Relative difference between average MFPT
in the unit disk for previously computed optimal configurations
$\overline{u}_{0}$, compared to that of the newly computed configurations
$\overline{u}_{0}^{\prime}$ according to
$(\overline{u}_{0}^{\prime}-\overline{u}_{0})/\overline{u}_{0}$. Here it can
be seen that the new values are consistently smaller, at most differing in the
third significant figure. A table of these values can be found in Table 1.

Figure 3: The putative optimal values of the merit function (2.3) for the
ellipse eccentricity values $\kappa=0$ (a), $\kappa=0.125$ (b), $\kappa=0.250$
(c), and $\kappa=0.500$ (d), as functions of the number of traps $N$. [The
corresponding numerical values are presented in Table 2.]

$N$ | Merit Value
---|---
| $\kappa=0$ | $\kappa=0.125$ | $\kappa=0.25$ | $\kappa=0.5$
1 | -0.0597 | -0.0594 | -0.0540 | 0.1730
2 | -0.0754 | -0.0792 | -0.0854 | 0.0175
3 | -0.0969 | -0.0967 | -0.0959 | -0.0452
4 | -0.1112 | -0.1113 | -0.1115 | -0.0793
5 | -0.1207 | -0.1207 | -0.1200 | -0.1007
6 | -0.1272 | -0.1274 | -0.1289 | -0.1154
7 | -0.1348 | -0.1347 | -0.1342 | -0.1261
8 | -0.1409 | -0.1408 | -0.1393 | -0.1343
9 | -0.1451 | -0.1451 | -0.1447 | -0.1407
10 | -0.1489 | -0.1494 | -0.1492 | -0.1457
11 | -0.1526 | -0.1532 | -0.1533 | -0.1498
12 | -0.1567 | -0.1566 | -0.1569 | -0.1530
13 | -0.1599 | -0.1598 | -0.1603 | -0.1559
14 | -0.1632 | -0.1632 | -0.1632 | -0.1587
15 | -0.1659 | -0.1660 | -0.1657 | -0.1614
16 | -0.1685 | -0.1686 | -0.1683 | -0.1642
17 | -0.1708 | -0.1705 | -0.1708 | -0.1668
18 | -0.1731 | -0.1731 | -0.1733 | -0.1693
19 | -0.1756 | -0.1755 | -0.1756 | -0.1718
20 | -0.1777 | -0.1776 | -0.1775 | -0.1741
21 | -0.1798 | -0.1797 | -0.1796 | -0.1756
22 | -0.1815 | -0.1815 | -0.1816 | -0.1768
23 | -0.1831 | -0.1831 | -0.1833 | -0.1800
24 | -0.1848 | -0.1851 | -0.1848 | -0.1820
25 | -0.1864 | -0.1867 | -0.1864 | -0.1834

$N$ | Merit Value
---|---
| $\kappa=0$ | $\kappa=0.125$ | $\kappa=0.25$ | $\kappa=0.5$
26 | -0.1880 | -0.1882 | -0.1880 | -0.1858
27 | -0.1897 | -0.1896 | -0.1896 | -0.1875
28 | -0.1911 | -0.1910 | -0.1911 | -0.1893
29 | -0.1925 | -0.1925 | -0.1925 | -0.1909
30 | -0.1940 | -0.1940 | -0.1938 | -0.1920
31 | -0.1953 | -0.1953 | -0.1952 | -0.1941
32 | -0.1964 | -0.1964 | -0.1965 | -0.1953
33 | -0.1977 | -0.1978 | -0.1978 | -0.1958
34 | -0.1989 | -0.1989 | -0.1990 | -0.1975
35 | -0.2000 | -0.2001 | -0.2000 | -0.1987
36 | -0.2012 | -0.2012 | -0.2012 | -0.2003
37 | -0.2025 | -0.2024 | -0.2023 | -0.2005
38 | -0.2035 | -0.2035 | -0.2033 | -0.2022
39 | -0.2045 | -0.2044 | -0.2043 | -0.2028
40 | -0.2056 | -0.2055 | -0.2053 | -0.2037
41 | -0.2065 | -0.2065 | -0.2064 | -0.2046
42 | -0.2074 | -0.2073 | -0.2073 | -0.2049
43 | -0.2083 | -0.2084 | -0.2083 | -0.2065
44 | -0.2093 | -0.2092 | -0.2092 | -0.2069
45 | -0.2102 | -0.2103 | -0.2102 | -0.2086
46 | -0.2110 | -0.2111 | -0.2111 | -0.2091
47 | -0.2119 | -0.2120 | -0.2119 | -0.2102
48 | -0.2127 | -0.2128 | -0.2128 | -0.2098
49 | -0.2136 | -0.2136 | -0.2136 | -0.2118
50 | -0.2144 | -0.2143 | -0.2143 | -0.2126

Table 2: Optimized values of the merit function (2.3), for each number of
traps $N$ and eccentricity $\kappa$ considered in this study. Plots of these
values are found in Figure 3.

Figure 4: Plots depicting optimal trap distribution for $N=5$, comparing
eccentricities of (a) $\kappa=0$, (b) $\kappa=0.125$, (c) $\kappa=0.250$, (d)
$\kappa=0.500$. Upper plots show positions of traps along with a crude
visualization of nearest-neighbour pairs calculated using Delaunay
triangulation. Lower plots show the scaling factor applied to the domain
boundary such that the scaled version would pass through each trap, numbered
from 1 to $N$.

Figure 5: Plots depicting optimal trap distribution for $N=10$, comparing
eccentricities of (a) $\kappa=0$, (b) $\kappa=0.125$, (c) $\kappa=0.250$, (d)
$\kappa=0.500$. Upper plots show positions of traps along with a crude
visualization of nearest-neighbour pairs calculated using Delaunay
triangulation. Lower plots show the scaling factor applied to the domain
boundary such that the scaled version would pass through a given trap,
numbered from 1 to $N$.

Figure 6: Plots depicting optimal trap distribution for $N=25$, comparing
eccentricities of (a) $\kappa=0$, (b) $\kappa=0.125$, (c) $\kappa=0.250$, (d)
$\kappa=0.500$. Upper plots show positions of traps along with a crude
visualization of nearest-neighbour pairs calculated using Delaunay
triangulation. Lower plots show the scaling factor applied to the domain
boundary such that the scaled version would pass through a given trap,
numbered from 1 to $N$.

Figure 7: Plots depicting optimal trap distribution for $N=40$, comparing
eccentricities of (a) $\kappa=0$, (b) $\kappa=0.125$, (c) $\kappa=0.250$, (d)
$\kappa=0.500$. Upper plots show positions of traps along with a crude
visualization of nearest-neighbour pairs calculated using Delaunay
triangulation. Lower plots show the scaling factor applied to the domain
boundary such that the scaled version would pass through a given trap,
numbered from 1 to $N$.

Figure 8: Plots depicting local properties of trap distribution, for domain
eccentricities $\kappa=0$ (a), $\kappa=0.125$ (b), $\kappa=0.250$ (c), and
$\kappa=0.500$ (d). The curve entitled “Measure of Area per Trap” shows the
distance $\langle d\rangle$ computed using the “area per trap” argument and
the resulting formula (4.3).

##  5 Discussion

At this point some interpretation of the previously stated results will be
presented. This discussion will concern the putative optimal values of the
average MFPT (2.1) in elliptic domains with internal traps, values of the
related merit function (2.3), the positions of traps within the domain, and
the bulk measures of trap distribution which were employed. To begin, the
method of study will be briefly reiterated. In order to study the dependence
of optimal trap configurations on both the number of traps, and the
eccentricity of the elliptic domain, merit functions for the average MFPT were
minimized for $N\leq 50$ and eccentricity (3.1) values of 0, 0.125, 0.25, and
0.5, while the area of the ellipse was kept constant, $|\Omega|=\pi$. In the
search for a global optimum, an iterative approach, which switched between
global and local searches, was used. The method used here was similar to that
used in Ref [9], though a different algorithm was used for the local search,
as well as in Ref [11] which used a different algorithm for both searches. A
somewhat different approach was employed in Ref [17], which made use of
numerical solutions to the Poisson problem.

In the case of the unit disk, comparing the results of this study to those of
a previous study [8] demonstrated that the optimums reported here are
consistently an improvement on previous work. This is due to the use of a more
accurate asymptotic expression [9], as well as removing the constraint that
all traps be located on rings within the domain.

Plots of the putative globally optimal merit function values as functions of
$N$, for each eccentricity, can be found in Figure 3. From these plots it
seems that eccentricity is an important factor when there are few traps,
around $N<10$, but each function behaves similarly as $N$ increases. In
particular, as the number of traps $N$ increases, it is natural to expect that
the average MFPT $\overline{u}_{0}$ (2.1) approaches zero; the merit function
$q(\mathbf{x})$ therefore must, as $N\to\infty$, approach from above the value
$|\Omega|/(4\pi^{2}D\nu)\simeq-0.238$, which agrees with the plots in Figure
3.

Examination of the positions of traps in the optimized configurations, both
visually and in terms of their radial coordinates, gives the impression that
the optimal configuration is one which consists of traps placed on the
vertices of nested polygons. These polygons, while irregular, seem to possess
some consistent structure, including being convex. (It is interesting to note
that the optimal configurations of confined interacting points often take
similar forms, both in two and three dimensions [18, 19, 20, 21].) Due to the
geometrical symmetries of the ellipse, the optimal configurations are defined
uniquely modulo the group $C^{2}\times C^{2}$ of reflections with respect to
both axes, which includes the rotation by $\pi$. The numerical algorithms,
however, choose a single specific representative of the equivalent putative
globally optimal configurations. For example, for non-symmetric numerically
optimal configuration, several traps may be found along the midline of one
half (right or left) of the domain. Optimal trap configurations with the same
symmetry group as the ellipse were also observed (see, e.g., Figure 5(c)).

In addition to the examination of individual trap positions in each optimized
configuration, quantities were calculated using the distances between
neighbouring traps, defined according to a Delaunay triangulation of the trap
coordinates, which served as bulk measures of the distribution of traps in
each configuration. Plots of these measures, shown in Figure 8, illustrate
that the mean distance between neighbouring traps tends to be close to the
diameter of a circle which would occupy the average area of the domain per
trap, as in equation (4.3). Additionally, the minimum distance between any two
traps tends to be twice the minimum distance between a trap and the domain
boundary, which supports the intuitive reasoning that for the boundary value
problem (1.1) with interior traps, the Neumann boundary condition on
$\partial\Omega$ “reflects” each trap, so that under the average MFPT
optimization, every trap tends to “repel” from its reflection in the boundary
the same way as it is repelled from other traps.

In future work, it would be of interest to address two related problems. The
first is to carry out similar investigations for the near-disk domains
considered in Ref. [9]. Another interesting direction is the development of a
scaling law which would predict the behaviour of the MFPT as the number of
traps increases with their positions defined according to a specific
distribution, in particular, for distributions that globally or locally
minimize MFPT, or other distributions of practical significance. A similar
problem, along with the dilute trap fraction limit of homogenization theory,
was addressed in Refs. [22, 23] for the narrow escape problem involving
boundary traps located on the surface of a sphere in three dimensions.

### Acknowledgements

A.C. is thankful to NSERC of Canada for research support through the Discovery
grant RGPIN-2019-05570.

## References

  * [1] Sidney Redner.  A Guide to First-Passage Processes.  Cambridge University Press, 2001.
  * [2] David Holcman and Zeev Schuss.  The narrow escape problem.  SIAM Review, 56(2):213–257, 2014.
  * [3] PG Saffman and M Delbrück.  Brownian motion in biological membranes.  Proceedings of the National Academy of Sciences, 72(8):3111–3113, 1975.
  * [4] Metzler Ralf, Redner Sidney, and Oshanin Gleb.  First-passage Phenomena and their Applications, volume 35.  World Scientific, 2014.
  * [5] Paul C Bressloff, Berton A Earnshaw, and Michael J Ward.  Diffusion of protein receptors on a cylindrical dendritic membrane with partially absorbing traps.  SIAM Journal on Applied Mathematics, 68(5):1223–1246, 2008.
  * [6] Paul C Bressloff and Jay M Newby.  Stochastic models of intracellular transport.  Reviews of Modern Physics, 85(1):135, 2013.
  * [7] D Holcman and Z Schuss.  Escape through a small opening: receptor trafficking in a synaptic membrane.  Journal of Statistical Physics, 117(5-6):975–1014, 2004.
  * [8] Theodore Kolokolnikov, Michele S Titcombe, and Michael J Ward.  Optimizing the fundamental neumann eigenvalue for the laplacian in a domain with small traps.  European Journal of Applied Mathematics, 16(2):161–200, 2005.
  * [9] Sarafa A Iyaniwura, Tony Wong, Colin B Macdonald, and Micheal J Ward.  Optimization of the mean first passage time in near-disk and elliptical domains in 2-D with small absorbing traps (in press).  arXiv preprint arXiv:2006.12722, 2020.
  * [10] Alexei F Cheviakov and Michael J Ward.  Optimizing the principal eigenvalue of the Laplacian in a sphere with interior traps.  Mathematical and Computer Modelling, 53(7-8):1394–1409, 2011.
  * [11] Jason Gilbert and Alexei Cheviakov.  Globally optimal volume-trap arrangements for the narrow-capture problem inside a unit sphere.  Physical Review E, 99(1):012109, 2019.
  * [12] Wesley JM Ridgway and Alexei F Cheviakov.  Locally and globally optimal configurations of N particles on the sphere with applications in the narrow escape and narrow capture problems.  Physical Review E, 100(4):042413, 2019.
  * [13] A Ismael F Vaz and Luís N Vicente.  A particle swarm pattern search method for bound constrained global optimization.  Journal of Global Optimization, 39(2):197–219, 2007.
  * [14] Jonathan Currie, David I Wilson, Nick Sahinidis, and Jose Pinto.  Opti: Lowering the barrier between open source optimizers and the industrial matlab user.  Foundations of computer-aided process operations, 24:32, 2012.
  * [15] Jeffrey C Lagarias, James A Reeds, Margaret H Wright, and Paul E Wright.  Convergence properties of the nelder–mead simplex method in low dimensions.  SIAM Journal on optimization, 9(1):112–147, 1998.
  * [16] Der-Tsai Lee and Bruce J Schachter.  Two algorithms for constructing a delaunay triangulation.  International Journal of Computer & Information Sciences, 9(3):219–242, 1980.
  * [17] Sarafa Iyaniwura, Tony Wong, Michael J Ward, and Colin B Macdonald.  Simulation and optimization of mean first passage time problems in 2-d using numerical embedded methods and perturbation theory.  arXiv preprint arXiv:1911.07842, 2019.
  * [18] Neil JA Sloane, Ronald H Hardin, TDS Duff, and John H Conway.  Minimal-energy clusters of hard spheres.  Discrete & Computational Geometry, 14(3):237–259, 1995.
  * [19] Vinothan N Manoharan, Mark T Elsesser, and David J Pine.  Dense packing and symmetry in small clusters of microspheres.  Science, 301(5632):483–487, 2003.
  * [20] M Saint Jean, C Even, and C Guthmann.  Macroscopic 2D Wigner islands.  EPL (Europhysics Letters), 55(1):45, 2001.
  * [21] MR Hoare and P Pal.  Physical cluster mechanics: Statics and energy surfaces for monatomic systems.  Advances in Physics, 20(84):161–196, 1971.
  * [22] Alexei F Cheviakov, Michael J Ward, and Ronny Straube.  An asymptotic analysis of the mean first passage time for narrow escape problems: Part ii: The sphere.  Multiscale Modeling & Simulation, 8(3):836–870, 2010.
  * [23] Alexei F Cheviakov and Daniel Zawada.  Narrow-escape problem for the unit sphere: Homogenization limit, optimal arrangements of large numbers of traps, and the n 2 conjecture.  Physical Review E, 87(4):042118, 2013.