# Mean First Passage Times for Transport Equations

Thomas Hillen Department of Mathematics, University of Alberta, Canada.
<EMAIL_ADDRESS>Maria R. D’Orsogna Department of Mathematics, California
State University at Northridge, Los Angeles, CA, 91330, USA<EMAIL_ADDRESS>Jacob C. Mantooth Deptartment of Applied and Computational Mathematics and
Statistics, University of Notre Dame, Notre Dame, IN, 46656, USA.
<EMAIL_ADDRESS>Alan E. Lindsay Department of Applied and Computational
Mathematics and Statistics, University of Notre Dame, Notre Dame, IN, 46656,
USA<EMAIL_ADDRESS>
###### Abstract

Many transport processes in ecology, physics and biochemistry can be described
by the average time to first find a site or exit a region, starting from an
initial position. Typical mathematical treatments are based on formulations
that allow for various diffusive forms and geometries but where only initial
and final positions are taken into account. Here, we develop a general theory
for the mean first passage time (MFPT) for velocity jump processes. For random
walkers, both position and velocity are tracked and the resulting Fokker-
Planck equation takes the form of a kinetic transport equation. Starting from
the forward and backward formulations we derive a general elliptic integro-PDE
for the MFPT of a random walker starting at a given location with a given
velocity. We focus on two scenarios that are relevant to biological modelling;
the diffusive case and the anisotropic case. For the anisotropic case we also
perform a parabolic scaling, leading to a well known anisotropic MFPT
equation. To illustrate the results we consider a two-dimensional circular
domain under radial symmetry, where the MFPT equations can be solved
explicitly. Furthermore, we consider the MFPT of a random walker in an
ecological habitat that is perturbed by linear features, such as wolf movement
in a forest habitat that is crossed by seismic lines.

###### keywords:

Random walks, Mean First Passage Time, Anisotropic diffusion, Homogenization.

35M13, 92B99, 47G20, 60G40, 60G50.

##  1 Introduction

A key issue in the study of species movement is determining the time necessary
for individuals to find a target, such as a food source, shelter, or mate, for
the first time [33, 60, 52]. In ecology, the time required for an animal to
first reach (or exit from) a region can be recorded via satellite tracking and
the data used to understand how the environment, road layout, seasonality,
climate and the presence of other species affect habitat selection, foraging
and migration patterns. This information may be useful in ecosystem management
[19, 41, 37, 36]. Within cellular biology, identifying the time to first reach
a specific site is also of primary interest, since first arrivals may trigger
irreversible on-site biochemical transformations that indicate the onset of
disease, the initiation of repair, or the completion of extracellular
signaling [20, 5, 55, 35, 38]. Examples include the first time a T-cell finds
an antigen presenting cell during immune response [43], the first time a
fibroblast reaches a wound to initiate healing [3], the first time a base
excision repair enzyme reaches a lesion on a DNA strand via charge transport
[21]. The first arrival time to a target is often a proxy for the efficiency
of the search and a common goal is to expedite, hinder or otherwise control
the search process itself [1, 14, 6, 39]. For example, in receptor-mediated
viral entry a virus can enter a cell only after binding to a critical number
of surface receptors as in the case of HIV, the SARS coronavirus, and
hepatitis C [23, 34, 6]. The timing of gene expression also depends on the
first time a specific set of intracellular proteins reach a threshold
concentration [57, 22]. Related questions arise in chemistry, material and
polymer science, and pertain to identifying the first assembly time of a
cluster starting from individual subunits, determining how fast protein,
filaments or other molecular aggregates reach a predetermined size, how fast
bacteriophage and viral capsids are assembled [64, 65, 18].

Figure 1: Schematic of anisotropic transport in a bounded planar region
$\Omega\in\mathbb{R}^{2}$ in the presence of one-dimensional substructures.
Diffusion is one-dimensional along the linear features and two-dimensional
away from them. This schematic can represent, for example, organelle transport
within a two-dimensional cell, grown on a flat surface, whereby particles
transition from performing linear random walks along microtubules and/or other
filaments emanating from the nucleus, to undergoing two-dimensional planar
diffusion in the cytoplasm. Similarly, it can represent animal motion in a
given environment with roads or other one-dimensional features that direct the
animal’s movement. Due to the stochastic nature of the transport process, each
individual trajectory starting at a given position $x\in\Omega$ and with a
given velocity $v$ will reach the target area, shown in red, at a specific
time. Our goal is to determine the mean first arrival time to the target. We
can model this transport scenario using a kinetic transport equation with
turning kernel given by the bimodal von-Mises distribution in a circular
geometry and assuming that motion is biased along the radial direction.

First arrival times distributions can be quite broad, especially if outliers
or rare events are possible. For some state space trajectories finding the
target may occur quickly, for others the search may take longer, and in some
cases it may never be completed, being halted by degradation or adsorption
into trapped states. A useful quantity is the mean first passage time (MFPT),
the average of all first arrival times of the underlying stochastic process
[56]. Note that since the MFPT is an average over all completion times, it is
enough for one trajectory to never reach the target site or threshold for it
to diverge. In transport phenomena the MFPT depends on dimensionality,
geometry, the type of motion involved, heterogeneity of the environment; in
threshold phenomena it depends on the reversibility of events that drive the
process, cooperativity, and ordering effects [16]. Multiple theories have been
developed to evaluate the MFPT in a variety of spatio-temporal scales and
geometries, mostly (but not all) in the context of diffusion type models [54,
8, 31, 40, 33, 24, 9]. The most common forms of movement are Brownian motion,
random walks, Levy flights, ballistic motion, and their combinations to
represent run-and-tumble motility in bacteria for example, or to model
alternating periods of movement and rest [4, 48, 44].

For kinetic transport equations however a MFPT theory is missing. These
equations were introduced in biological modeling in the 1980s, and have become
a powerful tool in the study of ecological and cell movement [2, 28, 26, 53,
29, 50, 52]. They are particularly useful when movement characteristics of
single individuals can be measured such as the velocities, turning rates and
directional preferences of migrating cells or animals. Recent advances in
microscopy and communication systems have facilitated the merging of
theoretical results with actual data, as it is now possible to observe
molecular or cellular motion at high resolution and to follow in detail the
movement of animals over large distances or underwater. It is therefore of
great interest to develop a MFPT theory for kinetic equations as first passage
computations may serve as a bridge between theoretic estimates that depend on
microscopic features of the transport process and the average timescales
observed experimentally or in the field.

Here, we are interested in animal or cellular species that orient their
movement along directional features of the environment. We derive expressions
for the mean first arrival time to a target as a function of the initial
position and velocity. Broadly speaking, since directional motion is
anisotropic, the mean first arrival time will depend on the microscopic
features of the transport process, such as turning rates, initial speed and
direction of movement. For the most general case we are not able to find a
closed-form equation for the MFPT. However, for diffusive and anisotropic
transport we are able to write integral equations of the first kind for the
MFPT starting from an initial position and velocity. In special cases, for
example for anisotropic, symmetric motion under the bimodal von-Mises
distribution, this general expression leads to a differential equation for the
MFPT that can be solved analytically and that can be readily applied to glioma
movement [51, 32, 62, 29, 17, 15], sea turtle orientation [52] and hill
topping of butterflies [50]. Our main results are contained in (31) and (37);
the analytically solvable variant of (37) for the special case of the bimodal
von-Mises distribution is shown in (51).

The rest of this paper is organized as follows. In Section 2 we introduce
kinetic transport equations for biological modelling, and we categorize them
as of Boltzmann type, diffusive type and anisotropic type. We consider the
forward and backward formulations, define the survival probability, and write
its dynamics using the backward equation. We thus derive the MFPT equation for
the expected exit time $\Theta(x,v)$ of a random walker starting at $(x,v)$
for the diffusive type in Section 3 and for the anisotropic case in Section 4.
The main results are given in formulas (31) and (37), respectively. The
parabolic scaling in Section 5 allows us to consider the limit in cases of
macroscopic spatio-temporal scales. In this scaling we derive an anisotropic
MFPT for the mean exit time $\Theta(x,v)=T(x)$ which no longer depends on the
initial velocity as shown in (51). In Section 6 we apply our results to
various cases. For anisotropic movement in a radially symmetric circular
domain we evaluate the MFPT to reach given boundaries in the disk and in the
annulus. A second application considers oriented movement in a habitat that is
criss-crossed by linear features. This scenario has been used to explain wolf
movement in presence of seismic lines in a forest landscape in Northern Canada
by [41]. Here we show how the analysis carried out by [41] fits within our
more general framework. We close with a conclusion Section 7 in which we
relate our results to existing methods and open the door to interesting
forthcoming problems.

##  2 Transport Equations

Transport equations describe the time evolution of a particle density
$p(x,v,t)$ at time $t\geq 0$, location $x\in\Omega\subset\mathbb{R}^{n}$ and
velocity $v\in V\subset\mathbb{R}^{n}$. For simplicity, we take
$V=[\sigma_{1},\sigma_{2}]\times\mathbb{S}^{n-1}$, where the speeds
$\sigma_{1},\sigma_{2}$ obey $0\leq\sigma_{1}\leq\sigma_{2}<\infty$. The time
evolution of $p(x,v,t)$ is described by the forward transport equation

(1) $p_{t}(x,v,t)+v\cdot\nabla_{x}p(x,v,t)={\cal L}p(x,v,t)\,,$

where the subscript $t$ denotes the partial time derivative, the subscript $x$
denotes the spatial derivative, and ${\cal L}$ is the turning operator, which
describes the specific directional changes of the particles. Typically, ${\cal
L}$ is defined via an integral representation

(2) ${\cal
L}\varphi(x,v,t)=-\mu(x)\varphi(x,v,t)+\mu(x)\int_{V}K(x,v,v^{\prime})\varphi(x,v^{\prime},t)dv^{\prime}\,$

where $\varphi(x,v,t)$ is a test function. The first term on the right hand
side of (2) is the rate at which particles at position $x$ change their
velocity from $v$ to any other velocity $v^{\prime}\in V$; the second term is
the rate at which particles at location $x$ switch into velocity $v$ from any
other velocity $v^{\prime}$, The quantity $\mu(x)$ is the turning rate and its
inverse is typically identified as the mean run time at position $x$. Although
general formulations that include spatially dependent turning rates $\mu(x)$
are possible, for simplicity we assume $\mu(x)=\mu$ to be spatially uniform in
this paper. The turning kernel $K(x,v,v^{\prime})$ denotes the probability
density of switching velocity from $v^{\prime}$ to $v$, given that a turn
occurs at location $x$. The properties of $K(x,v,v^{\prime})$ are key to the
analysis presented in this work. The minimal assumptions on $K$ are

(3) $K(x,v,v^{\prime})\geq 0,\quad K(x,\cdot,\cdot)\in L^{2}(V\times
V),\quad\int_{V}K(x,v,v^{\prime})dv=1.$

The second condition implies that the integral operator with kernel $K$ is a
compact Hilbert-Schmidt operator in $L^{2}(V)$ [27], and the last assumption
ensures that during velocity changes no particle is lost. Among the various
possible forms of $K(x,v,v^{\prime})$ are those that allow to represent
collisions between particles, periods of straight motion alternating with
abrupt directional changes (run-and-tumble models) and general jump velocity
processes where a particle’s velocity undergoes a series of discrete changes
[28, 46, 53, 10, 13]. Roughly, we can classify the most important cases into
three categories as shown below.

###  2.1 Boltzmann equations:

Transport equations, such as (1) have a long history in the study of the
thermodynamics of gases which include the effects of collisions between
particles [12]. In general formulations of the Boltzmann equation the operator
${\cal L}$ contains a non-linear collision kernel and the set of velocities
$V$ is unbounded. A linear kernel such as (2) arises as linearization of the
collision kernel at a Maxwellian equilibrium distribution [12]. An important
difference with respect to biological applications is that the Boltzmann
equation conserves mass, momentum and energy. i.e. the kernel of the operator
${\cal L}$ is five dimensional. In biological applications there is only one
conserved quantity, mass, and only in cases where there is no particle birth
or death.

###  2.2 Diffusive transport equations:

The diffusive transport equation is characterized by a turning operator ${\cal
L}$ whose null space consists of functions that are constant in velocity. This
means that at equilibrium there are no preferred directions and that over long
time scales the dynamics becomes diffusion-like. Full theories of diffusive
transport equations have been developed and applied to biological processes
such as chemotaxis [28, 46]. We derive the MFPT equation for diffusive
transport in Section 3.

###  2.3 Anisotropic transport equations:

Here, particles have strong orientation guidance from the underlying
environment. Mathematically this occurs via a non-trivial one-dimensional
kernel of the turning operator ${\cal L}$. The MFPT equation for anisotropic
transport is discussed in Section 3. We discuss this scenario in detail since
it is one of the key cases in this work. Under anisotropic transport it is
assumed that the turning kernel $K$ does not depend on the incoming velocity
$v^{\prime}$, i.e.

(4) $K(x,v,v^{\prime})=q(x,v),$

