Numerical analysis of homogeneous and inhomogeneous intermittent search strategies

Numerical analysis of homogeneous and inhomogeneous intermittent search strategies

Abstract

A random search is a stochastic process representing the random motion of a particle (denoted as the searcher) that is terminated when it reaches (detects) a target particle or area the first time. In intermittent search the random motion alternates between two or more motility modes, one of which is non-detecting. An example is the slow diffusive motion as the detecting mode and fast, directed ballistic motion as the non-detecting mode, which can lead to much faster detection than a purely diffusive search. The transition rate between the diffusive and the ballistic mode (and back) together with the probability distribution of directions for the ballistic motion defines a search strategy. If these transition rates and/or probability distributions depend on the spatial coordinates within the search domain it is a spatially inhomogeneous search strategy, if both are constant, it is a homogeneous one. Here we study the efficiency, measured in terms of the mean first-passage time, of spatially homogeneous and inhomogeneous search strategies for three paradigmatic search problems: 1) the narrow escape problem, where the searcher has to find a small area on the boundary of the search domain, 2) reaction kinetics, which involves the detection of an immobile target in the interior of a search domain, and 3) the reaction-escape problem, where the searcher first needs to find a diffusive target before it can escape through a narrow region on the boundary. Using families of spatially inhomogeneous search strategies, partially motivated by the spatial organization of the cytoskeleton in living cells with a centrosome, we show that they can be made almost always more efficient than homogeneous strategies.

pacs:

I Introduction

The successful usage of efficient search strategies is one of the most important needs in biology and human behavior. It can be observed on all length scales of life and in all kinds of complexity. Just to mention a few examples, humans use them for pattern recognition Najemnik2008 (). Predators apply certain strategies for hunting their moving prey Humphries2010 (). Ants use special techniques to find each other after being separated while being on a tandem run Franks2010 (). Some eukaryotic cells improve their chance to find a target by performing random walks with characteristic persistent time and persistent lengths, even in the absence of external signals Li2008 (). And there are many more observed examples in biological literature.
Although all these examples are quite different and seem to have nothing in common, they can commonly be described by first-passage processes Redner2001book (), which are stochastic processes that end if a certain event happens for the first time . The probability density for the time contains all temporal information about the efficiency of the search strategy. In Mattos2012 () it is shown, that one sometimes has to be careful with the reduction of this information to only one value, the so called mean first-passage time (MFPT)

Nevertheless, this is in most cases the only property which is used to classify the efficiency of the search strategy. Apart from the obvious reason of simplification for comparison, there is a second reason for this reduction: Often it is very hard or even impossible to calculate the whole first-passage time density function as a function of the initial conditions, but it is much easier to solve the time-independent differential equation system for its first moment, which is derived with the help of corresponding backward equations Risken1996book (); Redner2001book ().
The MFPT is a function of the tunable and the non-tunable parameters of the stochastic first-passage process. Typical tunable parameters are for example the persistence length in random walks Tejedor2012 (), the desorption rate in surface mediated diffusion Benichou2010 (); Calandre2014 () or the resetting rate in random motion with stochastic resetting Evans2011 (); Kusmierz2014 (). Typical non-tunable parameters of the search problem are for example the target size, the detection rate, the size and shape of the searching domain and constants of motion (velocity, diffusivity). A complete set of tunable parameters defines a search strategy for the problem which is defined via the non-tunable parameters. Hence, the best strategy is the set of tunable parameters which minimizes the MFPT .

A frequently used way of modeling real search is a so called intermittent search Benichou2005 (); Benichou2005a (); Benichou2007 (); Benichou2011 (); Loverdo2008 (); Bressloff2009 (); Loverdo2009 (); Smith2001 (); Bressloff2012 (). The searcher switches between phases of fast directed ballistic motion, during which it cannot recognize a target and phases of slow diffusion for detecting a target.

For a given size and shape of the search domain and the target, the efficiency, i.e. the MFPT , of an intermittent search still depends on a number of parameters. Since increasing the diffusion constant for the diffusive mode or increasing the velocity modulus for the ballistic mode always decreases the MFPT, even if done only locally, both are assumed to be fixed in the following. Then the MFPT is a function of the switching rates between both motility modes and a functional of distribution of the directions into which the searcher moves after a switch to the ballistic mode. If the searcher does not have a knowledge about his position in the search domain and the search domain is homogeneous such that at no position in the search domain certain directions for ballistic motion are preferred the directional distribution can be assumed to be uniform over all solid angles - as was done in Benichou2005 (); Benichou2005a (); Benichou2007 (); Benichou2011 (); Loverdo2008 (); Loverdo2009 (). This we denote as a spatially homogeneous (and isotropic) intermittent search strategy.
If on the other hand ballistic motion is only possible along predefined tracks, like in molecular motor assisted intracellular transport along the filaments of the cytoskeleton Alberts2014book (), or in cases the searcher utilizes any other transport network, the directional distribution for the ballistic motion should be described by a spatially inhomogeneous direction distribution, which then must represent the spatial organization of the tracks. Also in cases when the searcher does have knowledge about its position in the search domain and about its shape it might be more efficient to move in certain regions of the search domain preferentially into other directions than in other regions. An intermittent search strategy with a spatially varying direction distribution we denote as a spatially inhomogeneous (and non-isotropic) search strategy. In a recent letter schwarz2016 () we introduced the concept of spatially inhomogeneous intermittent search strategies and we presented results that showed that their optimum is in general more efficient than the optimum of homogeneous search strategies. In this paper we elaborate these and more results in detail, explain the computational techniques and show all computations explicitly.

