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 nondetecting. An example is the slow diffusive motion as the detecting mode and fast, directed ballistic motion as the nondetecting 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 firstpassage 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 reactionescape 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 firstpassage 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 firstpassage 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 firstpassage time density function as a function of the initial conditions, but it
is much easier to solve the timeindependent 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 nontunable parameters of the stochastic firstpassage 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 nontunable 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 nontunable 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 nonisotropic) 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 reactionescape 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 Tcells 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 (oneparameter) family is specially designed for solving the narrow escape problem efficiently and only investigated in that scenario. The second (twoparameter) 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 reactionescape problem for two particles, i.e an intermittent searching predatorparticle is looking
for a mobile preyparticle. After having found the prey, the particlecomplex 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 FokkerPlanck 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 (BallisticBallistic) and BD (BallisticDiffusive), 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) 
Nondimensionalisation
In order to reduce the number of parameters to a minimal independent set, characteristic length and timescales 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 timeindependent 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 gaussianlike 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.
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 .
radialperipheral 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.
Iii The algorithm
Due to the integro type of Eq. (3) and/or the large number of spatial coordinates in twoparticle problems, it is not possible to solve the complete FokkerPlanck 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 firstpassage kinetic Monte Carlo methods oppelstrup2006 (); oppelstrup2009 (); donev2010 ()
are currently the most powerful tools for simulating diluted reactiondiffusion 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 firstpassage 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 noninteracting particles in a sphere of radius
Algorithmically, the case of several noninteracting 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.
Based on the solution of the diffusive initial value problem in (appendix Eq. (36)) it is possible to generate stochastically a firstpassage time to the boundary of , where is sampled according to the corresponding firstpassage 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.
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 diffusionannihilation 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 onetoone relation between time and space in the case of ballistic motion.
For a better understanding, the pseudocode details are shown in Algorithm 1, exemplarily for the BD boundary condition.
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 pseudocode 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 masterequation
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.
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.
If the initial position of the particle is equally distributed within the sphere, and exactly
decrease by for all .
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 .
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.
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 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.
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.
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.
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 .
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 ”breakeven“ 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 xaxis 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 .