where $q\geq 0$. This assumption may seem restrictive, since many species
would show some form of persistence in movement direction and it is entirely
possible to keep the dependence on $v^{\prime}$ in anisotropic movement as
well. However, it has been shown in many applications that the simplifying
assumption (4) is extremely useful in the modelling process and it can be
justified in many cases, for example, when particles change direction using a
well-defined underlying network structure. The anisotropic framework was
developed in [26] and extended in [49]. To build $q(x,v)$ we assume that a
distribution of preferred directions $\tilde{q}(x,\theta)$ is given for every
spatial position $x$, where $\theta\in\mathbb{S}^{n-1}$ is a unit vector. We
impose

$\tilde{q}(x,\theta)\geq
0,\qquad\mbox{and}\qquad\int_{\mathbb{S}^{n-1}}\tilde{q}(x,\theta)d\theta=1$

and assume $q(x,v)$ is proportional to $\tilde{q}(x,\hat{v})$, where
$\hat{v}=v/\|v\|$ denotes the corresponding unit vector. Since $\tilde{q}$ is
a probability distribution on $\mathbb{S}^{n-1}$, and $q$ a probability
distribution on the set of velocities $V$, we need to normalize appropriately
leading to

(5) $q(x,v)=\frac{\tilde{q}(x,\hat{v})}{\omega},\quad\mbox{ with
}\quad\omega\equiv\int_{V}\tilde{q}(x,\hat{v})dv=\left\\{\begin{array}[]{ll}\frac{1}{n}(\sigma_{2}^{n}-\sigma_{1}^{n})&\mbox{for
}\sigma_{1}<\sigma_{2};\\\\[4.0pt] \sigma^{n-1}&\mbox{for
}\sigma_{1}=\sigma_{2}=\sigma.\end{array}\right.$

We refer to the rescaled quantity $q(x,v)$ as the directional distribution of
the underlying environment, although the true directional distribution is
$\tilde{q}(x,\hat{v})$. For this choice of turning kernel, equation (2)
simplifies to

${\cal
L}\varphi(v)=\mu(x)(q(x,v)\hat{\varphi}-\varphi(v)),\qquad\mbox{with}\quad\hat{\varphi}\equiv\int_{V}\varphi(v^{\prime})dv^{\prime}.$

To summarize, the anisotropic transport equation has the form

(6) $p_{t}+v\cdot\nabla
p=\mu(x)(q(x,v)\hat{p}-p)\qquad\mbox{with}\quad\hat{p}\equiv\int_{V}p(t,x,v^{\prime})dv^{\prime}.$

For the following analysis it is useful to consider two statistical
quantities, the expectation $\mathbb{E}_{q}(x)$ and the variance-covariance
matrix $\mathbb{V}_{q}(x)$ of $q(x,v)$ on $V$:

(7)
$\mathbb{E}_{q}(x)=\int_{V}vq(x,v)dv,\qquad\mathbb{V}_{q}(x)=\int_{V}(v-\mathbb{E}_{q}(x))(v-\mathbb{E}_{q}(x))^{T}q(x,v)dv.$

We now discuss a few examples of directional distributions $q$ by assuming
constant particle speed $\sigma$, i.e. $V=\sigma\mathbb{S}^{n-1}$. The
extension to $V=[\sigma_{1},\sigma_{2}]\times\mathbb{S}^{n-1}$ is
straightforward, but it introduces some extra parameters that need to be
carried through the integrals.

  1. 2.3.1

Uniform distribution: The simplest case arises if there is no directional bias
at all. This, of course, is also the most fundamental type of motion in
diffusive transport. We write

$\tilde{q}(x,\theta)=\frac{1}{|\mathbb{S}^{n-1}|},\qquad
q(x,v)=\frac{1}{\omega|\mathbb{S}^{n-1}|}.$

leading to

$\mathbb{E}_{q}(x)=0,\qquad\mathbb{V}_{q}(x)=\frac{\sigma^{2}}{\omega
n}{\mathbb{I}},$

where $\mathbb{I}$ denotes the identity matrix.

  2. 2.3.2

Strict alignment: If the underlying medium favors motion along a given
direction $\gamma(x)\in\mathbb{S}^{n-1}$, as for example along a fiber,
tissue, microtubule or track, then

$\tilde{q}(x,\theta)=\delta(\gamma(x)-\theta),\qquad
q(x,v)=\frac{1}{\omega}\delta(\gamma(x)-\hat{v}).$

In this case

$\mathbb{E}_{q}(x)=\sigma\gamma(x),\qquad\mathbb{V}_{q}(x)=0.$

  3. 2.3.3

von-Mises distribution: The von-Mises distribution is the analog of the normal
distribution on the unit sphere $\mathbb{S}^{n-1}$. Strictly speaking, it is
called von-Mises distribution only on $\mathbb{S}^{1}$. In higher dimensions
it is called the Fisher distribution. In two-dimensions the von-Mises
distribution is [30]

(8) $\tilde{q}(x,\theta)=\frac{1}{2\pi
I_{0}(k(x))}e^{k(x)\gamma(x)\cdot\theta},$

where $k(x)\geq 0$ is a measure of concentration of the distribution,
$\gamma(x)\in\mathbb{S}^{1}$ is a given preferred direction, and $I_{0}(k(x))$
is the modified Bessel function of first kind of order $0$. As $k(x)\to 0$ the
von-Mises distribution becomes uniform (as in 2.3.1) whereas for
$k(x)\to\infty$ it becomes singular (as in 2.3.2). The expectation and
variance of the two-dimensional von-Mises distribution are [30]

(9a) $\displaystyle\mathbb{E}_{q}(x)$ $\displaystyle=$
$\displaystyle\sigma\frac{I_{1}(k(x))}{I_{0}(k(x))}\gamma(x),$ (9b)
$\displaystyle\mathbb{V}_{q}(x)$ $\displaystyle=$
$\displaystyle\sigma^{2}\left[\frac{1}{2}\left(1-\frac{I_{2}(k(x))}{I_{0}(k(x))}\right)\mathbb{I}_{2}+\left(\frac{I_{2}(k(x))}{I_{0}(k(x))}-\left(\frac{I_{1}(k(x))}{I_{0}(k(x))}\right)^{2}\right)\gamma(x)\gamma(x)^{T}\right]\,,$

where $I_{0},I_{1},I_{2}$ denote the modified Bessel functions of first kind
of order 0,1 and 2, respectively.

  4. 2.3.4

Bimodal von-Mises distribution: In many cases $q$ is symmetric, for example
when guided motion along a two-dimensional fiber or a track has the same
probability along each direction. In this case $q(x,-v)=q(x,v)$ and the
bimodal von-Mises distribution can be written as the average of two single
direction terms

(10) $\tilde{q}(x,\theta)=\frac{1}{4\pi
I_{0}(k(x))}\Bigl{(}e^{k(x)\gamma(x)\cdot\theta}+e^{-k(x)\gamma(x)\cdot\theta}\Bigr{)},$

with expectation and variance

(11a) $\displaystyle\mathbb{E}_{q}(x)$ $\displaystyle=$ $\displaystyle 0\,,$
(11b) $\displaystyle\mathbb{V}_{q}(x)$ $\displaystyle=$
$\displaystyle\sigma^{2}\left[\frac{1}{2}\left(1-\frac{I_{2}(k(x))}{I_{0}(k(x))}\right)\mathbb{I}_{2}+\frac{I_{2}(k(x))}{I_{0}(k(x))}\gamma(x)\gamma(x)^{T}\right]\,.$

Note that instead of using a dyadic product $\gamma\gamma^{T}$ some authors
prefer to use a tensor product $\gamma\otimes\gamma$. We use the two
interchangeably in this paper.

###  2.4 Backward Transport Equation

The Kolmogorov backward transport equation is the foundation from which the
mean first passage time equation is derived. For clarity we distinguish
initial location and velocity $(x,v)$ from terminal location and velocity
$(\bar{x},\bar{v})$. Thus $p(\bar{x},\bar{v},t|x,v)$, denotes the probability
density of a random walker starting at location
$x\in\Omega\subset\mathbb{R}^{n}$ with initial velocity $v\in V$ to be at
location $\bar{x}$ with velocity $\bar{v}\in V$ at time $t\geq 0$. Let us
rewrite the forward equation

(12)
$p_{t}+v\cdot\nabla_{\bar{x}}p=-\mu\,p+\mu\int_{V}K(\bar{x},\bar{v},\bar{v}^{\prime})p(x,\bar{v}^{\prime},t|x,v)d\bar{v}^{\prime},$

where the gradient subscript allows to distinguish between the spatial
derivative with respect to the current location ($\nabla_{\bar{x}}$) and the
derivative with respect to the starting location ($\nabla$ without index). The
model (12) can be seen as a forward equation of the stochastic process of a
velocity jump process [45]. To derive the corresponding Kolmogorov backward
equation, we consider the formal adjoint. We write (12) as

$p_{t}(\bar{x},\bar{v},t|x,v)={\cal A}p(\bar{x},\bar{v},t|x,v).$

The operator $\cal A$ acts on functions $\varphi(\bar{x},\bar{v})$ of the
variables $(\bar{x},\bar{v})$ in $\Omega\times V$ and is defined as

${\cal
A}\varphi(\bar{x},\bar{v})=-\bar{v}\cdot\nabla_{\bar{x}}\varphi(\bar{x},\bar{v})-\mu\,\varphi(\bar{x},\bar{v})+\mu\int_{V}K(\bar{x},\bar{v},\bar{v}^{\prime})\varphi(\bar{x},\bar{v}^{\prime})d\bar{v}^{\prime}.$

To proceed, recall that the Chapman-Kolmogorov equation [63], under the
assumption of time translational invariance, yields

$p(\bar{x},\bar{v},t|x,v)=\int_{\Omega}\int_{V}p(\bar{x},\bar{v},\tau|x^{\prime},v^{\prime})p(x^{\prime},v^{\prime},t-\tau|x,v)dv^{\prime}dx^{\prime},$

for any $\tau\geq 0$. Taking the time derivative on both sides we find

$\displaystyle p_{t}(\bar{x},\bar{v},t|x,v)$ $\displaystyle=$
$\displaystyle\int_{\Omega}\int_{V}p(\bar{x},\bar{v},\tau|x_{1},v_{1}){\cal
A}p(x^{\prime},v^{\prime},t-\tau|x,v)dv^{\prime}dx^{\prime}$ $\displaystyle=$
$\displaystyle\int_{V}\int_{\Omega}{\cal
A^{*}}p(\bar{x},\bar{v},\tau|x^{\prime},v^{\prime})p(x^{\prime},v^{\prime},t-\tau|x,v)dx^{\prime}dv^{\prime},$

where the operator $\cal A^{*}$ is the adjoint of $\cal A$ and acts on the
variables $(x^{\prime},v^{\prime})$ in $\Omega\times V$. Notice that if we let
$\tau\to t$ then $\lim_{\tau\to
t}p(x^{\prime},v^{\prime},t-\tau|x,v)=p(0,x^{\prime},v^{\prime}|x,v)=\delta(x^{\prime}-x)\delta(v^{\prime}-v)$
leading to

(13) $p_{t}(\bar{x},\bar{v},t|x,v)={\cal A^{*}}p(\bar{x},\bar{v},t|x,v),$

where the adjoint $\cal A^{*}$ acts on the $(x,v)$ values in $\Omega\times V$.
This is the Kolmogorov backward equation. To find $\cal A^{*}$ explicitly,
given $\cal A$ on $\Omega\times V$, we consider two functions
$\phi(x,v),\psi(x,v)$ and compute

$\displaystyle\int_{\Omega}\int_{V}\phi{\cal A}\psi\,dv\,dx$ $\displaystyle=$
$\displaystyle-\int_{\Omega}\int_{V}\phi(x,v)(v\cdot\nabla_{x})\psi(x,v)\,dv\,dx$
$\displaystyle+\int_{\Omega}\int_{V}\phi(x,v)\left(-\mu\psi(x,v)+\mu\int_{V}K(x,v,v^{\prime})\psi(x,v^{\prime})dv^{\prime}\right)dv\,dx$
$\displaystyle=$
$\displaystyle\int_{V}\int_{\Omega}\psi(x,v)(v\cdot\nabla_{x})\phi(x,v)dx\,dv-\mu\int_{V}\int_{\Omega}\phi(x,v)\psi(x,v)dx\,dv$
$\displaystyle+\mu\int_{V}\int_{\Omega}\left(\int_{V}K(x,v,v^{\prime})\phi(x,v)dv^{\prime}\right)\psi(x,v^{\prime})\,dx\,dv.$

The last identity follows from integration by parts and through the boundary
condition $\phi(x,v)=\psi(x,v)=0$ for $x\in\partial\Omega$. If we now denote
the adjoint integral kernel as

$K^{*}(x,v,v^{\prime})=K(x,v^{\prime},v),$

then $\cal A^{*}$, the adjoint operator of $\cal A$, is given by

${\cal
A}^{*}\phi(x,v)=(v\cdot\nabla)\phi(x,v)-\mu\phi(x,v)+\mu\int_{V}K^{*}(x,v,v^{\prime})\phi(x,v^{\prime})dv^{\prime}.$

The Kolmogorov backward equation in explicit form becomes

(14) $p_{t}-v\cdot\nabla
p=-\mu\,p+\mu\int_{V}K^{*}(x,v,v^{\prime})p(\bar{x},\bar{v},t|x,v^{\prime})dv^{\prime},$

This equation was previously derived by Stroock [61]. We can write (14) to
mirror (1) by introducing the turning operator ${\cal L}^{*}$

(15) $p_{t}(t,x,v)-v\cdot\nabla p(t,x,v)={\cal L}^{*}p(t,x,v),$

where ${\cal L}^{*}={\cal A}^{*}-(v\cdot\nabla)$ is given by

(16) ${\cal
L}^{*}\varphi(t,x,v)=-\mu\varphi(t,x,v)+\mu\int_{V}K^{*}(x,v,v^{\prime})\varphi(t,x,v^{\prime})dv^{\prime}.\,$

It is straightforward to show that ${\cal L}^{*}$ is the formal adjoint of
${\cal L}$ defined in (2).

###  2.5 The Survival Probability

To derive an expression for the mean first passage time of the transport
process we first introduce the survival probability $s(x,v,t)$. This quantity
is defined as the probability that having started from an initial position
$x\in\Omega$ with velocity $v\in V$ the final position $\bar{x}$ and velocity
$\bar{v}$ at time $t$ remain in the same domain $\Omega\times V$. We thus
assume $p(\bar{x},\bar{v},t=0|x,v)=\delta(\bar{x}-x)\delta(\bar{v}-v)$ with
$x\in\Omega$ and $v\in V$ and define

(17)
$s(x,v,t)=\int_{V}\int_{\Omega}p(\bar{x},\bar{v},t|x,v)\,d\bar{x}\,d\bar{v}\qquad\forall
x\in\Omega,\forall v\in V.$

In principle the domain for the forward and backward equations can be
different, as one can study the exit time from a subset of the domain on which
the forward equation is valid. Here, for simplicity we assume they are the
same and utilize absorbing boundary conditions on $\Omega$. Thus, we set
$s(x,v,t)=0$ for any $x\in\partial\Omega$ to indicate that a random walker
starting at the boundary of $\Omega$ does not survive at any time.
Furthermore, the initial condition
$p(\bar{x},\bar{v},t=0|x,v)=\delta(\bar{x}-x)\delta(\bar{v}-v)$ implies that
$s(x,v,t=0)=1$ for all $x\in\Omega$. We also define the probability flux
$j(x,v,t)$ as

(18)
$j(x,v,t)=\int_{V}\int_{\Omega}v\,p(\bar{x},\bar{v},t|x,v)\,d\bar{x}\,d\bar{v}=vs(x,v,t).$

We now integrate the backward equation (14) over the final position
$\bar{x}\in\Omega$ and $\bar{v}\in V$ to find a transport equation for
$s(x,v,t)$,

(19) $\displaystyle s_{t}-v\cdot\nabla s=-\mu
s+\mu\int_{V}K^{*}(x,v,v^{\prime})s(x,v^{\prime},t)dv^{\prime}.$

Given the above definitions of $s(x,v,t)$ and $j(x,v,t)$ we note that one is
most often concerned with the survival probability, and the mean first passage
time, of reaching the boundary $\partial\Omega$ regardless of the velocity.
Hence we introduce the marginal survival probability $S(x,t)$ defined as

(20) $S(x,t)=\frac{1}{|V|}\int_{V}s(x,v,t)\,dv\,,$

where $|V|=\int_{V}dv$. The corresponding flux is

(21) $J(x,t)=\frac{1}{|V|}\int_{V}j(x,v,t)\,dv.$

Note that $S(x,t)$ remains a probability and $S(x,t=0)=1$. The flux $J(x,t)$
still carries units of a velocity. For later use we collect the initial
conditions here:

(22) $p(\bar{x},\bar{v},0|x,v)=\delta(\bar{x}-x)\delta(\bar{v}-v),\qquad
s(x,v,0)=1,\qquad S(x,0)=1.$

To study the dynamics of $S$ and $J$ we must specify the properties of the
turning kernel $K(x,v,v^{\prime})$ and of its adjoint $K^{*}(x,v,v^{\prime})$;
we do so in the two cases below for diffusive and anisotropic transport.

##  3 MFPT for Diffusive Transport

We now write the equation for the mean first passage time under the diffusive
transport regime discussed in 2.2. In addition to (3) we assume symmetry and
particle conservation of the turning kernel $K$

(23)
$K(x,v,v^{\prime})=K(x,v,-v^{\prime}),\qquad\mbox{and}\qquad\int_{V}K(x,v,v^{\prime})dv^{\prime}=1.$

For the adjoint operator this implies

(24)
$K^{*}(x,v,v^{\prime})=K^{*}(x,-v,v^{\prime}),\qquad\mbox{and}\qquad\int_{V}K^{*}(x,v,v^{\prime})dv=1.$

We use assumption (24) and integrate equation (19) for the survival
probability $s(x,v,t)$ over velocity space to obtain a conservation law

(25) $\displaystyle S_{t}-\nabla\cdot J=0.$

We also multiply the backward equation (19) by $v$, integrate over velocity
space, and divide by the volume $|V|$ to obtain

(26) $\displaystyle J_{t}-\frac{1}{|V|}\int_{V}v(v\cdot\nabla)s\,dv=-\mu
J+\frac{\mu}{|V|}\int_{V}\int_{V}v\,K^{*}(x,v,v^{\prime})s(x,v^{\prime},t)\,dv^{\prime}\,dv.$

Since assumption (24) on the symmetry of $T$ implies that

$\int_{V}vK^{*}(x,v,v^{\prime})dv=0,$

the last integral vanishes. Finally, using tensor notation for a dyadic
product of two vectors we write

(27) $J_{t}-\frac{1}{|V|}\nabla\cdot\int_{V}\int_{V}v\otimes
v\,s(x,v,t)\,dv=-\mu J.$

Equations (25) and (27) do not not lead to a closed system of equations for
$S$ and $J$. We will discuss the closure problem in more detail later. For
now, we identify the complement of $S(x,t)$, the function $F(x,t)=1-S(x,t)$,
as the exit probability from $\Omega$. Its density is given by

$f(x,t)=F_{t}(x,t)=-S_{t}(x,t).$

The expected survival time $T(x)$ before exiting $\Omega$ given the random
walker started at the initial position $x\in\Omega$ can be now computed as

$T(x)\equiv\int_{0}^{\infty}tf(x,t)dt=\int_{0}^{\infty}S(x,t)dt.$

We refer to $T(x)$ as the mean first passage time. Integrating (25) for the
dynamics of $S$, and (26) for the dynamics of $J$, both over time, and using
the initial conditions (22), we obtain the following

(28a) $\displaystyle-1$ $\displaystyle=$ $\displaystyle\nabla\cdot\int Jdt$
(28b) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{1}{|V|}\nabla\cdot\int_{V}\int_{0}^{\infty}v\otimes
v\,s(x,v,t)\,dt\,dv-\mu\int Jdt.$

where it is assumed that $p(\bar{x},\bar{v},t\to\infty|x,v)=0$ and that
$\int_{V}vdv=0$, by symmetry of the $V$ domain. We can now take the divergence
of (28b) and substitute the resulting expression in (28a) to obtain

(29) $-1=\frac{1}{\mu|V|}\nabla\otimes\nabla:\int_{V}\int_{0}^{\infty}v\otimes
v\,s(x,v,t)\,dt\,dv.$

Here, we have used the colon notation $:$ to denote the tensor convolution of
two vectors $a$ and $b$ as

$a\otimes a:b\otimes b=\sum_{i,j=1}^{n}a_{i}a_{j}b_{i}b_{j}.$

It is useful to introduce $\Theta(x,v)$, the expected survival time before
leaving $\Omega$ given the random walker started at position $x$ with velocity
$v$. Thus, through $\Theta(x,v)$ we keep the initial velocity explicit in the
definition of the MFPT. We write

(30)
$\Theta(x,v)\equiv\int_{0}^{\infty}ts_{t}(x,v,t)dt=\int_{0}^{\infty}s(x,v,t)dt.$

Integrating (30) over $V$ we obtain

$T(x)=\frac{1}{|V|}\int_{V}\Theta(x,v)dv.$

Using equation (29) we find the main result of this section, which is the MFPT
equation for the case of diffusive transport

(31) $-1=\frac{1}{\mu|V|}\nabla\otimes\nabla:\int_{V}v\otimes
v\,\Theta(x,v)dv.$

###  3.1 Example:

We assume constant speed so that $V=\sigma{\mathbb{S}}^{n-1}$ and rotational
symmetry in the initial velocity dependence of $s(x,v,t)=s(x,t)$ so that

$\int_{V}v\otimes
v\,s(x,v,t)dv=\frac{\sigma^{2}}{n}{\mathbb{I}}\int_{V}s(x,t)dv,$

where $\mathbb{I}$ denotes the identity matrix. Substituting this expression
in (31) we find

(32) $-1=\frac{\sigma^{2}}{n\mu}\Delta T(x).$

This equation has the well known form of a MFPT for a diffusion process in an
$n$-dimensional system with diffusion constant
$D_{n}=\frac{\sigma^{2}}{n\mu}$. This coefficient also agrees with the
parabolic limit of the transport equation, as carried out in several papers,
for example [28].

##  4 MFPT for Anisotropic Transport

We now derive the equation for the mean first passage time under anisotropic
transport, where directional cues from the environment bias particle movement
along one or more select directions. We discussed anisotropic transport in
Section 2.3. For simplicity, in addition to (3) we also assume that
$K(x,v,v^{\prime})=q(x,v)$ is independent of the incoming velocity
$v^{\prime}$ and that $q(x,v)$ is symmetric in $v$. We summarize the
properties of $K(x,v,v^{\prime})$ and $q(x,v)$ below

(33a) $\displaystyle K(x,v,v^{\prime})=q(x,v),\qquad\mbox{i.e.}\qquad
K^{*}(x,v,v^{\prime})=q(x,v^{\prime}),$ (33b) $\displaystyle q\geq 0,\qquad
q\in L^{2}(V),\qquad\int_{V}q(x,v)dv=1,\qquad q(x,-v)=q(x,v).$

This choice of $K(x,v,v^{\prime})$ does not satisfy the conditions listed in
(24), since

$\int_{V}K^{*}(x,v,v^{\prime})dv=q(x,v^{\prime})|V|,$

and the right hand side equals $1$ if and only if
$q(x,v^{\prime})=\frac{1}{|V|}$. Hence, the results shown in Section 3 for
diffusive transport apply only to the simple case of a uniform turning kernel
but are not valid for a general, non uniform $q(x,v)$. To find an expression
for the MFPT under anisotropic transport we start with the backward equation
(19) for the survival probability $s(x,v,t)$ adapting it to the case at hand
and write

(34) $\displaystyle s_{t}-v\cdot\nabla s=-\mu
s+\mu\int_{V}q(x,v^{\prime})s(x,v^{\prime},t)dv^{\prime}.$

We now multiply equation (34) by $q(x,v)$ and integrate over the velocity
domain to obtain

$\displaystyle\int_{V}q(x,v)s_{t}(x,v,t)dv-\int_{V}vq(x,v)\cdot\nabla
s(x,v,t)dv$
$\displaystyle=-\mu\int_{V}q(x,v)s(x,v,t)dv+\mu\int_{V}\int_{V}q(x,v)q(x,v^{\prime})s(x,v^{\prime},t)dv^{\prime}dv=0.$

The last identity is derived from the property that $q(x,v)$ integrates to
unity in velocity space. To simplify the notation we introduce the following
two quantities

$n(x,t)=\int_{V}q(x,v)s(x,v,t)dv\quad\mbox{and}\quad
w(x,t)=\int_{V}vq(x,v)\cdot\nabla s(x,v,t)dv,$

which leads to

(35) $n_{t}-w=0.$

We now derive a differential equation for $w$, including dependencies only
where explicitly needed

$\displaystyle w_{t}$ $\displaystyle=$ $\displaystyle\int_{V}vq\cdot\nabla
s_{t}\,dv=\int_{V}vq\cdot\nabla\left(v\cdot\nabla s-\mu
s+\mu\int_{V}q(x,v^{\prime})s(x,v^{\prime},t)dv^{\prime}\right)dv$
$\displaystyle=$ $\displaystyle\int_{V}q\,(v\otimes
v:\nabla\otimes\nabla)\,s\,dv-\mu\int_{V}vq\cdot\nabla
s\,dv+\mu\int_{V}vq(x,v)dv\cdot\nabla\int_{V}q(x,v^{\prime})s(x,v^{\prime},t)\,dv^{\prime}.$

Since $\int_{V}vq(x,v)dv=0$ by the symmetry assumption in (33b), the last term
vanishes yielding

(36) $w_{t}=\int_{V}q(x,v)(v\otimes v:\nabla\otimes\nabla)s(x,v,t)dv-\mu w.$

We now integrate (35) and (36) over time and obtain boundary terms. Since no
particle survives forever, i.e. $s(x,v,t\to\infty)=0$, we can write
$n(x,\infty)=w(x,\infty)=0.$ Moreover, the initial condition $s(x,v,0)=1$
yields $n(x,0)=1$ and $w(x,0)=0.$ Using these identities we find

$\displaystyle-1$ $\displaystyle=$ $\displaystyle\int_{0}^{\infty}w(x,t)dt$
$\displaystyle 0$ $\displaystyle=$
$\displaystyle\int_{V}\int_{0}^{\infty}q(x,v)(v\otimes
v:\nabla\otimes\nabla)s(x,v,t)dt\,dv-\mu\int_{0}^{\infty}w(x,t)dt,$

which combine to yield

$-1=\frac{1}{\mu}\int_{V}\int_{0}^{\infty}q(x,v)(v\otimes
v:\nabla\otimes\nabla)s(x,v,t)dv\,dt.$

We can rewrite this expression by invoking $\Theta(x,v)$, the MFPT where both
the initial condition and velocity are specified, to finally obtain

(37) $-1=\frac{1}{\mu}\int_{V}q(x,v)(v\otimes
v:\nabla\otimes\nabla)\Theta(x,v)dv.$

Equation (37) is an interesting anisotropic integro-differential equation for
the mean first exit time $\Theta(x,v)$ when starting at $(x,v)$. In this
general setting, we cannot derive a simple equation for the mean exit time
$T(x)$, since the initial velocity cannot be neglected. Also, we cannot simply
integrate (37) over velocity space, since the velocity and derivative terms
are mixed. This issue is similar to a moment closure problem, which is well
known and well studied in the theory of transport equations [11, 25]. In order
to proceed we must find an approximation for the $q(x,v)$-modulated second
moment of $\Theta(x,v)$. One way to accomplish this is to consider parabolic
scaling, which leads to an equation for $T(x)$ alone. We do this in the next
section.

##  5 Parabolic Scaling for Anisotropic Transport

Parabolic scaling is a well known technique to study transport processes on
macroscopic time and space scales. While the scaling of the forward equation
is standard, the corresponding scaling for the backward equation is a new
result to the best of the authors’ knowledge. Here, for consistency, we
present both cases, forward and backward. For simplicity we assume that
$K(x,v,v^{\prime})=q(x,v)=q(x,-v)$ is a symmetric, anisotropic directional
distribution that does not depend on the incoming velocity $v^{\prime}$ and
whose properties are listed in (33b). We now collect some functional
analytical properties of $K(x,v,v^{\prime})=q(x,v)$. First, for the forward
transport equation in (1) we write the turning operator $\cal L$ defined via
the integral representation (2) as a mapping on $L^{2}(V)$ so that for a test
function $\varphi(x,v)$

$\mathcal{L}\varphi(x,v)=-\mu\varphi(x,v)+\mu
q(x,v)\hat{\varphi}(x),\qquad\hat{\varphi}(x)=\int_{V}\varphi(x,v^{\prime})dv^{\prime}.$

Similarly, for the backward Kolmogorov equation in (15) the turning operator
${\cal L}^{*}$ can be written as

$\mathcal{L}^{*}\varphi(x,v)=-\mu\varphi(x,v)+\mu\int_{V}q(x,v^{\prime})\varphi(x,v^{\prime})dv^{\prime}.$

To perform the parabolic scaling we need to find the null spaces (kernel) and
the pseudo-inverses of the ${\cal L,L^{*}}$ operators on the complement of
their null spaces. We collect these properties in the following Lemmas.

######  Lemma 5.1.

The kernels of the operators $\mathcal{L}$ and $\mathcal{L}^{*}$ are the one-
dimensional spaces spanned by the functions $q(x,v)$ and 1, respectively.

$\mbox{ker}\,\mathcal{L}=\langle
q(x,v)\rangle,\qquad\mbox{ker}\,\mathcal{L}^{*}=\langle 1\rangle.$

######  Proof 5.2.

To find the kernel of $\mathcal{L}$ we seek functions $\phi(x,v)$ that satisfy
$\mathcal{L}\phi=0$. This condition is

$-\mu\phi(x,v)+\mu
q(x,v)\hat{\phi}(x)=0,\quad\hat{\phi}(x)\equiv\int_{V}\phi(x,v)dv.$

To solve the above equations consistently, $\phi(x,v)$ must be proportional to
$q(x,v)$. To find the kernel of the adjoint ${\cal L}^{*}$ we must instead
solve $\mathcal{L}^{*}\phi=0$, or

$-\mu\phi(x,v)+\mu\int_{V}q(x,v^{\prime})\phi(x,v^{\prime})dv^{\prime}=0,$

which instead implies that $\phi(x,v)$ must be a constant.

######  Lemma 5.3.

For the ${\cal L}$ operator consider the inverse of the generating function
$q(x,v)$ as a weight in $L^{2}(V)$ and restrict ${\cal L}$ to the domain of
functions that are $L^{2}(V)$ when weighted by $q^{-1}$, which we denote by
$L^{2}_{q^{-1}}(V)$. A similar construct for ${\cal L}^{*}$ using the
corresponding generating function $1$ does not alter the domain of ${\cal
L}^{*}$. Then, the pseudo-inverse operators acting on

$\mathcal{L}:L^{2}_{q^{-1}}(V)\mapsto
L^{2}_{q^{-1}}(V),\qquad\mathcal{L}^{*}:L^{2}(V)\mapsto L^{2}(V)$

are, respectively, the multiplication operators

$\left(\mathcal{L}|_{\langle
q\rangle^{\perp}}\right)^{-1}=-\frac{1}{\mu},\qquad\left(\mathcal{L}^{*}|_{\langle
1\rangle^{\perp}}\right)^{-1}=-\frac{1}{\mu}.$

######  Proof 5.4.

To find ${\cal L}^{-1}$ the pseudo-inverse of $\mathcal{L}$ on
$\left(\mbox{ker}\,\mathcal{L}\right)^{\perp}$, the orthogonal complement of
$\mbox{ker}\,\mathcal{L}$, we consider a given
$\phi\in\mbox{ker}\,\mathcal{L}^{\perp}$ and find
$\psi\in\mbox{ker}\,\mathcal{L}^{\perp}$ that solves $\mathcal{L}\psi=\phi$,
so that $\psi={\cal L}^{-1}\phi$. The equation $\mathcal{L}\psi=\phi$ reads

(38) $-\mu\psi(x,v)+\mu q(x,v)\hat{\psi}(x)=\phi(x,v).$

The assumption that $\psi\in\langle q\rangle^{\perp}$ in $L^{2}_{q^{-1}}(V)$
implies that

$0=(\psi,q)_{q^{-1}}=\int_{V}\psi(x,v)q(x,v)\frac{dv}{q(x,v)}=\hat{\psi}(x)$

where $(\psi,q)_{q^{-1}}$ denotes the inner product of $\psi(x,v)$ and
$q(x,v)$ in $V$ space using $q(x,v)^{-1}$ as a weight function. Hence (38) is
solved as

$\psi(x,v)=-\frac{1}{\mu}\phi(x,v).$

Similarly, for $\mathcal{L}^{*}$ we consider
$\phi\in\mbox{ker}\,\mathcal{L^{*}}^{\perp}$ and find
$\psi\in\mbox{ker}\mathcal{L^{*}}^{\perp}$ that solves
$\mathcal{L^{*}}\psi=\phi$. We find

(39)
$-\mu\psi(x,v)+\mu\int_{V}q(x,v^{\prime})\psi(x,v^{\prime})dv^{\prime}=\phi(x,v).$

We integrate this equation over $V$ and use the fact that $\phi,\psi\in\langle
1\rangle^{\perp}$ to find

$0+\mu|V|\int_{V}q(x,v^{\prime})\psi(x,v^{\prime})dv^{\prime}=0.$

Hence equation (39) is solved as

$\psi(x,v)=-\frac{1}{\mu}\phi(x,v).$

implying that $\left(\mathcal{L}|_{\langle
q\rangle^{\perp}}\right)^{-1}=-1/\mu$.

We use these results in the next section where we derive expressions for the
mean first passage time through parabolic scaling in the case of anisotropic
transport where $K(x,v,v^{\prime})=q(x,v)$ and under symmetry conditions
$q(x,v)=q(x,-v)$.

###  5.1 Parabolic Scaling of the Forward Transport Equation

As mentioned above, parabolic scaling of the forward transport equation is
well established in the literature [7, 29]. In the anisotropic case the
forward transport equation (6) for $p(x,v,t)$ reads

(40) $p_{t}+v\cdot\nabla p=\mathcal{L}p.$

We now introduce a small parameter $\varepsilon>0$ and rescale macroscopic
time and space through the parabolic scaling variables $\tau$ and $\xi$
defined as

(41) $\tau=\varepsilon^{2}t,\qquad\xi=\varepsilon x.$

The forward transport equation now becomes

(42) $\varepsilon^{2}p_{\tau}+\varepsilon v\cdot\nabla_{\xi}p=\mathcal{L}p,$

where $p(\xi,v,\tau)$ is expressed in terms of the rescaled quantities. Note
that $\cal L$ contains no derivatives so the formal expression for ${\cal L}p$
in (42) is unaffected by the scaling procedure. We now expand (42) in powers
of $\varepsilon$

(43) $p(\xi,v,\tau)=p_{0}(\xi,v,\tau)+\varepsilon
p_{1}(\xi,v,\tau)+\varepsilon^{2}p_{2}(\xi,v,\tau)+\cdots$

where the functions $p_{0},p_{1},p_{2},\ldots$ are to be determined by
comparing orders of $\varepsilon$ in (42). Results are as follows

  * $\varepsilon^{0}:$

Collecting the zeroth-order in $\varepsilon$ terms in (42) leads to

$0=\mathcal{L}p_{0}.$

Hence, $p_{0}\in\mbox{ker}\,\mathcal{L}$ by Lemma 5.1 and

(44)
$p_{0}(\xi,v,\tau)=\hat{p}_{0}(\xi,\tau)q(\xi,v),\qquad\hat{p}_{0}(\xi,\tau)=\int_{V}p_{0}(\xi,v,\tau)dv.$

  * $\varepsilon^{1}:$

Collecting the first-order in $\varepsilon$ terms in (42) leads to

$v\cdot\nabla_{\xi}p_{0}=\mathcal{L}p_{1}.$

This equation can be solved if the left hand side is in the orthogonal
complement of the kernel of $\mathcal{L}$ in $L^{2}_{q^{-1}}(V)$. A direct
calculation yields

$(v\cdot\nabla_{\xi}p_{0},q)_{q^{-1}}\equiv\int_{V}v\cdot\nabla_{\xi}\hat{p}_{0}(\xi,\tau)q(\xi,v)\frac{dv}{q(\xi,v)}=\nabla_{\xi}\cdot\int_{V}vq(\xi,v)dv\
\hat{p}_{0}(\xi,\tau)=0,$

where $(v\cdot\nabla_{\xi}p_{0},q)_{q^{-1}}$ denotes the inner product of the
two arguments in $V$ space using $q^{-1}$ as a weight function. The last
identity comes from the symmetry condition $q(x,v)=q(x,-v)$, which implies
$q(\xi,v)=q(\xi,-v)$ and thus $\int_{V}vq(\xi,v)dv=0$. Then, according to
Lemma 5.3

(45)
$p_{1}(\xi,v,\tau)=-\frac{1}{\mu}v\cdot\nabla_{\xi}\left(\hat{p}_{0}(\xi,\tau)q(\xi,v)\right).$

  * $\varepsilon^{2}:$

Collecting the second-order in $\varepsilon$ terms in (42) leads to

$p_{0\tau}+v\cdot\nabla_{\xi}p_{1}=\mathcal{L}p_{2}.$

We integrate this equation over $V$ and use (44) and (45) to find

$\int_{V}\hat{p}_{0\tau}(\xi,\tau)q(\xi,v)\,dv+\int_{V}v\cdot\nabla_{\xi}\left(-\frac{1}{\mu}v\cdot\nabla_{\xi}(\hat{p}_{0}(\xi,\tau)q(\xi,v))\right)dv=0,$

which we can write using tensor notation as

$\hat{p}_{0\tau}=\frac{1}{\mu}\nabla_{\xi}\otimes\nabla_{\xi}:\left[\hat{p}_{0}(\tau,\xi)\int_{V}v\otimes
vq(\xi,v)dv\right].$

Finally, the parabolic limit of the forward transport equation can be written
through the introduction of a diffusion tensor $\mathbb{D}(\xi)$ as follows

(46)
$\hat{p}_{0\tau}(\xi,\tau)=\nabla_{\xi}\otimes\nabla_{\xi}:\left[\mathbb{D}(\xi)\hat{p}_{0}(\xi,\tau)\right],\quad\mbox{with}\quad\mathbb{D}(\xi)=\frac{1}{\mu}\int_{V}v\otimes
vq(\xi,v)dv.$

  * Equation (46) is a well studied anisotropic diffusion equation [58]. The mean first passage time $\tau(\xi)$ to the boundary of the spatial domain can be obtained by introducing the survival probability as shown in Section 2.5 and following the same procedures here described. We do this below.

###  5.2 Parabolic Scaling of the Backward Kolmogorov Equation

Under anisotropic, symmetric transport, the backward Kolmogorov equation for
the survival probability $s(x,v,t)$ in (34) can be written as

$s_{t}-v\cdot\nabla s=\mathcal{L}^{*}s.$

Similarly to the forward case, we rescale space and time as in (41) to obtain

(47) $\varepsilon^{2}s_{\tau}-\varepsilon
v\cdot\nabla_{\xi}s=\mathcal{L}^{*}s.$

We can now expand $s(\xi,v,\tau)$ in powers of $\varepsilon$ so that

$s(\xi,v,\tau)=s_{0}(\xi,v,\tau)+\varepsilon
s_{1}(\xi,v,\tau)+\varepsilon^{2}s_{2}(\xi,v,\tau)+\cdots$

and compare orders of $\varepsilon$ as follows

  * $\varepsilon^{0}:$

Collecting the zeroth-order in $\varepsilon$ terms in (47) leads to

$0=\mathcal{L}^{*}s_{0}.$

Hence, $s_{0}\in{\mbox{ker}}\,{\cal L}^{*}$ by Lemma 5.1 and is independent of
$v$:

(48) $s_{0}(\xi,v,\tau)=s_{0}(\xi,\tau).$

  * $\varepsilon^{1}:$

Collecting the first-order in $\varepsilon$ terms in (47) leads to

$-v\cdot\nabla_{\xi}s_{0}=\mathcal{L}^{*}s_{1}.$

We check the solvability condition on $\langle 1\rangle^{\perp}$:

$(v\cdot\nabla_{\xi}s_{0},1)=\nabla_{\xi}\cdot\int_{V}vdv\,s_{0}(\xi,\tau)=0.$

The last identity arises from the assumption of symmetry in the velocity space
$V$ so that if $v\in V$ then also $-v\in V$. Lemma 5.3 implies

(49) $s_{1}(\xi,v,\tau)=\frac{1}{\mu}v\cdot\nabla_{\xi}s_{0}(\xi,\tau).$

  * $\varepsilon^{2}:$

Collecting the second-order in $\varepsilon$ terms in (47) leads to

$s_{0\tau}-v\cdot\nabla_{\xi}s_{1}=\mathcal{L}^{*}s_{2}=-\mu
s_{2}+\mu\int_{V}q(\xi,v^{\prime})s_{2}(\xi,v^{\prime},\tau)dv^{\prime}.$

We multiply this equation by $q(\xi,v)$, integrate over $V$, and express
$s_{1}$ through (49) to find

$\displaystyle
s_{0\tau}-\int_{V}q(\xi,v)v\cdot\nabla_{\xi}\left(\frac{1}{\mu}v\cdot\nabla_{\xi}s_{0}(\xi,\tau)\right)dv$
$\displaystyle=$
$\displaystyle-\mu\int_{V}q(\xi,v)s_{2}(\xi,v,t)dv+\mu\int_{V}q(\xi,v)dv\int_{V}q(\xi,v^{\prime})s_{2}(\xi,v^{\prime}\tau)dv^{\prime}=0,$

where we use the condition $\int_{V}q(x,v)dv=1$ which also implies
$\int_{V}q(\xi,v)dv=1$. The above identity can be rewritten as

$s_{0\tau}=\frac{1}{\mu}\int_{V}q(\xi,v)v\otimes
vdv:\nabla_{\xi}\otimes\nabla_{\xi}s_{0}.$

Finally, utilizing the same diffusion tensor as in (46) we write

$s_{0\tau}=\mathbb{D}(\xi):\nabla_{\xi}\otimes\nabla_{\xi}s_{0},$

which is the parabolic limit of the backward equation (34). If we integrate
this equation over time from $0$ to $\infty$, and use the fact that
$s_{0}(\xi,0)=1$ then we obtain the MFPT equation for $\tau(\xi)$

(50) $-1=\mathbb{D}(\xi):\nabla_{\xi}\otimes\nabla_{\xi}\tau(\xi),$

Shifting back to the original spatio-temporal coordinates $(x,t)$, and using
$T(x)=\varepsilon^{2}\tau(\varepsilon x)$, we find an expression for the mean
first passage time $T(x)$ that is independent of the velocity $v$

(51)
$-1=\mathbb{D}(x):\nabla\otimes\nabla\,T(x),\qquad\mathbb{D}(x)=\frac{1}{\mu}\int_{V}v\otimes
v\,q(x,v)dv.$

##  6 Applications

To illustrate our results we consider several examples of anisotropic oriented
movement in two-dimensional regions. We begin with the bimodal von-Mises
distribution introduced in (10) with $n=2$ to express the orientation biases
of the random walker. Since in this case
$V=\sigma{\mathbb{S}^{n-1}}=\sigma\mathbb{S}^{1}$, and
$\omega=\sigma^{n-1}=\sigma$ from (5), we write

(52) $q(x,v)=\frac{\tilde{q}(x,\theta)}{\sigma}=\frac{1}{4\pi\sigma
I_{0}(k(x))}\left(e^{k(x)\gamma(x)\cdot\theta}+e^{-k(x)\gamma(x)\cdot\theta)}\right)=\frac{\cosh({k(x)\gamma(x)\cdot\theta})}{2\pi\sigma
I_{0}(k(x))},$