Thus the goal of this paper is to compare the efficiency of spatially homogeneous and inhomogeneous search strategies in spherical domains by determining, numerically, the optimal parameter for different setups: 1) the narrow escape problem , where a searcher has to find a small region on the boundary of the search area, 2) the reaction kinetics enhancement by ballistic motion, where the searcher has to find a immobile target particle within the search domain, and 3) the reaction-escape problem, which combines 1 and 2 such that a searcher has to find a mobile target particle first before it can escape through a narrow region on the boundary of the search domain. The latter example is motivated by a transport process within T-cells attached to a target cell that it is supposed to kill: vesicles loaded with cytotoxic proteins first have to attach to another vesicle containing receptor proteins before they can dock at the immunological synapse, a small region on the cell membrane in contact with the target cell, and release their content there.

Since determining the optimum of the MFPT as a functional of a space and angle dependent direction distribution is not feasible we confine ourselves to two different families of direction distributions. The first (one-parameter) family is specially designed for solving the narrow escape problem efficiently and only investigated in that scenario. The second (two-parameter) family is inspired by the spatial organization of the cytoskeleton of spherical cells with a centrosome. It will be studied for all the three scenarios.

In order to compare the gain of efficiency for different situations, we introduce the dimensionless time

(1)

which is the MFPT of the intermittent search strategy normalized by the MFPT for the purely diffusive searcher. Hence, for an intermittent searcher is more efficient and for a purely diffusive search is faster on average.

The paper is organized as follows: Section II introduces our model of intermittent search in the general case with space and time dependent transition rates. It explains the meaning of the occurring parameters exemplarily in the context of intracellular transport. In almost all cases, it is not possible to solve the differential equation system of the model in a straight forward way via finite element method (FEM).
In consequence, section III introduces the Green’s function method, which is used to solve the model stochastically.
Section IV faces the classical narrow escape problem, meaning, a particle looks for a certain region at the boundary. For the purely diffusive scenario the scaling of the MFPT as a function of the size and the position of the target area is understood for quite a large range of problems Krapivsky1996a (); Redner2001book (); Benichou2005b (); Singer2006a (); Singer2006b (); Singer2006c (); Schuss2007 (); Schuss2012 (); Benichou2008 (); Chevalier2011 (); Cheviakov2012 (). Even in the absence of analytic or asymptotic expressions, the purely diffusive MFPT problem can be solved fast and easily via FEM calculations. For spatial dimensions this is in most cases not possible for the master equation system of intermittent search (Eqs. (2)-(3)) due to the integro type of the partial differential equation. As far as we know, there are no studies on the intermittent search narrow escape problem in a sphere available. Hence, we start the numeric study of this problem in the case of a homogeneous velocity direction distribution. Afterwards we modify the velocity direction distribution to show, that there are more efficient strategies than a homogeneous one.
Section V asks for the best search strategy for a target located within the sphere. In the case of a homogeneously distributed velocity direction and a target which is centered in the middle of the sphere, there are studies on this problem Benichou2011 (); Loverdo2008 (); Loverdo2009 (). We numerically confirm their results, including the very weak dependence of the MFPT on the transition rate from diffusive to ballistic motion, but disprove their optimality assumption for . Furthermore, we study less homogeneous cases, for which there are no MFPT expressions available up to now.
Section VI finally faces a reaction-escape problem for two particles, i.e an intermittent searching predator-particle is looking for a mobile prey-particle. After having found the prey, the particle-complex has to find a small escape area at the boundary. Again, there are already some results for purely diffusive predators in different domains Redner2001book (); Krapivsky1996 (); Redner1999 (), but not for intermittent searching ones in a spherical domain.
Finally, appendix A introduces exact and very fast methods to sample the later defined probability densities of the algorithm of section III .

Ii The model

Intermittent search is generally based on (at least) two different phases for a searcher Benichou2011 (). On the one hand, there is a searching phase of slow (or none) motion, in which the searcher is able to detect a target. On the other hand, there is a relocation phase of directed fast motion without the ability of target detection. Commonly, and also in our case, the searching phase is modeled by pure diffusion with diffusivity . The probability density for being in the diffusive state at position at time will be called in the following, where denotes the search volume of the particle. The relocation phase is modeled by straight ballistic motion. The probability density for being in the ballistic state at position at time and moving with velocity is denoted in the following, where is the unity vector in direction of the solid angle .

In intracellular transport vesicles (proteins, organelles) switch between diffusion within the cytosol and almost ballistic motion by molecular motor assisted movement along cytoskeleton filaments. The density of these filaments in direction of the solid angle is generally very inhomogeneous in space: for instance in cells with a centrosome microtubules emanate radially from the centrosome towards the cell periphery, where the actin cortex, a thin sheet of actin filaments underneath the cell membrane, provides transport in random directions. Sometimes the filament density even varies over time (for instance during cell polarization). In consequence, the likelihood of a switch between the two phases and the choice of the ballistic direction generally depends on the position of the searcher. Formally we describe a spatially varying distribution of directions by the density . It is proportional to the rate of a switch from diffusive to ballistic motion in direction at position at time . In the context of intracellular transport it can be interpreted as the filament density of the cytoskeleton in direction .