where $k(x)\geq 0$ is a measure of concentration and the unit vector
$\gamma(x)$ defines the preferred orientation at location $x\in\Omega$. The
bimodal von-Mises distribution includes both directions $\pm\gamma(x)$ to
allow for symmetric movement along the preferred orientation. The functions
$I_{j}(k(x))$ are modified Bessel functions of the first kind of order $j\geq
0$. Note that if $k(x)=0$ then $q(x,v)$ becomes the uniform distribution on a
circle $q(x,v)=1/(2\pi\sigma)$ and that if we let $k(x)\to\infty$ the
distribution displays two $\delta$-singularities, one at $\theta=\gamma(x)$
and one at $\theta=-\gamma(x)$.

In this scenario, where particles move at fixed speed $\sigma$ and assuming
the turning kernel biases motion along a preferred orientation $\gamma(x)$
that does not depend on incoming velocity, we assume that we are close to the
parabolic limit, $\Theta(x,v)=T(x)$, and use (51). The diffusion tensor
$\mathbb{D}(x)$ that appears in (51) can be calculated using (11) leading to
[30]

(53)
$\mathbb{D}(x)=\frac{\sigma^{2}}{2\mu}\left(1-\frac{I_{2}(k(x))}{I_{0}(k(x))}\right)\mathbb{I}+\frac{\sigma^{2}}{\mu}\frac{I_{2}(k(x))}{I_{0}(k(x))}\gamma(x)\gamma^{T}(x).$