The master equation system of our model for one searching particle is given by the Fokker-Planck equation system:

(2)
(3)

where and are transition rates from diffusive to ballistic motion and vice versa. In the context of modeling intracellular transport they are the attachment and detachment rates (from cytoskeleton filaments).
In consequence, the diffusing searcher experiences a total annihilation rate

(4)

with which it is transformed into a ballistically moving particle with a randomly chosen direction (and velocity ) with probability

(5)

A ballistically moving particle switches back to diffusive motion with rate .

Within this article, a target shall always be detected immediately, when the diffusive searcher reaches the target area for the first time (in reaction kinetics this means reaction upon contact). One could also consider detection or reaction with a finite rate within the target area Benichou2011 (). But we restrict ourselves to the case , i.e. target detection is always modeled via the boundary condition

(6)

where either is the detection area within or is the detection area at the surface of (narrow escape problem).

In section VI, we consider the problem of two moving particle, which will react immediately if their distance becomes smaller than a certain value. Hence their probability distributions are not independent and the solution does not factorize. Consequently, the master equation system depends on spatial coordinates and coordinates for . As the exact notation of this master equation system and the corresponding boundary conditions is very lengthy but straightforward, we will skip it here.

Apart from the initial conditions , the equation system is augmented by boundary conditions at the boundary for and boundary conditions at for . Two different boundary conditions, in the following called BB (Ballistic-Ballistic) and BD (Ballistic-Diffusive), will be used within the studies of this article.

BB boundary condition

For BB boundary conditions we assume that a ballistically moving particle hitting the boundary is simply reflected and stays in the ballistic mode, no matter whether this happens at the target area or not (i.e. the target is not detected then). A particle in the diffusive mode is reflected at every point of the boundary , which does not belong to the target area and stays in the diffusive mode. Fig. 1a visualizes the BB condition in a sketch. Formally this boundary conditions are described by

(7)

where denotes the outward pointing unity vector perpendicular to the boundary at position and denotes the solid angle which belongs to the reflection of at the surface position .

BD boundary condition

For BD boundary conditions we assume that a ballistically moving particle hitting the boundary switches to the diffusive motion. If this part of the boundary belongs to the target area, the particle is immediately detected. A particle in the diffusive mode is reflected at every point of the boundary , which does not belong to the target area and stays in the diffusive mode. Fig. 1b visualizes the BD condition in a sketch. Formally this boundary conditions are described by

(8)
Figure 1: Sketch of boundary conditions in a search volume with escape area (dotted line). Grey wiggly lines represent diffusive motion, green lines ballistic motion. a) BB: a ballistically moving particle is reflected at the boundary and stays in the ballistic mode, i.e is only detected if it is reached diffusively. b) BD: a ballistically moving particle switches to diffusion at the boundary, i.e will be detected if it reaches ballistically (trajectory, starting at ) or diffusively (trajectory, starting at ).

Nondimensionalisation

In order to reduce the number of parameters to a minimal independent set, characteristic length- and time-scales where chosen by introducing the dimensionless spatial and temporal coordinates

(9)

In consequence, Eqs. (2) and (3) are always solved in the unit sphere and look the following way:

(10)
(11)

with

(12)

Apart from the sphere radius , the absolute value of the velocity also vanished in the dimensionless coordinates, as holds. Furthermore, is not changed by the dimensionless units, i.e. .

Models for the direction distribution

Eq. (4) introduced the total transition rate for a switch from diffusive to ballistic motion at position at time . Although it is numerically possible to handle this most general scenario (Algorithm 1 in section III) this rate will be constant in time and space in the investigated models, i.e. without further loss of generality we set and in consequence Eq. (5) simplifies to

(13)

Within the studies of this paper, two different families of time-independent inhomogeneous distributions will be compared to the homogeneous distribution

(14)

Both will be rotational symmetric, i.e. depends only on the radius and the angle

(15)

between the vectors and .
This symmetry also holds for the homogeneous case of , where the probability density for the angle is independent of and given by

(16)

varying Gaussian distribution

The distribution for the angle (Eq. (15)) , introduced now, will only be applied to the narrow escape problem. The principle idea is to find the probability density , which minimizes the MFPT of the narrow escape problem. Mathematically, this is a variational problem. In consequence, a numeric solution requires an apriori assumption for a class of density functions, which is motivated now:
If the particle is close to the center of the simulation sphere, a mainly radially outward pointing velocity direction is for sure the best strategy, as it is the fastest way to reach the sphere’s boundary. At the boundary this distribution is not optimal any more, as there is no velocity component in parallel to the boundary. Without this parallel component, the searcher gets stuck at a relatively small part of the boundary.
In consequence, the spread of the distribution should increase with . Following this argumentation, the gaussian-like probability density , illustrated in Fig. 2, was chosen for our simulations:

(17)

where

(18)

denotes the spreading of the gaussian and

(19)

is the normalization of the distribution.

Figure 2: The class of probability densities (Eq. (17)) as a function of the variable and the spreading parameter .

The class parameter controls the speed of the increase of the distribution spreading. For , the velocity direction points radially outwards for all as tends to zero. The spread (Eq. (18)) increases monotonically in and in . For , we are dealing with the totally homogeneous velocity direction distribution .

radial-peripheral distribution

The second investigated distribution is inspired by the spatial organization of the cytoskeleton of spherical cells with a centrosome and was introduced in schwarz2016 (), see Fig. 3a for a sketch. It contains two parameters:

(20)

The parameter is the probability to move radially outwards, and the probability to move inwards inside the inner spherical region with radius . represents the width of the outer shell in which the homogeneous strategy is applied, hence represents the totally homogeneous searching strategy. A ballistically moving particle switches to the diffusive state when it reaches the radius and . The distribution will only be investigated for the boundary condition BD. A sketch of the resulting stochastic processes is given in Fig. 3b.

Figure 3: a) The cytoskeleton transport network of a spherical cell with a centrosome: Microtubules (green lines), emanating from the MTOC of the cell close to the nucleus, are orientated predominantly radially to the cell membrane. Kinesin and Dynein motor proteins transport cargo along them. The actin cortex (red lines) close to membrane is built by isotropically orientated actin filaments. Myosin motors transport cargo along. b) Sketch of the stochastic process with the direction distribution . is the thickness of the outer region. is the starting point of the particle, Grey wiggly lines represent diffusive motion, green lines ballistic radial motion (for , outward with probability , inward with probability ), red lines ballistic motion in random directions(for ).

Iii The algorithm

Due to the integro type of Eq. (3) and/or the large number of spatial coordinates in two-particle problems, it is not possible to solve the complete Fokker-Planck equation system (Eqs. (2)-(3)) via FEM. Only the purely diffusive case of one particle is always solvable. Hence, the numerical results of this article were mostly derived with Monte Carlo techniques, which will be explained in this section.

Green’s function reaction dynamics Zon2005a (); Zon2005b () and first-passage kinetic Monte Carlo methods oppelstrup2006 (); oppelstrup2009 (); donev2010 () are currently the most powerful tools for simulating diluted reaction-diffusion processes. In contrast to the traditional way of simulating diffusion by an enormous number of very small (compared to the system size) random hops, they propagate diffusing particles randomly within so called protective domains over rather long distances. The core of these methods are Green’s functions, the solution of the initial value diffusion problem within the protective domains. In essence, these methods work the following way:
For a given starting configuration of interacting diffusing particles within a domain , a protective domain is assigned to each particle with for . A necessary restriction for the choice of each domain is the knowledge of an analytic expression for the Green’s function for the initial value diffusion problem according to absorbing boundary conditions at the interior of and the boundary conditions of at common boundaries of and (as far as they exist). Based on these Green’s functions it is possible to sample for the particle which will leave its domain first and a corresponding time for this first-passage event. Finally, the exit position is sampled depending on . If the distance of to the protective domains of all other particles is larger than a given threshold, we look for a new protective domain for the particle . Otherwise, we have to sample new positions for all particles, whose protection domains are too close to . In the end, a new protective domain has to be assigned to all these particles.
In schwarz2013 () we developed an improvement of these routines for a wider range of applications including external space and time depending transition rates.
For a more detailed general explanation of these methods and for proofs of their correctness, the reader is referred to the original articles Zon2005a (); Zon2005b (); oppelstrup2006 (); oppelstrup2009 (); donev2010 (); schwarz2013 (). The rest of this methodical chapter only focuses on describing the concrete algorithm for particles in a sphere, switching between ballistic and diffusive motion according to the model definition of section II. The method will be explained in the most general context of spatially and temporally varying rates, see Eqs. (2 -3).

iii.1 non-interacting particles in a sphere of radius

Algorithmically, the case of several non-interacting particles is identical to the case of only one particle. Consequently, we restrict the following algorithm description to only one particle.
For a diffusive particle being at position at time we use two different types of domains for propagating the particle within the simulation sphere of radius . If the distance to the boundary of the sphere is larger than a very small threshold value , a sphere with radius , centered around , will be assigned to the particle. An example for such a situation is the blue particle in Fig. 4.

Figure 4: Illustration for the choices of protecting spheres/cones: The blue particle on the left is closer to the boundary of the simulation sphere than to any other particle, in consequence the radius of its protecting sphere is limited by the distance to the boundary. The green and the red particle in the middle are closer to each other than to the cell boundary. If they can react, their radius is limited by their distance. The gray particle has reached the boundary. Hence it is propagated in a cone, which is illustrated in the 2d projection.

Based on the solution of the diffusive initial value problem in (appendix Eq. (36)) it is possible to generate stochastically a first-passage time to the boundary of , where is sampled according to the corresponding first-passage time probability (appendix Eq. (A.2)). If the particle does not switch to ballistic movement before time , a random position update to the boundary of the sphere is done and a new protective domain must be assigned to the particle afterwards. Otherwise a new particle position within the sphere is sampled by using the radial probability density (appendix Eq. (A.2))