In the following examples, we show solutions of (51) using various choices for
the domain $\Omega$ and the fields $\gamma(x)$ and $k(x)$ in (53).

###  6.1 MFPT for Anisotropic Transport on a circular domain

Bica et al [7] studied anisotropic transport on a two-dimensional disk where
the bimodal von-Mises distribution was used to describe biological particles
moving in radially symmetric environments. Applications include intra-cellular
transport along microtubules, cancer cell migration along collagen fibres, or
animals aggregating near watering holes. Here we similarly assume a two-
dimensional circular domain $\Omega$ with radius $R_{0}$, and employ planar
polar coordinates $(r,\varphi)$. We now need to specify $k(x)$, $\gamma(x)$.
For most of this section we let the parameter of concentration be uniform,
$k(x)=k\geq 0$, however we keep general $k(x)$ in our analytical derivations.
We also assume that the preferred orientation $\gamma(x)$ is either radial,
i.e. as if along the spokes of a bicycle wheel, or circular, along the
direction perpendicular to the radial direction, following concentric circles
inside the disk. If we now write the radial unit vector as
$\hat{x}=x/||x||=(\cos\varphi,\sin\varphi)$ and its orthogonal as
$\hat{x}^{\perp}=(-\sin\varphi,\cos\varphi)$, then $\gamma(x)$ can be
represented as