Due to the fact, that the boundaries of the sphere and the protecting sphere have always only one point in common, it does not work to use only spheres for the protecting domains. With probability one, the particle will touch the boundary of the simulation sphere, i.e. we would end up with an infinite sequence of protecting spheres, whose radii converge to zero. The best possibility to overcome this problem would be the usage of protection domains, whose boundaries coincide locally with the boundary of the simulation sphere in an area and not just in one point. Due to the missing knowledge of corresponding Green’s functions and/or the ability to sample efficiently within these domains, this is not possible. In consequence, for , we locally approximate the boundary of the simulation sphere by a suitable geometry, which is a spherical cone with a reflecting conical and an absorbing spherical boundary (appendix Eq. (39)). An example for such a situation is the gray particle in Fig. 4. If the distance to the boundary is larger than after being propagated within the cone, we again go on with a protecting sphere, otherwise, we use again a cone. The accuracy of this method is tunable via the two geometry boundary approximation parameters and the maximum radius of the protecting cone. It is important to mention, that is just an upper limit for the cone’s radius. If the center position of the cone is closer than to the target area , is chosen to be the minimal distance of to .

In order to demonstrate the high accuracy, we compared our Monte Carlo method with the solution of a commercial FEM solver for a purely diffusive narrow escape problem. Starting at , the searcher has the find the escape area with polar angle (bright area at the top in Fig. 4). The FEM simulation was done on a very fine triangulation ( 200000 elements) using the rotation symmetry of the problem and yields the expectation value for the needed search time. The Monte Carlo simulation with samples was done for the geometry approximating parameters , and yields the almost perfectly matching value . A much stronger criterion than the comparison of expectation values is the equality of the survival probability (probability of not having reached the escape area until ) for all . The again almost perfectly matching result is shown in Fig. 5.

Figure 5: Comparison of the Monte Carlo method and the solution of a FEM solver on the basis of the survival probability for a purely diffusive narrow escape process with and . The inset shows the very small relative deviation of these distributions.

All numerical results of this article are expected to be in the same numerical exactness (expect the sampling deviation in the case of a smaller number of samples), as the values of and where always chosen to be on the save side according to the smallest occurring length scale. However, the accuracy was successfully checked wherever this was possible (either by analytic values or FEM values).

For a diffusive particle at position , the total rate for a switch to a ballistic movement in an arbitrary direction is given by (see Eq. (4)). If this rate is spatially inhomogeneous, the methods of Zon2005a (); Zon2005b (); oppelstrup2006 (); oppelstrup2009 (); donev2010 () will fail, as there is in general no analytic solution (Green’s function) to the diffusion-annihilation equation available. The algorithm, presented in schwarz2013 (), overcomes this problem by using a spatially maximal rate

(21)

in order to sample a candidate time for a switch from diffusive to ballistic motion according to the probability density

(22)

A new position is assigned to the particle with the help of (appendix Eq. (A.2)). With probability the particle moves on diffusively. With probability it switches to ballistic motion with velocity , sampled according to the probability density (see Eq. (5)).
For a back switch to diffusive motion only a corresponding time must be sampled, as there is a one-to-one relation between time and space in the case of ballistic motion.
For a better understanding, the pseudo-code details are shown in Algorithm 1, exemplarily for the BD boundary condition.