(54) $\gamma(x)=\left\\{\begin{array}[]{ll}\hat{x}&\mbox{for radial
orientation,}\\\\[5.0pt] \hat{x}^{\perp}&\mbox{for circular
orientation.}\end{array}\right.$

These two choices can be combined into a common notation by introducing the
anisotropy indicator $\alpha(x)$ as follows [7]

(55a)
$\mathbb{D}(x)=\frac{\sigma^{2}}{\mu}\left\\{\begin{array}[]{lcl}\frac{1}{2}(1-\alpha(x))\mathbb{I}+\alpha(x)\hat{x}\hat{x}^{T}&&x\neq
0\\\\[4.30554pt] \frac{1}{2}\mathbb{I}&&x=0.\end{array}\right.$ where
$\alpha(x)$ for general $k(x)$ is given by (55b)
$\alpha(x)=\frac{I_{2}(k(x))}{I_{0}(k(x))}\times\left\\{\begin{array}[]{rl}+1&\mbox{for
radial orientation,}\\\\[5.0pt] -1&\mbox{for circular
orientation.}\end{array}\right.$

If we set $k(x)=k$ then $\alpha(x)=\alpha$ is also uniform, and since
$I_{2}(k)\leq I_{0}(k)$ for all $k\geq 0$ it is also true that
$\alpha\in(-1,1)$. Note that the uniform distribution that arises from setting
$k(x)=0$ in (52) corresponds to $\alpha(x)=0$ due to (55b) and the properties
of the Bessel functions. The limit $k(x)\to\infty$, which leads to a
singularity along $\pm\gamma(x)$ as discussed above, results in $\alpha(x)\to
1$ for radial orientation and $\alpha(x)\to-1$ for circular orientation. This
implies that motion is singularly biased along the radial direction for
$\alpha\to 1^{-}$ and along the tangential direction for $\alpha\to-1^{+}$.
Finally, in (55a) $\mathbb{I}$ is the identity matrix and the tensor product
$\hat{x}\hat{x}^{T}$ is given by

$\hat{x}\hat{x}^{T}=\left(\begin{array}[]{cc}\cos^{2}\varphi&\sin\varphi\cos\varphi\\\
\sin\varphi\cos\varphi&\sin^{2}\varphi\end{array}\right).$

We can now solve the anisotropic mean first passage time equation (51). Given
the radial symmetry, we use the notation $T_{\alpha}(r)$ to indicate the MFPT
for the random walker starting at $r$ under the anisotropy indicator $\alpha$.
A direct calculation of (51) and (55) assuming $\alpha(x)=\alpha(r)$ is only
radially dependent, yields

(56)
$\frac{1+\alpha(r)}{2r}\,\Big{(}rT^{\prime}_{\alpha}(r)\Big{)}^{\prime}-\frac{\alpha(r)}{r}T^{\prime}_{\alpha}(r)=-\frac{1}{2D}$

where the diffusion constant $D=\sigma^{2}/(2\mu)$ is consistent with the
$n=2$ diffusivity introduced in (32). Note that the case $\alpha=0$
corresponds to the standard diffusion process in two-dimensions. The general
solution of (56) is

(57)
$T_{\alpha}(r)=-\frac{r^{2}}{4D}+H_{1}+H_{2}\int^{r}\exp\left[-\int^{\eta}\frac{1}{s}\frac{1-\alpha(s)}{1+\alpha(s)}ds\right]d\eta,$

where $H_{1},H_{2}$ are constants that depend on the chosen boundary
conditions. For uniform $\alpha(s)=\alpha$, (56) reduces to

(58)
$T_{\alpha}(r)=-\frac{r^{2}}{4D}+H_{1}+H_{2}r^{\frac{2\alpha}{1+\alpha}}.$

Below, we solve for $H_{1},H_{2}$ in given geometric scenarios and for fixed
$\alpha\in(-1,1)$. Note that since (58) contains two unknowns, two conditions
are necessary to identify them.

####  6.1.1 Exit from a disk:

In this scenario we assume the random walker can exit the disk only through
the boundary at $r=R_{0}$. This implies that the mean first exit time at
$r=R_{0}$ must be zero. We also impose smoothness at the origin so that all
derivatives at $r=0$ be zero. Thus, we solve (58) in the domain $0<r<R_{0}$
with boundary conditions

(59) $T^{\prime}_{\alpha}(0)=0,\quad T_{\alpha}(R_{0})=0,$

to obtain the solution

(60) $T_{\alpha}(r)=\frac{1}{4D}\big{(}R_{0}^{2}-r^{2}\big{)}.$

Remarkably, (60) is independent of $\alpha$. Thus, the two-dimensional MFPT to
reach the boundary of a disk for a particle moving under the von-Mises jump
distribution (52) with a uniform concentration parameter $k(x)=k$ has the same
form as in the case of a random, diffusive process.

(a) Exit from inner circle.

(b) Exit from inner or outer circle.

Figure 2: The MFPT (56) for $D=0.5$ on an annulus with domain $\rho<r<R_{0}$,
inner circle radius $\rho=0.5$, outer circle radius $R_{0}=3$ and several
values of $\alpha\in(-1,1)$. Panel (a): Particles exit the annulus only
through the inner circle; the corresponding boundary conditions under which to
solve (58) are $T_{\alpha}(\rho)=T^{\prime}_{\alpha}(R_{0})=0$. Note that
$T_{\alpha}(r)$ increases with $\alpha$ for given $r>\rho$ and
$T_{\alpha}(r>\rho)\to\infty$ as $\alpha\to-1^{+}$ as shown in (63e). Under
purely circular motion particles move tangentially and cannot decrease their
distance from the origin, hence they will never exit the annulus. Panel (b):
Particles may exit the annulus through both the inner and outer circles; the
corresponding boundary conditions under which to solve (58) are
$T_{\alpha}(\rho)=T_{\alpha}(R_{0})=0$. In the limit $\alpha\to-1^{+}$ the
solution $T_{\alpha}(r)$ displays a discontinuity at $r=\rho$ as particles
located at $r\to\rho^{+}$ will be able to exit the annulus only through the
outer circle at $r=R_{0}$. Particle simulations (black squares) based on
$N=10^{4}$ independent stochastic simulations of (1) with $\mu=10^{4}$,
$\sigma=\sqrt{\mu}$ agree with theoretical predictions. Specialized numerical
algorithms are needed to simulate the exit from the inner circle as
$\alpha\to{-1}^{+}$ since $T_{\alpha}\gg 1$. Indeed, we note the difference in
the scale of the vertical axes between the two panels. All quantities are in
arbitrary units.

####  6.1.2 Exit from the inner circle of an annulus:

We now consider motion constrained to an annular region so that $\rho<r<R_{0}$
for all trajectories. As shown below, the effects of anisotropy can be more
pronounced here compared to motion in a disk. We specifically assume that the
random walker can only exit the annulus via the inner circle at $r=\rho$ and
impose reflecting (Neumann) boundary conditions at the outer circle at
$r=R_{0}$. Thus we solve (58) in the domain $\rho<r<R_{0}$ with boundary
conditions

(61) $T_{\alpha}(\rho)=0,\quad T^{\prime}_{\alpha}(R_{0})=0.$

The general solution of (56) under uniform $\alpha$ is

(62) $T_{\alpha}(r)=\frac{1}{4D}(\rho^{2}-r^{2})+\frac{1}{2\beta
D}R_{0}^{\frac{\beta}{\alpha}}(r^{\beta}-\rho^{\beta}),\qquad{\mbox{with}}\qquad\beta\equiv\frac{2\alpha}{1+\alpha}.$

We then derive the limits for the MFPT corresponding to motion that is purely
radial ($\alpha\to 1^{-}$), isotropic ($\alpha\to 0$), or purely circular
($\alpha\to-1^{+}$)

(63a) $\displaystyle\lim_{\alpha\to{1^{-}}}T_{\alpha}(r)$
$\displaystyle=\frac{1}{4D}(\rho^{2}-r^{2})+\frac{R_{0}}{2D}(r-\rho),$ (63b)
$\displaystyle\lim_{\alpha\to 0}T_{\alpha}(r)$
$\displaystyle=\frac{1}{4D}(\rho^{2}-r^{2})+\frac{R_{0}^{2}}{2D}\log\Big{(}\frac{r}{\rho}\Big{)},$
(63e) $\displaystyle\lim_{\alpha\to-1^{+}}T_{\alpha}(r)$
$\displaystyle=\left\\{\begin{array}[]{cl}+\infty&r>\rho;\\\\[4.0pt]
0&r=\rho.\end{array}\right.$

Since the function $r-\rho-R_{0}\log(r/\rho)$ is a negative, decreasing
function of $r$ for $\rho<r<R_{0}$ the MFPT in the purely radial case
$(\alpha\to 1^{-}$ in (63a)) is less than in the isotropic case $(\alpha\to 0$
in (63b)). This implies that exiting the annulus from its inner circle is
faster if the motion is constrained along the radial direction rather than
being isotropic. Also note that under purely circular motion
($\alpha\to-1^{+}$ in (63e)) the MFPT $T_{\alpha}(r>\rho)=+\infty$ for any
starting position $r$. Since the preferred motion is along the tangent
direction, a particle that started at $r>\rho$ can never randomly turn to
smaller values of $r$ to reach the only exit boundary at $r=\rho$. Hence, it
will remain trapped in the annulus for infinite time.

####  6.1.3 Exit from the outer circle of an annulus:

We now consider the opposite scenario where the random walker can only exit
the annulus via the outer circle at $r=R_{0}$ and impose reflecting (Neumann)
boundary conditions at the inner circle at $r=\rho$. Thus we solve (58) in the
domain $\rho<r<R_{0}$ with boundary conditions

(64) $T^{\prime}_{\alpha}(\rho)=0,\quad T_{\alpha}(R_{0})=0.$

The general solution of (56) under uniform $\alpha$ is given by replacing
$\rho$ and $R_{0}$ in (62)

(65) $T_{\alpha}(r)=\frac{1}{4D}(R_{0}^{2}-r^{2})+\frac{1}{2\beta
D}\rho^{\frac{\beta}{\alpha}}(r^{\beta}-R_{0}^{\beta}),\qquad{\mbox{with}}\qquad\beta\equiv\frac{2\alpha}{1+\alpha}.$

The limits for the MFPT corresponding to motion that is purely radial
($\alpha\to 1^{-}$), isotropic ($\alpha\to 0$), or purely circular
($\alpha\to-1^{+}$) are

(66a) $\displaystyle\lim_{\alpha\to{1^{-}}}T_{\alpha}(r)$
$\displaystyle=\frac{1}{4D}(R_{0}^{2}-r^{2})+\frac{\rho}{2D}(r-R_{0}),$ (66b)
$\displaystyle\lim_{\alpha\to 0}T_{\alpha}(r)$
$\displaystyle=\frac{1}{4D}(R_{0}^{2}-r^{2})+\frac{\rho^{2}}{2D}\log\Big{(}\frac{r}{R_{0}}\Big{)},$
(66c) $\displaystyle\lim_{\alpha\to-1^{+}}T_{\alpha}(r)$
$\displaystyle=\frac{1}{4D}(R_{0}^{2}-r^{2}).$

In this case, the function $r-R_{0}-\rho\log(r/R_{0})$ is a negative,
increasing function of $r$ for $\rho<r<R_{0}$ so that the MFPT in the purely
radial case $(\alpha\to 1^{-}$ in (66a)) is still less than in the isotropic
case $(\alpha\to 0$ in (66b)). Thus, exiting the annulus from its outer circle
is faster if particle motion is constrained along the radial direction rather
than being isotropic. Since $r<R_{0}$, comparing (66c) to (66a) and (66b)
shows that purely circular motion yields the longest MFPT to exiting the
annulus. Furthermore for $\alpha\to-1^{+}$ the MFPT does not depend on $\rho$
since the motion is purely tangential and, as discussed above, a particle
starting at $r>\rho$ will always randomly increase its distance from origin.

Upon comparing (63a) and (66a), one can show that under purely radial motion,
the MFPT to exit the annulus from the outer or inner circle is the same if
trajectories start at $r=r_{\rm c}=(\rho+R_{0})/2$, the exact midpoint. For
$r>r_{\rm c}$ the MFPT is shortest if exiting from the outer circle, for
$r<r_{\rm c}$ the MFPT is shortest if exiting from the inner circle.

####  6.1.4 Exit from the inner or outer circle of an annulus:

Finally, we solve (58) in the case where a particle can exit the annulus from
both inner and outer circles using the following boundary conditions

$T_{\alpha}(\rho)=T_{\alpha}(R_{0})=0.$

We find

(67)
$T_{\alpha}(r)=-\frac{1}{4D}r^{2}+\frac{1}{4D}\left(\frac{R_{0}^{\beta}\rho^{2}-R_{0}^{2}\rho^{\beta}}{R_{0}^{\beta}-\rho^{\beta}}\right)+\frac{r^{\beta}}{4D}\left(\frac{R_{0}^{2}-\rho^{2}}{R_{0}^{\beta}-\rho^{\beta}}\right),\qquad\beta=\frac{2\alpha}{1+\alpha}.$

As done above, we highlight the limiting cases of purely radial ($\alpha\to
1^{-}$), isotropic ($\alpha=0$), and purely circular ($\alpha\to-1^{+}$)
motion to find

(68a) $\displaystyle\lim_{\alpha\to{1^{-}}}T_{\alpha}(r)$
$\displaystyle=\frac{1}{4D}(\rho-r)(r-R_{0}),$ (68b)
$\displaystyle\lim_{\alpha\to 0}T_{\alpha}(r)$
$\displaystyle=\frac{1}{4D}\frac{(R_{0}^{2}-\rho^{2})\log
r-(r^{2}-\rho^{2})\log R_{0}+(r^{2}-R_{0}^{2})\log\rho}{\log R_{0}-\log\rho},$
(68e) $\displaystyle\lim_{\alpha\to-1^{+}}T_{\alpha}(r)$
$\displaystyle=\left\\{\begin{array}[]{cl}\displaystyle{\frac{1}{4D}}(R_{0}^{2}-r^{2})&r>\rho;\\\\[4.0pt]
0&r=\rho.\end{array}\right.$

By comparing (68a) with (63a) and (66a), one can show that under purely radial
motion the MFPT to exit the annulus from either the inner or outer circle is
always less compared to when exit is allowed from only one of the two
boundaries. Note the discontinuity in (68e): under purely circular motion
particles starting at $r>\rho$ can exit the annulus through the outer circle
so we recover the same result as in (66c).

###  6.2 MFPT for anisotropic transport with linear features.

Figure 3: Top row: The MFPT to the boundary of the domain $\Omega=[-1,1]^{2}$
in the presence of one or multiple linear segments orienting particle motion.
Results are obtained by numerically solving (51) for $T(x)$ in the domain
$\Omega=[-1,1]^{2}$ under the bimodal von-Mises distribution (52) modulated by
the linear structures shown in the inset as described in (69). The cutoff
distance $d_{0}$ for particles to move along the linear structures is set at
$d_{0}=0.02$ and the corresponding measure of concentration is $k_{0}=25$. The
anisotropy profiles from left to right are a vertical line, a slant line, 3
parallel vertical lines and 10 randomly intersecting lines. Bottom row: The
difference between the MFPT under anisotropic transport and the MFPT for a
random walker in the same two-dimensional domain. The scenario of 10 randomly
intersecting segments in the right-most panels yields the smallest difference
between the two MFPTs. All quantities are in arbitrary units.

As discussed in (2.3), anisotropic transport is often facilitated by the
presence of linear structures in the environment. For example, roads and
seismic lines modulate animal movement in ecological landscapes while intra-
cellular transport occurs along actin filaments or microtubules. Multiple
segments (of roads, or filaments) carry random relative orientations and
occasionally intersect. As such, we study the MFPT in the presence of one or
several preferential orientations. For concreteness, we use a square domain
$\Omega$ and assume that particles can exit the domain through any of its
edges. We thus solve (51) with absorbing (Dirichlet) boundary conditions. To
model the presence of linear features we introduce a bimodal von-Mises
distribution constructed using (52) and a tessellation process. For each point
$x$ we find the Euclidean distance $d(x)$ to the closest linear feature and
evaluate the anisotropic distribution in (52) by setting $k(x),\alpha(x)$ as
follows

(69) $k(x)=\left\\{\begin{array}[]{rl}k_{0},&d(x)<d_{0}\\\\[4.0pt] \qquad
0,&\mbox{otherwise}\end{array}\right.,\qquad\alpha(x)=\frac{I_{2}(k(x))}{I_{0}(k(x))}.$

The directional field $\gamma(x)=(\cos\varphi,\sin\varphi)$ is taken as the
slope of the nearest line to point $x$. These choices imply that once a
particle is at a distance greater than a cutoff value $d_{0}$ from any linear
feature, i.e. for $d(x)\geq d_{0}$, movement is an isotropic random walk since
in this case $k(x)=\alpha(x)=0$. Once the particle is closer to a linear
feature, for $d(x)<d_{0}$, movement is biased towards a bidirectional random
walk along the $\pm\gamma(x)$ directions corresponding to the specific linear
feature as per (52). The larger the value of $k_{0}$, the larger the bias
along $\pm\gamma(x)$. Note that setting $\alpha(x)$ positive in (69)
corresponds to radial orientation in (55b) and that $k_{0}\to\infty$ yields
$\alpha(x:d(x)<d_{0})\to 1$. The contour plots shown in the top row of Fig. 3
are the numerical solutions of (51) and (69) using linear segments given by a
vertical line, a slant line, 3 parallel vertical lines and 10 randomly
intersecting lines as shown in the inset and using $k_{0}=25$ and $d_{0}=0.02$
which correspond to $\alpha(x:d(x)<d_{0})=0.922$. For these parameters the
effective motion is a combination of isotropic two-dimensional ($d(x)\geq
d_{0}$) and quasi one-dimensional random walks $(d(x)<d_{0})$. The modulation
of the MFPT due to the simple linear features is clearly distinguishable in
all four scenarios shown in Fig. 3. We also numerically evaluate the MFPT in
the same square geometry but under isotropic diffusion, which corresponds to
setting $d_{0}\to 0$ in (69), and show the difference between the two MFPTs in
the bottom row of Fig. 3. The right-most panels in Fig. 3 with 10 randomly
oriented linear features yield an MFPT that is closest to an isotropic random
walker.

###  6.3 MFPT on an ecological landscape

Finally, we consider an actual ecological landscape with linear features to
specifically model animal movement using the results and methods of Section
6.2. McKenzie et al. [41, 42] have previously used the anisotropic model (51)
and fitted it to oriented movement of red foxes in Minnesota in the presence
of prey and a den site [41], and to wolf movement data in forest areas of the
Canadian Rocky Mountains, where forests are disturbed by seismic lines [42].
They find that preferred movement of wolves along seismic lines increases the
encounter rate of predator and prey, thereby reducing prey density.

Here we revisit this case as an illustration and use an aerial image that was
not used in [42]. Fig. 4(a) shows a landscape of intersecting roads and
seismic lines in the boreal forest of Western Canada. In [30, 29] the forward
transport equation (6) was solved on the rectangular domain of this image.
Here we use the same image to evaluate the MFPT to the boundary of the domain.
The image was first digitized and thresholded such that landscape features are
represented by a boolean variable. The binarized image is then processed into
two fields, $\gamma(x)$ and $d(x)$, representing the direction of, and
distance to, the closest linear feature for each $x\in\Omega$, where $\Omega$
is the rectangular domain of the image. The resulting linear features and
their associated directions are shown in Fig. 4(b). The directional field is
incorporated into the bimodal Von-Mises distribution (52) so that motion near
linear features is biased along the two directions $\pm\gamma(x)$. In the
blank space of Fig. 4(b), where directional information is not shown, the
diffusion is isotropic.

(a) Forest image.

(b) Linear features with orientation.

(c) Normalized MFPT $T(x)/T_{\rm max}.$

Figure 4: Calculation of the MFPT on an ecological landscape. Panel (a):
Aerial snapshot of the boreal forest in Western Canada. Panel (b): A digitized
and thresholded version of the same image highlighting roads and seismic lines
with bidirectional information. Panel (c): Contour plot of the MFPT as derived
from (51), (53) and (69) using the method described in Section 6.2. Also shown
are three single trajectories from particle simulations of (1), each
initialized at the origin with parameters $\mu=10^{4}$, $\sigma=10^{4}$ and
bimodal von-Mises turning kernel (52). Particle trajectories remain closely
anchored to the linear segments before eventual absorption at different edges
of $\partial\Omega$. All quantities are in arbitrary units.

Similarly to what done in Section 6.2, the MFPT equation (51) is solved on the
image’s rectangular domain with $k(x)$ given in (69) and the direction field
shown in Fig. 4(b). We allow the random walker to exit the domain through any
of the edges of the image and apply absorbing (Dirichlet) boundary conditions
on $\partial\Omega$. The contour plot of the resulting numerical solution of
(51) is shown in Fig. 4(c). We observe a similar phenomenon to that observed
in Fig. 3 where the superposition of linear features with different angles
yields an apparently isotropic MFPT. Fig. 4(c) also displays three single
trajectories obtained from stochastic simulations of (1) which shows strong
alignment on the linear features. Hence, the occupation density of particles
is strongly localized on the anisotropic features, while the MFPT is modulated
to by these features to a much lesser degree. This, of course, is because the
MFPT is the average of the mean first passage time derived from all possible
trajectories.

##  7 Conclusions

Kinetic transport equations form a large class of models that are used in many
applications from physics, engineering, material sciences. Here we focus on
biological applications, which include a wide variety species and spatio-
temporal scales such as E. coli run and tumble movement [47], brain tumor
invasion [62], sea turtle orientation [52], DNA repair mechanisms [21], and
many more.

As it turns out, the general transport equation (1) is too general to derive a
meaningful mean first passage time equation. Hence, we focused on three cases:
the diffusive, anisotropic, and parabolic limit cases, all of which are
relevant to biology. Each of these cases translates into a set of assumptions
on the model parameters and in particular on the integral kernel
$K(x,v,v^{\prime})$. We summarize these assumptions in Table 1 for accessible
comparison. The diffusive and anisotropic cases are not disjoint, since, for
example, the constant kernel $K(x,v,v^{\prime})=|V|^{-1}$ is contained in
both.

Case | Reference | Assumptions
---|---|---
General assumption |  (3) | $K\geq 0,\quad K\in L^{2}(V\times V),\quad\int_{V}K(x,v,v^{\prime})dv=1.$
Diffusive case |  (3) + (23) |  $K(x,v,v^{\prime})=K(x,v,-v^{\prime}),\quad\int_{V}K(x,v,v^{\prime})dv^{\prime}=1$.
Anisotropic case |  (3) + (33) |  $K(x,v,v^{\prime})=q(x,v),\quad q(x,-v)=q(x,v)$.
Parabolic scaling |  (3) + (33) | $\tau=\varepsilon^{2}t,\quad\xi=\varepsilon x.$
|  \+ (41,43) | $p(\xi,v,\tau)=p_{0}(\xi,v,\tau)+\varepsilon p_{1}(\xi,v,\tau)+\varepsilon^{2}p_{2}(\xi,v,\tau)+\cdots$
Table 1: Summary of the basic assumptions for the various cases studied in
this manuscript.

In our analysis we find two types of MFPT equations. For the diffusive and
anisotropic cases we find an integro-PDE for the MFPT $\Theta(x,v)$, (31) and
(37), respectively. The expected exit time $\Theta(x,v)$ depends on the
initial location $x$ and the initial velocity $v$, and the explicit dependence
on $v$ makes this equation so interesting. A systematic study of equations of
the type (31) and (37) has not been done yet. In the parabolic scaling case,
the $v$-dependence disappears, and we obtain the anisotropic MFPT equation
(51) for $T(x)$. This model fits within the classical theory of MFPT
equations, as, for example, illustrated in [64, 58, 63, 56], and much is known
about its properties and its solutions.

Solutions to MFPT equations depend, of course, on the domain $\Omega$ and the
boundary conditions on $\partial\Omega$. For the parabolic limit case (51),
boundary conditions can be included in a standard way, using $T(x)=0$ on the
absorbing part of the boundary and $\hat{n}\cdot\mathbb{D}(x)\nabla T(x)=0$ on
the no-flux part of the boundary, where $\hat{n}$ denotes an outward normal
vector on the boundary. The domain can be connected or disconnected, it can
have holes and other special shapes, and the literature on various domains is
vast (for example [56, 24, 14, 54, 33]).

For the diffusive or the anisotropic case, however, boundary conditions need
to include the velocity $v$. A reasonable choice to describe absorbing
boundaries is to assume that a random walker on the boundary $\partial\Omega$
that starts in an outward direction, is lost immediately. In mathematical
terms, for any $x\in\partial\Omega$ such that $\hat{n}(x)$ denotes the unit
outward normal vector, we have

(70) $\Theta(x,v)=0\quad\mbox{for all}\quad v\in V:\;v\cdot\hat{n}>0.$

On the reflecting parts of the boundary various conditions can be stipulated.
They often take the form of an integral relation [59]. Given
$x\in\partial\Omega$, for each $v\in V$ with $v\cdot\hat{n}(x)\leq 0$, we
define

(71)
$\Theta(x,v)=\int_{\\{v^{\prime}\cdot\hat{n}>0\\}}I(x,v,v^{\prime})\Theta(x,v^{\prime})dv^{\prime},$

with an appropriate integral kernel $I$. Combining the MFPT equations (31) and
(37) with these type of boundary conditions (70) and (71) provides a
formidable mathematical and numerical challenge. It will require substantial
new ideas to develop theories for existence, uniqueness, boundedness etc. We
hope our work stimulates further research in this direction.

It should be noted that the MFPT equations derived here have no drift term,
i.e. no first order term. This is by design, since we assume symmetry of the
underlying kernels, $K(x,v,v^{\prime})=K(x,v,-v^{\prime})$ and
$q(x,v)=q(x,-v)$, see Table 1 above. These symmetry conditions ensure that
certain integrals in our derivation are zero. If these conditions are relaxed,
those integrals stay, and we expect them to lead to drift terms in the end
result. We did not work out the details here, since our focus was on
developing the overall framework to derive the MFPT equation for transport
equations. By relaxing the above assumptions, more general models can be
derived.

The bimodal von-Mises distribution is an excellent non-trivial example to
study the models developed here. It obeys the conditions listed above and
allowed us to analytically solve (51) for several movement types and
geometries, including radial and circular motion on a disk and on an annulus.
We also found numerical solutions for a specific preferred orientation (or
orientations) as provided by the underlying environment. Finally, we compared
the MFPTs in these various settings to the isotropic case and discussed under
what conditions anisotropy yields longer or shorter MFPTs to the boundary of
given domains. Our work provides a framework to quantify important timescales
associated to biological or animal movement processes without resolving the
entire forward process.

There are many avenues for future work arising from this study. Our examples
of MFPT calculation in sections 6.2 and 6.3 suggest that the presence of
multiple linear anisotropic features produces a superposition that is
independent of individual heterogeneities. Hence, we hypothesize that a
homogenization limit can be obtained to reflect an effective description of
diffusion processes in crowded domains. Our principal consideration of the
bimodal von-Mises produces unbiased one-dimensional diffusion along linear
features. It is natural to extend the formulations obtained herein towards
more general classes of distributions including single modal von-Mises and
multi- modal distributions such as Kent and Bingham distributions in two-
dimensions or higher. These may be more difficult to evaluate as, for example,
the symmetry condition $q(x,v)=q(x,-v)$ does not hold for the single modal
von-Mises distribution in (8). Other avenues of exploration include deeper
investigations of ecological scenarios, for example, by studying how the MFPT
changes upon allowing animals to exit only from select portions of the
boundary, modifying the parameter $k_{0}$ or allowing it to have different
values along different segments for the purpose of optimally steering herds of
animals to given locations, where, for example, foraging conditions are more
suitable.

## Acknowledgments

We acknowledge fruitful discussions with Kevin Painter and Stuart Johnson. TH
is supported through a discovery grant of the Natural Science and Engineering
Research Council of Canada (NSERC), RGPIN-2023-04269. AEL acknowledges support
from NSF grant DMS-2052636. MRD acknowledges support from ARO grant
W911NF-23-1-0129.

##  Appendix A Discretization of the anisotropic MFPT equation

Our numerical simulation of (37) with diffusivity tensor (55a) is based on
finite difference simulation in Matlab. We write the following for
$x=(x_{1},x_{2})$ and $\nabla=(\partial_{x_{1}},\partial_{x_{2}})$

(72) $\mathbb{D}(x):\nabla\otimes\nabla
T=D_{11}\partial^{2}_{x_{1}x_{1}}T+2D_{12}\partial^{2}_{x_{1}x_{2}}T+D_{22}\partial^{2}_{x_{2}x_{2}}T=-1,\qquad\mathbb{D}(x)=\begin{pmatrix}D_{11}&D_{12}\\\
D_{12}&D_{22}\end{pmatrix}.$

In the formulation above, the tensor coefficients $D_{11}$, $D_{12}$ and
$D_{22}$ are spatially varying functions defined in (55). A mesh for the
rectangular domain $\Omega=(a,b)\times(c,d)$ is formed at points
$x_{1}^{(j)}=a+(j-1)\Delta x_{1}$ for $j=1,\ldots,N_{1}$ and
$x_{2}^{(k)}=c+(k-1)\Delta x_{2}$ with $k=1,\ldots,N_{2}$, such that

$\frac{b-a}{N_{1}-1}=\Delta x_{1},\qquad\frac{d-c}{N_{2}-1}=\Delta x_{2}.$

We define $T^{(j,k)}=T(x_{1}^{(j)},x_{2}^{(k)})$ and utilize the standard
finite difference approximations

(73a)
$\displaystyle\partial^{2}_{x_{1}x_{1}}T^{(j,k)}=\frac{T^{(j+1,k)}-2T^{(j,k)}+T^{(j-1,k)}}{\Delta
x_{1}^{2}},\qquad\partial^{2}_{x_{2}x_{2}}T^{(j,k)}=\frac{T^{(j,k+1)}-2T^{(j,k)}+T^{(j,k-1)}}{\Delta
x_{2}^{2}},$ (73b)
$\displaystyle\partial^{2}_{x_{1}x_{2}}T^{(j,k)}=\frac{T^{(j+1,k+1)}-T^{(j+1,k-1)}-T^{(j-1,k+1)}+T^{(j-1,k-1)}}{4\Delta
x_{1}\Delta x_{2}}.$

Hence, combining the discretizations (73) with the equation (72), we arrive at
the discrete linear system

$\displaystyle D_{11}^{(j,k)}\frac{T^{(j+1,k)}-2T^{(j,k)}+T^{(j-1,k)}}{\Delta
x_{1}^{2}}+D_{12}^{(j,k)}\frac{T^{(j+1,k+1)}-T^{(j+1,k-1)}-T^{(j-1,k+1)}+T^{(j-1,k-1)}}{2\Delta
x_{1}\Delta x_{2}}$ $\displaystyle+\
D_{22}^{(j,k)}\frac{T^{(j,k+1)}-2T^{(j,k)}+T^{(j,k-1)}}{\Delta x_{2}^{2}}=-1.$

In our simulation, we apply the Dirichlet boundary so that $T^{(j,k)}=0$ for
any $j\in\\{1,N_{1}\\}$ and any $k\in\\{1,N_{2}\\}$.

## References

  * [1] LJS. Allen, An Introduction to Stochastic Processes with Applications to Biology, Second Edition, Chapman & Hall/CRC Mathematical and Computational Biology, Boca Raton, FL, 2010.
  * [2] W. Alt, Biased random walk model for chemotaxis and related diffusion approximation, J. Math. Biol., 9 (1980), pp. 147–177.
  * [3] D. Ambrosi, A. Gamba, and G. Serini, Cell directional and chemotaxis in vascular morphogenesis, Bull. Math. Biol., 66 (2004), pp. 1851–1873.
  * [4] L. Angelani, R. Di Leonardo, and M. Paoluzzi, First-passage time of run-and-tumble particles, Eur. Phys. J. E, 37 (2014), pp. 1–6.
  * [5] AJ. Bernoff, A. Jilkine, A. Navarro Hernández, and AE. Lindsay, Single-cell directional sensing from just a few receptor binding events, Biophys. J., 122 (2023), pp. 3108–3116.
  * [6] AJ. Bernoff and AE. Lindsay, Numerical approximation of diffusive capture rates by planar and spherical surfaces with absorbing pores, SIAM J. Appl. Math., 78 (2018), pp. 266–290.
  * [7] I. Bica, T. Hillen, and KJ. Painter, Aggregation of biological particles under radial directional guidance, J. Theor. Biol., 427 (2017), pp. 77–89.
  * [8] PC. Bressloff and JM. Newby, Stochastic models of intracellular transport, Rev. Mod. Phys., 85 (2013), pp. 135–196.
  * [9] PC. Bressloff and RD. Schumm, The narrow capture problem with partially absorbing targets and stochastic resetting, Multiscale Model. Simul., 20 (2022), pp. 857–881.
  * [10] JA. Carrillo, M. Fornasier, G. Toscani, and F. Vecil, Particle, kinetic, and hydrodynamic models of swarming, Springer, 2010, pp. 297–336.
  * [11] C. Cercignani, The Boltzmann equation, in The Boltzmann equation and its applications, Springer, 1988, pp. 40–103.
  * [12] C. Cercignani, R. Illner, and M. Pulvirenti, The mathematical theory of diluted gases, Springer, New York, 1994.
  * [13] MAJ. Chaplain, M. Lachowicz, Z. Szymańska, and D. Wrzosek, Mathematical modelling of cancer invasion: The importance of cell–cell adhesion and cell–matrix adhesion, Math. Models Methods in Appl. Sci., 21 (2011), pp. 719–743.
  * [14] T. Chou and MR. D’Orsogna, First passage problems in biology, in First-passage phenomena and their application, World Scientific, Singapore, 2014, pp. 306–345.
  * [15] M. Conte, Y. Dzierma, S. Knobe, and C. Surulescu, Mathematical modeling of glioma invasion and therapy approaches via kinetic theory of active particles, Math. Models Methods in Appl. Sci., 33 (2023), pp. 1009–1051.
  * [16] MR. D’Orsogna and T. Chou, First passage and cooperativity of queuing kinetics, Phys. Rev. Lett., 95 (2005), p. 170603.
  * [17] C. Engwer, T. Hillen, M. Knappitsch, and C. Surulescu, Glioma follow white matter tracts: A multiscale DTI-based model, J. Math. Biol., 71 (2015), pp. 551–582.
  * [18] M. Faran and G. Bisker, Nonequilibrium self-assembly time forecasting by the stochastic landscape method, J. Phys. Chem. B, 127 (2023), pp. 6113–6124.
  * [19] P. Fauchald and T. Torkild, Using first-passage time in the analysis of area-restricted search and habitat selection, Ecology, 84 (2003), pp. 282–288.
  * [20] PW. Fok and T. Chou, Accelerated search kinetics mediated by redox reactions of DNA repair enzymes, Biophys. J., 96 (2009), pp. 3949––3958.
  * [21] PW. Fok, CL. Guo, and T. Chou, Charge-transport-mediated recruitment of DNA repair enzymes, J. Chem. Phys., 129 (2008), p. 235101.
  * [22] KR. Ghusinga, JJ. Dennehy, and A. Singh, First-passage time approach to controlling noise in the timing of intracellular events, Proc. Natl. Acad. Sci., 114 (2017), pp. 693–698.
  * [23] MM. Gibbons, T. Chou, and MR. D’Orsogna, Diffusion-dependent mechanisms of receptor engagement and viral entry, J. Phys. Chem. B, 114 (2010), pp. 15403–15412.
  * [24] DS. Grebenkov, Universal formula for the mean first passage time in planar domains, Phys. Rev. Lett., 117 (2016), p. 260201.
  * [25] T. Hillen, On the ${L}^{2}$-moment closure of transport equations: The general case, Discrete & Contin. Dyn. Syst. -B, 5 (2005), pp. 299–318.
  * [26] , ${M}^{5}$ mesoscopic and macroscopic models for mesenchymal motion, J. Math. Biol., 53 (2006), pp. 585–616.
  * [27] , Elements of Applied Functional Analysis, AMS Open Notes, Providence, 2023.
  * [28] T. Hillen and H. Othmer, The diffusion limit of transport equations derived from velocity jump processes, SIAM J. Appl. Math., 61 (2000), pp. 751–775.
  * [29] T. Hillen and KJ. Painter, Transport and anisotropic diffusion models for movement in oriented habitats, in Dispersal, individual movement and spatial ecology, Springer, 2013, pp. 177–222.
  * [30] T. Hillen, KJ. Painter, AC. Swan, and A. Murtha, Moments of von-Mises and Fisher distributions and applications, Math. Biosci. Eng., 14 (2017), pp. 673–694.
  * [31] D. Holcman and Z. Schuss, The narrow escape problem, SIAM Review, 56 (2014), pp. 213–257.
  * [32] P. Kumar, J. Li, and C. Surulescu, Multiscale modeling of glioma pseudopalisades: Contributions from the tumor microenvironment, J. Math. Biol., 82 (2021), pp. 1–45.
  * [33] V. Kurella, JC. Tzou, D. Coombs, and MJ. Ward, Asymptotic analysis of first passage time problems inspired by ecology, Bull. Math. Biol., 77 (2015), pp. 83–125.
  * [34] T. Lagache, C. Sieben, T. Meyer, A. Herrmann, and D. Holcman, Stochastic model of acidification, activation of hemagglutinin and escape of influenza viruses from an endosome, Front. Phys., 5 (2017), p. 25.
  * [35] SD. Lawley, AE. Lindsay, and CE. Miles, Receptor organization determines the limits of single-cell source location detection, Phys. Rev. Lett., 125 (2020), p. 018102.
  * [36] M. Le Corre, C. Dussault, and SD. Côté, Detecting changes in the annual movements of terrestrial migratory species: Using the first-passage time to document the spring migration of caribou, Mov. Ecol., 19 (2014), pp. 1–11.
  * [37] M. Le Corre, M. Pellerin, D. Pinaud, G. Van Laere, H. Fritz, and S. Saïd, A multi-patch use of the habitat: Testing the first-passage time analysis on roe deer Capreolus capreolus paths, Wildl. Biol., 14 (2008), pp. 339–349.
  * [38] AE. Lindsay, AJ. Bernoff, and A. Navarro Hernández, Short-time diffusive fluxes over membrane receptors yields the direction of a signalling source, R. Soc. Open Sci., 10 (2023), p. 221619.
  * [39] AE. Lindsay, AJ. Bernoff, and MJ. Ward, First passage statistics for the capture of a brownian particle by a structured spherical target with multiple surface traps, Multiscale Model. Simul., (2017), pp. 74–109.
  * [40] AE. Lindsay, T. Kolokolnikov, and JC. Tzou, Narrow escape problem with a mixed trap and the effect of orientation, Phys. Rev. E, 91 (2015), p. 032111.
  * [41] HW. McKenzie, MA. Lewis, and EH. Merrill, First passage time analysis of animal movement and insights into the functional response, Bull. Math. Biol., 71 (2009), pp. 107–129.
  * [42] HW. McKenzie, EH. Merrill, RJ. Spiteri, and MA. Lewis, How linear features alter predator movement and the functional response, Interface Focus, 2 (2012), pp. 205–216.
  * [43] J. Morgan and AE. Lindsay, Modulation of antigen discrimination by duration of immune contacts in a kinetic proofreading model of T cell activation with extreme statistics, PLOS Comp. Biol., 19 (2023), pp. 1–17.
  * [44] NM. Mutothya and X. Yong, Mean first passage time for diffuse and rest search in a confined spherical domain, Physica A Stat. Mech. Appl., 567 (2021), p. 125667.
  * [45] HG. Othmer, S.R. Dunbar, and W. Alt, Models of dispersal in biological systems, J. Math. Biol., 26 (1988), pp. 263–298.
  * [46] HG. Othmer and T. Hillen, The diffusion limit of transport equations ii: Chemotaxis equations, SIAM J. Appl. Math., 62 (2002), pp. 1222–1250.
  * [47] HG. Othmer and C. Xue, Multiscale models of taxis-driven patterning in bacterial populations, SIAM Appl. Math., 70 (2009), pp. 133–167.
  * [48] A. Padash, AV. Chechkin, B. Dybiec, M. Magdziarz, B. Shokri, and R. Metzler, First passage time moments of asymmetric Lévy flights, J. Phys. A Math. Theor., 53 (2020), p. 275002.
  * [49] KJ. Painter, Modelling migration strategies in the extracellular matrix, J. Math. Biol., 58 (2009), pp. 511–543.
  * [50] KJ. Painter, Multiscale models for movement in oriented environments and their application to hilltopping in butterflies, Theor. Ecol., 7 (2014), pp. 53–75.
  * [51] KJ. Painter and T. Hillen, Mathematical modelling of glioma growth: The use of Diffusion Tensor Imaging (DTI) data to predict the anisotropic pathways of cancer invasion, J. Theor. Biol., 323 (2013), pp. 25–39.
  * [52] KJ. Painter and T. Hillen, Navigating the flow: Individual and continuum models for homing in flowing environments, J. R. Soc. Interface, 12 (2015), p. 20150647.
  * [53] B. Perthame, Transport Equations in Biology, Frontiers in Mathematics, Birkhäuser Verlag, Basel, Switzerland, 2007.
  * [54] S. Pillay, MJ. Ward, A. Peirce, and T. Kolokolnikov, An asymptotic analysis of the mean first passage time for narrow escape problems: Part I: Two-dimensional domains, Multiscale Model. Simul., 8 (2010), pp. 803–835.
  * [55] CE. Plunkett and SD. Lawley, Boundary homogenization for patchy surfaces trapping patchy particles, J. Chem. Phys., 158 (2023), p. 094104.
  * [56] S. Redner, A guide to first-passage processes, Cambridge University Press, Cambridge, UK, 2001.
  * [57] K. Rijal, A. Prasad, and D. Das, Protein hourglass: Exact first passage time distributions for protein thresholds, Phys. Rev. E, 102 (2020), p. 052413.
  * [58] H. Risken, The Fokker-Planck equation: Methods of solution and applications, Springer-Verlag, Berlin, DE, 2nd ed., 1992.
  * [59] H. Schwetlick, Limit sets for multidimensional nonlinear transport equations, J. Differ. Equations, 179 (2002), pp. 356–368.
  * [60] TL. Stepien, C. Zmurchok, JB. Hengenius, RM. Caja Rivera, MR. D’Orsogna, and AE. Lindsay, Moth mating: Modeling female pheromone calling and male navigational strategies to optimize reproductive success, Appl. Sci., 10 (2020), p. 6543.
  * [61] DW. Stroock, Some stochastic processes which arise from a model of the motion of a bacterium, Prob. Theory Rel., 28 (1974), pp. 305–315.
  * [62] AC. Swan, T. Hillen, J. Bowman, and A. Murtha, An anisotropic model for glioma spread, Bull. Math. Biol., 80 (2017), pp. 1259––1291.
  * [63] NG. Van Kampen, Stochastic processes in physics and chemistry, North Holland, Amsterdam, NL, 3rd ed., 2007.
  * [64] G. Weiss, First passage time problems in chemical physics, Adv. Chem. Phys., 13 (1967), pp. 1–18.
  * [65] R. Yvinec, MR. D’Orsogna, and T. Chou, First passage times in homogeneous nucleation and self-assembly, J. Chem. Phys., 137 (2012), p. 244107.