1:Input:
2:Output:
3:;
4:;
5:;
6:repeat
7:     if ( then
8:         if (then
9:               random number according to ;
10:         end if
11:         Choose the protecting sphere/cone with            maximal radius as a function of ;
12:          + random number according to            for ;
13:         if (then
14:               rand. position at absorbing part of ;
15:              ;
16:         else
17:               rand. position update within according                   to ;
18:              ;
19:              if  then
20:                  ;
21:                   random velocity according to ;
22:              end if
23:         end if
24:     else
25:          time when ballistic particle hits boundary;                  
26:          random exponentially distributed number                   with rate ;
27:         if (then
28:              ;
29:              ;
30:         else
31:              ;
32:              ;
33:         end if
34:         ;
35:     end if
36:until  (distance to absorbing part of sphere threshold)
37:return
Algorithm 1 one particle

iii.2 interacting particle in a sphere of Radius

If there are at least two particles in the simulation sphere, which are able to react, the choice of the protection boxes does not only depend on the position of the particle, but also on the distance between these reacting particles. In general, protecting spheres/cones of reacting particles are not allowed to be closer to each other than the interaction distance . An example for such a situation is the red and green particle in the middle of Fig. 4. A similar problem as the boundary approximating problem in the subsection before has to be solved here. If we choose the protecting spheres/cones of interacting particles always in a way, that the boundaries of these protection boxes have their minimal distance in only one point, we will for sure end up in an infinite sequence of protection spheres/cones, whose radii tend to zero. In general there are two ways to overcome this problem. The first one is discussed in Zon2005b () and the problem is solved via a coordinate transformation for the two particle positions to the difference vector and the mass point vector. As the problem factorizes in these coordinates, one ends up with two independent problems. Although a position update takes more time in these situations due to the fact, that radial symmetry is lost within the protection boxes in these coordinates, this is a very powerful tool for particles, which are far away (compared to their distance) from the boundary of the simulation sphere. But for particles, whose distance to the boundary is only a little bit larger than their distance to each other, this method does not work well. Hence, we decided to use a second tunable approximation by defining a parameter : If the distance between two reactive particles is less than these particles react. If we choose , we are numerically for sure on the safe side, as all results look totally the same as in the case . A comparison to the solution of a FEM solver is not possible anymore, even for purely diffusive particles, due to the high spatial dimension () of the problem. A pseudo-code description would be quite large and the general idea is the same as in Algorithm 1. The interested reader is again referred to Zon2005a (); Zon2005b (); oppelstrup2006 (); oppelstrup2009 (); donev2010 (); schwarz2013 ().

Iv Narrow escape problem

The narrow escape problem for a purely diffusive particle in a sphere (and other simple domains) has already been studied in several publications. A nice overview, containing analytic asymptotic expressions, is given in Cheviakov2012 () and Schuss2012 (). Within this section, we consider the problem of a particle, which moves according to an intermittent search strategy, meaning, the master-equation system of its movement is given by the Eqs. (10) and (11) until it reaches the absorbing part of the boundary of the simulation sphere for the first time. This escape area is given by a spherical cab with polar angle , like it is shown at the north pole of Fig. 4. The position of this cap is of course not known by the particle.
The MFPT to the absorbing cap is a function of and the velocity direction distribution . Furthermore, it depends on the initial position of the diffusively starting particle. But in the case of small target areas, the relative influence of the initial position totally vanishes. Depending on the diffusivity and , we study the optimal solution to the escape problem, i.e. we always look for the transition rates which minimize (and simultaneously ).

purely diffusive search

The reference time for a purely diffusive searcher ( and/or ) is inversely proportional to the diffusivity . Among others, Cheviakov2012 () has derived a very exact analytic approximation of the problem for small for arbitrary starting positions . For ,

(23)

holds. has been calculated via - Monte Carlo samples for each and compared to the analytic approximation in Eq. (23). The result is shown in Fig. 6.

Figure 6: and as a function of for : Each red dot is the average of Monte Carlo samples. It coincides very well with the analytic approximation of Cheviakov2012 () (green line), given in Eq. (23). The inset shows the relative difference between the curves.

For small values of the relative deviation between the simulated value of and is extremely small and only based on stochastic fluctuations (inset of Fig. 6 ). For larger values of it slightly increases, which is not based on a drop of exactness in our numerical routines, but on the fact that the approximation becomes worse for larger opening angles. If the initial position of the particle is equally distributed within the sphere, and exactly decrease by for all , which has also been checked numerically.

random velocity model

Before studying intermittent strategies, it is insightful to have a look at the opposite choice of transitions rates and , which is a random velocity model, given by the limit and .
For the BB boundary condition the corresponding MFPT tends trivially to infinity for all , as the ballistically moving particle is reflected at the boundary without target area detection and never switches to diffusive mode (see Fig. 1b).

For the BD boundary condition, this is not the case. The resulting random velocity model is given by a ballistically moving particle, which detects the escape area at the boundary when reaching it and randomly chooses a new direction for the ballistic motion when reaching a part of the sphere’s boundary which does not belong to the target area. In consequence, the corresponding MFPT depends on the velocity direction distribution , the opening angle and slightly on the initial position . For and the case of a homogeneous velocity direction density (Eq. 14) we derived an approximating expression for :

(24)

Fig. 7 shows ( samples) and in a logscale plot. The relative deviation of and vanishes for , which is shown in the inset.

Figure 7: and as a function of for : Each red dot is the average of Monte Carlo samples. It coincides very well with the analytic approximation (green line), given in Eq. (24). The inset shows the relative difference between the curves.

If the initial position of the particle is equally distributed within the sphere, and exactly decrease by for all .

Figure 8: narrow escape, BB, ; The normalized MFPT (Eq.1) is color-coded as a function of the rates and for different values of and (interpolation from a grid of data points each). samples have been done for each pair ( , ), which leads to relative stochastic fluctuations of , which are smaller than 0.2%. The position of the minimum is always shown with a red dot.

A comparison of and points out an important difference in the behavior of divergence of and for small escape areas:

After having studied the two possible extreme cases in search behavior, which is necessary for understanding the later discussed dependence, we now face intermittent strategies and analyze their efficiency.

In subsection IV.1 the condition BB is studied, i.e. a ballistically moving particle, which hits the boundary of the simulation sphere, stays in its ballistic mode with the reflected velocity direction. The arrival at the escape area of the sphere will only be detected if the particle is in the diffusive mode, otherwise it is reflected. We compare the problem of the homogeneously distributed direction density to the inhomogeneous scenario of .
In subsection IV.2 the condition BD is studied, i.e. a ballistically moving particle, which hits the boundary of the simulation sphere, switches immediately to the diffusive mode, i.e. if this switch happens at the escape area, the particle immediately recognizes the exit. Here, the homogeneous case is compared to the inhomogeneous scenarios of and .

iv.1 Bb

For the BB condition, the searcher will start in the center of the sphere and the escape area is given by a spherical cab with angle within this subsection, i.e. the radius of the absorbing spherical cab is seven times smaller then the radius of the sphere. In consequence, the area the particle searches, is about of the total spherical surface, i.e. we are in the limit of a small escape area. In this setup, the reference time (taken from the MC data of Fig. 6) is given by

(25)

homogeneous distribution

For different values of , we look for the best strategy to search for the absorbing area as a function of the switching parameters and . as a function of and is shown in Fig. 8 for four different examples of . For larger than about 0.025 there is no benefit in a mixed strategy. Here, a purely diffusive particle is on average the better searcher as diffusive motion is faster on these scales. As decreases, phases of ballistic displacement become more and more efficient, as the diffusive displacement per time unit shrinks. Hence, a global minimum occurs in the () space, i.e. there is a benefit in an intermittent search strategy. As expected, this benefit further increases with decreasing , i.e. increases monotonically and

(26)

holds, although . Surprisingly, the efficiency of the strategy changes only very little in a quite large (relative to the absolute values) surrounding of the optimal solution ) for all diffusivities . This can be seen by having a closer look to the values of the isolines in Fig. 8. In consequence, due to stochastic fluctuations, the relative error in the optimal values of and is much larger than the relative error in the value of and . Fig. 9 shows , and (in the inset) as a function of the diffusivity .

Figure 9: narrow escape, BB, ; The optimal transition rates , and the resulting normalized MFPT (inset) as a function of . For the diffusion coefficients shown in Fig. 8, the data coincides with the coordinates of the red dots and its value of .

The corresponding values of and are also listed in Table 1 and plotted in Fig. 11 for a comparison to the later treated inhomogeneous search scenarios. and decrease monotonically in . As this happens faster for than for , the fraction of time spend in the diffusive mode increases with .

Due to the enormous numerical effort, it is not possible to vary systematically here. Nevertheless, we exemplarily investigated also some values of for smaller and larger values of . Similar to the results in the following chapters, we found, that a decrease in target size results in an increase in both transition rates.

inhomogeneous distribution

Within this subsubsection, we study the influence of on the search strategy and the transition rates. Fig. 10 shows the dependence of on for representatively selected values of and the corresponding optimal parameters of the homogeneous scenario, shown in Fig. 9.

The global minimum for each will be denoted in the following.

Figure 10: narrow escape, BB, ; as a function of the spreading parameter for the transition rates of the optimal solution of Fig. 9. The minimum position (gray bar) of all curves is almost identical. The colored numbers on the right show the ratio .

At the values of coincide with the corresponding of Fig. 9. In none of the cases, the value at is the minimum. It follows, that an anisotropic velocity direction distribution increases the efficiency of the search strategy significantly. For small values of , is much smaller than , which can also be seen by comparing the blue and the green curves of Fig. 11 and the corresponding values in Table 1. As increases the benefit of an inhomogeneous strategy becomes less pronounced. It is remarkable, that the degree of inhomogeneity is constant for all .

Nevertheless it is even possible to decrease further: In Fig. 10 the transition rates were chosen as the optimal solution for the homogeneous case. There is no reason, that this is also the optimal choice in the inhomogeneous case. In consequence, we varied and simultaneously for finding the optimal parameters and for the MFPT (be aware of the different meaning of the index ”opt“ and ”OPT“). The results are shown in Table 1 and is plotted in Fig. 11.

0.02 386 371 307 238 0.8 3.4 11.5 8 0.35 0.325
0.015 514 465 337 264 1.5 4.1 18 8.5 0.35 0.325
0.01 771 610 349 297 4 6 25 9.5 0.35 0.325
0.0075 1028 720 377 321 5 6.3 30 10 0.35 0.325
0.005 1542 888 398 353 8.4 7.9 36 11 0.35 0.325
1/300 2314 1071 433 386 12 9.6 42 12 0.35 0.325
3085 1211 448 410 15 10.75 48 13 0.35 0.325
0.002 3856 1326 466 429 17.5 11.5 50 13.5 0.35 0.325
0.001 7712 1727 530 492 26 15 60 15 0.35 0.3
15420 2217 618 562 38 20 75 18 0.35 0.3
38560 3026 740 670 55 25.75 95 21 0.35 0.3
Table 1: : purely diffusive MFPT; : optimized intermittent MFPT for with optimal rates and ; : optimized intermittent MFPT for with optimal inhomogeneity coefficient and fixed rates and ; : optimized intermittent MFPT for with optimal inhomogeneity coefficient and corresponding optimal rates and .

Table 1 delivers some remarkable results:

  • The optimal value of seems to be almost constant in all cases. For the rates of the homogeneous optimization and for the rates of the inhomogeneous optimization the best solution is always given by . Hence, the degree of inhomoegeneity for an optimal solution does not seem to depend much on the diffusion coefficient and the transition rates, which is quite surprising.

  • Comparing the values of with , one recognizes remarkable changes in the transition rates. For large values of , the change is more than a factor of 10.

  • Like in the homogeneous case, the efficiency of the inhomogeneous strategy changes only very little in a quite large surrounding of the optimal solution .

For concluding this subsection, Fig. 11 shows the optimal MFPTs for the different discussed scenarios.

Figure 11: narrow escape, BB; MFPT for purely diffusive search (Eq. 25, red line); optimal search with homogeneously distributed velocity direction (green); optimal search with inhomogeneously distributed velocity direction for the fixed rates , for the homogeneous scenario (blue) ; optimal search with inhomogeneously distributed velocity direction for rates , (black).

Compared to a purely diffusive searcher (red), an intermittent search strategy with a homogeneous velocity direction distribution (green) optimizes the search process especially for small significantly, which has already been shown in the inset of Fig. 9. In the next step, we introduced an inhomogeneity in the velocity direction distribution (blue), but kept the optimal rates of the homogeneous case. Again, the largest benefit can be seen for small (see Fig. 10). In the last step, we varied the transition rates and the degree of inhomogeneity simultaneously (black). Although the optimal rates changed dramatically, the additional benefit is much smaller than in the optimization steps before. But this time it increases with .

iv.2 Bd

For all investigated direction distributions ( and both inhomogeneous scenarios , ) in this subsection the optimal search strategy is either a purely diffusive one (for large) or the simulations yield . Exemplarily, this is shown in Fig. 12. For and different values of the figure shows as a function of the transition rates for the case of and the initial position in the origin.

Figure 12: narrow escape, BD, ; The normalized MFPT (Eq.1) is color-coded as a function of the parameters and for different values of and (interpolation from a grid of data points each). Between and samples have been done for each pair ( , ), which leads to relative stochastic fluctuations of which are smaller than 0.2% for. The position of the minimum is always shown with a red dot.

A comparison to Fig. 8 shows the different behaviour of the optimal solution for the two boundary conditions. We verified also for smaller values of and larger ones (). Consequently, the numerical effort of finding the best strategy is dramatically reduced, as there is one parameter less to vary. Due to this reduced effort, the variation of the absorbing angle will also be studied in the case of a .
Apart from this additional study, the beginning of the subsection is organized identically to the one before: We start with the case of , followed by the inhomogeneous scenario for . In both cases the initial position is the origin. Afterwards we study the case for a homogeneously distributed initial position and .

homogeneous distribution

At first, we study whether an intermittent search strategy or a purely diffusive strategy is better for a given pair of parameters (). For the reason of completeness we faced this question for all values of and not only for a ”narrow” escape area. The result for is shown in Fig. 13.

Figure 13: narrow escape, BD, ; Diagram for the choice of the best search strategy as a function of and : In the red (a) and the green (b) domain, an intermittent search strategy is preferable, whereas in the white domain (c) pure diffusion is the best strategy. For the construction of the diagram, the behaviour at the position of the dots was investigated.
Figure 14: narrow escape, BD, ; Examples of the function ) for the three different colored areas in Fig. 13.

In the red (a) domain an intermittent search strategy is preferable. starts monotonically decreasing at . It follows the global optimum at . An example for this behavior for , is given in Fig. 14.
In the green domain (b), intermittent search is also preferable. Although starts monotonically increasing, it decreases to for some values of . Again, an example for this behavior for , is given in Fig. 14.
Finally, in the white domain (c) holds, hence a diffusive search is the best strategy. An example for this behavior for , is also given in Fig. 14.

Fig. 13 only answers the question about the best strategy in principle, it is neither quantifying the transition rate nor the MFPTs and . A quantification has only been done in the case of small escape areas () due to the following reasons:
If the escape area is large, the searcher will find it soon, hence there is no need for a special strategy. The largest impact of on the efficiency of the strategy is given for small values of , i.e. for large either a purely diffusive searcher or a random velocity model () is always close to the optimal strategy. Additionally, for small values of the optimal strategy is almost independent of the starting position of the searcher, i.e. the shown results for a searcher starting at the origin will remain true in the more general context of an arbitrary initial position. For the angle this independence is shown explicitly later.

Fig. 15 quantifies the values of , and for .

Figure 15: narrow escape, BD, ; , and the corresponding as a function of and (interpolation from the non equidistant grid shown at the bottom of Fig. 13). The thick black line in each subfigure separates the area of intermittent search and purely diffusive search. It coincides with the boundary line between the area a and c in Fig. 13. top: in a logscale color plot, the dashed lines with transparent label show isolines of the purely diffusive search scenario (Eq. (23)) for the reason of comparison. middle: is color-coded plotted bottom: in a logscale color plot.

The corresponding curves, from which the optimal values of and were taken for each data point, qualitatively all look like the red curve in Fig. 14. Depending on and , up to samples have been performed for each parameter triple (, , ). As the depth of the minimum position is differently strong pronounced this is necessary to control the stochastic fluctuations in the value of .
For the optimal strategy is trivially given by for all , i.e. the optimal strategy is the random velocity model with MFPT , which is very well approximated by in Eq. (24), shown in Fig. 7. For small diffusivities the transition rate is finite. Its value strongly depends on the size of the escape area, i.e. on the value of . The thick black line in Fig. 15 shows the ”break-even“ diffusivity , where the optimal strategy changes from intermittent search to purely diffusive search. increases monotonically in . It rises the interesting question about the limit of for (be aware of in Fig. 15). If held, for every there would be a threshold value below which pure diffusion would be the best strategy. In the opposite case of a positive limit , i.e. , intermittent search would be more efficient for all , no matter how small becomes. Due to the divergence of the MFPT for it is not possible to face this limit numerically for the reason of running time. Nevertheless, there are clear arguments for a limit : The second derivatives of the isolines of in the second subfigure of Fig. 15 seems to vanish for small . Hence, they were expected to reach the x-axis in a straight line at positions larger than zero. Due to the enormous running time for very small angles we verified this hypothesis of affine extrapolation only partially at for some .

For comparison to the boundary condition BB and the later investigated inhomogeneous search scenarios, the angle is shown separately in Fig. 16 and the corresponding values of are shown in Fig. 18 and Table 2.

Figure 16: narrow escape, BD, ; and (inset) as a function of for : The corresponding curves from which the minima are taken qualitatively all belong to case (a) in the diagram. For samples for each investigated the position of the minimum