2.1 Model formulation

Let $p = p(t,x,w)$ denote the cell density distribution, defined at time $t\geq0$ , position $x \in \Omega \subseteq \mathbb{R}^d$ and velocity $w \in V$ . We typically assume V to be a compact set as in [Reference Hillen and Othmer26], and in particular consider $V=[s_1,s_2]\times\mathbb{S}^{d-1} $ ([Reference Hillen24]), where $\mathbb{S}^{d-1}$ is the set of all possible directions $\hat{w} \in \mathbb{R}^d$ and the hat-symbol indicates unit vectors. The limits $s_1$ and $s_2$ denote the minimal and maximal speed ( $|w|$ ) of the cells, with $0\leq s_1 \le s_2 <\infty$ . Note that if the speed is approximately constant, we simply set $V=s\mathbb{S}^{d-1}$ .

The governing transport equation for describing cell movement is

(2.1) \begin{equation}\frac{\partial p(t,x,w)}{\partial t} + w\cdot \nabla p(t,x,w) = \mathcal{L} p(t,x,w) \, ,\end{equation}

where the operator $\nabla$ denotes the spatial gradient. The turning operator $\mathcal{L}$ is a linear operator that models the change in velocity of individuals per unit of time at (x,w) that is not due to the free particle drift. $\mathcal{L}$ is generally defined as an integral operator on $L^2$ spaces [Reference Hillen and Othmer26]:

\begin{equation*}\begin{array}{lr}\displaystyle \mathcal{L} : \displaystyle L^2(V) \longmapsto L^2(V) \, , \\[10pt]\displaystyle p(t,x,w)\longmapsto \mathcal{L}p(t,x,w) \, ,\end{array}\end{equation*}

where (t,x) are independent parameters and

(2.2) \begin{equation}\mathcal{L}p(t,x,w)=-\mu(t,x,w) p(t,x,w)+\displaystyle \int_{V} \mu(t,x,w') q(t,x,w,w') p(t,x,w') dw' \, .\end{equation}

This operator describes the velocity scattering. As noted earlier, the key determinants of the migration path are the turning rate function, $\mu(t,x,w)$ , and the turning kernel, q(t,x,w,w’). Viewed in this light, the first term on the right-hand side of (2.2) models particles switching away from velocity w and the second one takes into account the particles switching into velocity w from all other velocities.

The turning kernel, q(t,x,w,w’), denotes the probability measure of switching velocity from w’ to w, given that a turn occurs at location x and time t. Here we adopt the same reasoning as in [Reference Hillen24], by assuming that reorientation is dominated by the fibrous/anisotropic environmental structure, such as collagen matrix fibres or white matter tracts. The choice of new direction is therefore derived directly from this structure (rather than, say, incoming velocity) and, for simplicity, we assume that q is directly proportional to the distribution of the oriented environmental fibres. As such, we interchangeably use ‘turning kernel’ or ‘fibre distribution’ when referring to q throughout this manuscript. Based on the above, for q we assume

The simple assumption adopted in [Reference Hillen24] was to directly link a probability measure describing the directional distribution of fibres, $\tilde{q}(t,x,\hat{w})$ , defined on $\mathbb{R}_+\times \Omega\times \mathbb{S}^{d-1}$ and satisfying $\tilde{q}\geq 0, \displaystyle \int_{\mathbb{S}^{d-1}}\tilde{q}d\hat{w}=1$ , to the turning kernel:

\begin{align*}q(t,x,w)=\dfrac{\tilde{q}(t,x,\hat{w})}{\omega}, \qquad \omega=\int_V \tilde{q}(t,x,\hat{w}) dw.\end{align*}

Note that q(t,x,w) assumes $w\in V$ while $\tilde q(t,x,\hat w)$ is defined for unit vectors only, but their only difference lies in a constant scaling factor that accounts for the difference between V and $\mathbb{S}^{d-1}$ . Given this, we will interchangeably call q the turning kernel or the distribution of fibre orientations. As a consequence of A1– A3, the operator (2.2) simplifies to

(2.3) \begin{equation}\mathcal{L}p(t,x,w)=-\mu(t,x,w) p(t,x,w)+q(t,x,w)\displaystyle \int_{V} \mu(t,x,w') p(t,x,w') \, dw' \, .\end{equation}

The turning rate function, $\mu(t,x,w)$ , gives the rate at which velocity switches are made for a particle located at x at time t, moving in direction w. It is at this point where we substantially diverge from [Reference Hillen24], lifting the assumptions on $\mu$ and allowing w-, x-, and t-dependence. Significantly, this allows the turning rate $\mu$ to depend directly on the fibre orientation q, for example, allowing a cell to continue movement with the same velocity if it is moving in the direction of highly aligned fibres. Note that as $\mu=\mu(t,x,w)$ is a turning rate, $1/\mu(t,x,w)$ defines the mean time spent by a cell running along a linear tract with velocity w between two consecutive turns performed at time t, location x. We assume

2.2 Statistical properties

To analyse (2.1) under (2.3), we make use of a number of statistical properties of the corresponding fibre and turning distributions, such as expectations and variances. A summary of these expressions is given in Table 1.

1. (1) Distribution of new directions. We consider the distribution of newly chosen directions, q, with expectation:

   (2.4) \begin{equation}{\bf E}_{q}(t,x)=\displaystyle \int_{V} q(t,x,w)w \,dw.\end{equation}
   This is also the the mean new velocity after a turn and has the variance–covariance matrix:
   (2.5) \begin{equation}\mathbb{V}_q(t,x) = \displaystyle \int_{V} q(t,x,w)\,(w-{\bf E}_q)\otimes(w-{\bf E}_q)\, dw\,.\end{equation}

2. (2) Cell mean velocity and variance. We introduce similar macroscopic quantities for the cell population, although we stress p is not itself a probability measure. First, we define the macroscopic density of the population p at time t and position x as:

   (2.6) \begin{equation}\bar{p}(t,x)=\displaystyle \int_{V}p(t,x,w) dw \,,\end{equation}
   and the total mass of the population in $\Omega$ ,
   \begin{align*}m(t)=\int_{\Omega}\bar{p}(t,x)dx\,.\end{align*}
   Note that, with no population kinetics and assuming suitably lossless boundary conditions, the total mass will be conserved in time. With these definitions in place, we can introduce the moments of the normalised cell distribution $\hat p(t,x,w)=\dfrac{p(t,x,w)}{\bar p(t,x)}$ , which will be a probability measure for all t and x. In particular, we can introduce the expectation:
   \begin{align*}{\bf E}_{\hat{p}}(t,x)=\int_V \hat{p}(t,x,w)w dw\,,\end{align*}
   which is the mean velocity of the normalised population, and the variance (variance–covariance matrix):
   \begin{align*}\mathbb{V}_{\hat{p}}(t,x) = \displaystyle \int_{V} \hat{p}(t,x,w)\,(w-{\bf E}_{\hat{p}})\otimes(w-{\bf E}_{\hat{p}})\, dw\,.\end{align*}
   The latter provides information on the width of the distribution $\hat{p}$ in different directions. This tensor is symmetric, but can be anisotropic, that is, the level sets of $\hat w\mapsto \hat w^T \mathbb{V}_{\hat p} \hat w$ are ellipsoids. With this we can identify the mean velocity of the cell population as:
   \begin{align*}\displaystyle \int_{V} p(t,x,w) w dw =\bar{p}(t,x){\bf E}_{\hat{p}}(t,x)\end{align*}
   and the variance–covariance of the population velocity as:
   \begin{align*}\displaystyle \int_{V} p(t,x,w)\,(w-{\bf E}_{p})\otimes(w-{\bf E}_{p})\, dw = \bar{p}(t,x)\mathbb{V}_{\hat{p}}(t,x).\end{align*}

3. (3) Turning part of the population. The turning operator definition (2.3) reveals a new macroscopic quantity:

   (2.7) \begin{equation}p_\mu (t,x,w) = \mu(t,x,w) p(t,x,w).\end{equation}
   $p_\mu$ can be interpreted as the part of the cell population that moves in direction w and is currently turning. Then, the total turning population per unit time is
   (2.8) \begin{equation}\bar{p}_{\mu}(t,x)=\displaystyle \int_V \mu(t,x,w) p(t,x,w) dw,\end{equation}
   which is the expression in (2.3). The turning operator (2.3) can then be re-written as:
   (2.9) \begin{equation}\mathcal{L}p(t,x,w)=\bar{p}_{\mu}(t,x) q(x,w) -\mu(t,x,w)p(t,x,w).\end{equation}
   By normalising $p_\mu$ , we can define the mean incoming velocity of the turning population as:
   \begin{align*}{\bf E}_{p_{\mu}}(t,x)=\dfrac{\displaystyle \int_V \mu(t,x,w) p(t,x,w) w dw}{\bar{p}_{\mu}(t,x) },\end{align*}
   and its variance–covariance matrix accordingly.

4. (4) Run times along the fibres. We discover later that the stationary distributions are proportional to the ratio $q/\mu$ . In fact, the corresponding distribution arises in many forthcoming calculations. Hence, we introduce

   \begin{align*} C(t,x) = \int_V\frac{q(t,x,w)}{\mu(t,x,w)} dw \end{align*}
   and the normalised distribution:
   (2.10) \begin{equation} T(t,x,w) = \frac{1}{C(t,x)} \frac{q(t,x,w)}{\mu(t,x,w)}. \end{equation}
   Since $\mu(t,x,w)$ is a turning rate, the function:
   \begin{align*} \tau(t,x,w)=\frac{1}{\mu(t,x,w)} \end{align*}
   is the mean time spent moving in direction w. Then, T(t,x,w) is the distribution of run times along the fibres of the network and its expectation:
   (2.11) \begin{equation}{\bf E}_{T}(t,x)=\displaystyle \dfrac{1}{C(t,x)} \int_{V} \dfrac{q(t,x,w)}{\mu(t,x,w)}w \,dw \,,\end{equation}
   is the average velocity along the fibre distribution. Note that this quantity discriminates between the choice of a parabolic scaling (leading to a diffusive limit, for ${\bf E}_T\approx 0$ ) or a hyperbolic scaling (for ${\bf E}_T\neq 0$ ). With this interpretation, we can regard the following normalisation constant:
   \begin{align*} C(t,x) = \int_V \tau(t,x,w) q(t,x,w) dw \end{align*}
   as the mean run time between turns. We can further define the displacement vector and the mean displacement vector, respectively
   (2.12) \begin{equation} \chi(t,x,w) = w\tau(t,x,w)\,, \qquad \bar \chi(t,x)= \int_V \chi(t,x,w) q(t,x,w) dw\,, \end{equation}
   such that
   \begin{align*} {\bf E}_{T}(t,x) = \frac{\bar \chi(t,x)}{C(t,x)} \end{align*}
   becomes the ratio of the mean displacement vector over the mean run time on the fibre network.

5. (5) Persistence. Persistence, $\psi_d$ , is a measure of a random walker’s tendency to maintain direction during directional changes. It has values $\psi_d\in [-1,1]$ , where $\psi_d=1$ denotes perfect persistence (continuing with the previous direction), $\psi_d=0$ denotes uniform turning and $\psi_d=-1$ indicates a switch into the opposite direction [Reference Othmer, Dunbar and Alt34,Reference Othmer and Hillen35]. Persistence is often computed as the mean cosine along a particle trajectory; however, in our abstract framework, we define it as the mean velocity of the equilibrium distribution T, that is,

   (2.13) \begin{equation}\psi_d(w)=\frac{{\bf E}_T \cdot w}{|{\bf E}_T| |w|}\,.\end{equation}
   Explicit use of the vector product shows that $\psi_d$ does indeed arise as a mean cosine, as:
   \begin{align*} \psi_d = \cos(\sphericalangle({\bf E}_T, w)),\end{align*}
   where $\sphericalangle({\bf E}_T, w)$ denotes the angle between ${\bf E}_T$ and w. We have not yet shown that T is the equilibrium distribution, but do so in the next section. Observe that if $\mu$ does not depend on w, equation (2.13) becomes
   (2.14) \begin{equation}\psi_d(w)=\frac{{\bf E}_q \cdot w}{|{\bf E}_q| |w|},\end{equation}
   as defined in [Reference Othmer, Dunbar and Alt34]. Hence, for turning rates that do not depend on the direction, the persistence will vanish for a bidirectional tissue (i.e. ${\bf E}_q=0$ ). Our extension to w-dependence in $\mu$ lifts this limitation, allowing non-vanishing persistence even for a bidirectional tissue with ${\bf E}_T\neq 0$ . In other words, equation (2.1) with (2.3) is capable of generating a persistent random walk even under a fully symmetric configuration.

Table 1.

Summary of the key probability distributions used during our analysis

We return to the turning operator (2.9). As expected, due to assumption A2, we observe that the total cell density during turning will be conserved:

\begin{eqnarray*}\int_{V} \mathcal{L}p(t,x,w) dw &=& \int_{V} \bar p_{\mu}(t,x) q(t,x,w) dw - \int_{V} \mu(t,x,w) p(t,x,w) dw\\&=& \bar p_\mu(t,x) - \bar p_{\mu}(t,x) \\&=& 0.\end{eqnarray*}

The average outgoing velocity of the total turning population, meanwhile, changes

(2.15) \begin{eqnarray}\int_{V} \mathcal{L}p(t,x,w) w \, dw&=& \int_V\bar p_\mu(t,x) q(t,x,w) w \, dw - \int_{V} \mu(t,x,w) p(t,x,w) w \, dw \nonumber \\&=& \bar{p}_{\mu}(t,x){\bf E}_q(t,x) - \bar{p}_{\mu}(t,x){\bf E}_{p_{\mu}}(t,x) \nonumber\\&=& \bar{p}_{\mu}(t,x)\big( {\bf E}_q(t,x)-{\bf E}_{p_{\mu}}(t,x)\big).\end{eqnarray}

Therefore, the mean velocity of the turning population, ${\bf E}_{p_{\mu}}$ , relaxes towards the average post-turning velocity, ${\bf E}_q$ , imposed by the fibre network.

For much of the analysis, we assume fibre orientations and turning rates are time-independent, that is, q(x,w) and $\mu(x,w)$ . Similar arguments apply for the time-dependent case, yet the notation becomes clumsier.

2.3 Equilibrium state of the transport equation

As the first step towards finding equilibrium distributions, we compute the kernel of the turning operator $\mathcal{L}$ . A function $\phi(w)$ belongs to $\mbox{ker}(\mathcal{L})$ if and only if it satisfies

\begin{align*}\begin{array}{lr}\mathcal{L}\phi(w)=0\\[8pt]\Longleftrightarrow -\mu(x,w) \phi(w)+q(x,w)\displaystyle \int_{V} \mu(x,w') \phi(w') dw'=0\\[8pt]\Longleftrightarrow \mu(x,w) \phi(w)= q(x,w)\displaystyle \int_{V} \mu(x,w') \phi(w') dw'.\end{array}\end{align*}

We can write

\begin{align*}\phi(w)=\frac{q(x,w)}{\mu(x,w')} \bar{\phi}_{\mu}(x)\,,\end{align*}

where $\bar{\phi}_{\mu}$ is a w-independent function:

\begin{align*}\bar{\phi}_{\mu}(x)=\displaystyle \int_{V} \mu(x,w) \phi(w) dw,\end{align*}

which is the fraction of turning cells of the population $\bar{\phi}$ . The w-dependence of elements of $\mbox{ker}(\mathcal{L})$ is given by the ratio $q(x,w)/\mu(x,w)$ , that is, it is given by the run-time distribution along the fibres, T(x,w), which we introduced in (2.10). Then

\begin{align*} \phi(x,w) = \bar\phi(x) T(x,w), \qquad \mbox{and} \qquad \mbox{ker}(\mathcal{L}) = \langle T (x,w) \rangle. \end{align*}

Hence, a stationary state of the equations (2.1)–(2.3) has the form:

(2.16) \begin{equation}M(x,w)= \beta(x) T(x,w)\end{equation}

which we call the Maxwellian. We remark that (2.1) with (2.3) is the linear Boltzmann equation. In particular, the quantity:

\begin{align*}\sigma(t,x,w,w')=q(t,x,w)\mu(t,x,w')\end{align*}

is the cross section and the equilibrium probability density (normalised to 1) given by $\mathcal{L}(T)=0$ is (2.10). The function T will be non-negative, since q and $\mu$ are non-negative, and it is an $L^1(V)$ function via assumption ${\bf M2}$ . Therefore, the stochastic process ruled by $\sigma$ satisfies micro-reversibility, that is,

\begin{align*}\sigma(x,w,w')T(x,w')=\sigma(x,w',w)T(x,w).\end{align*}

Consequently, existence and uniqueness of a non-negative solution $p\in L^1(\Omega\times V)$ of (2.1),(2.3) with initial condition $p^0 \in L^1(\Omega \times V)$ and non-absorbing boundary conditions [Reference Cercignani6] is a classical result of kinetic theory (see, e.g. [Reference Petterson40]). We consider, in particular, a special class of non-absorbing boundary conditions, given by no-flux boundary conditions [Reference Plaza42].

Arguing as in [Reference Bisi, Carrillo and Lods3,Reference Pettersson41], we can prove a linear version of the classical H-Theorem for the linear Boltzmann equations (2.1), (2.3) with initial condition $p^0=p(0,x,w) \in \, L^1(\Omega \times V)$ . Let us introduce the entropy $S=-H$ where, for any given convex function $\Phi:\mathbb{R}_+\rightarrow \mathbb{R}_+$ :

\begin{align*}\ H_{\Phi}[p|M](t)= \int_\Omega \int_V p(t,x,w) \Phi\left( \dfrac{p(t,x,w)}{M(x,w)}\right) \, dw \, dx , \qquad p(t,\cdot, \cdot) \in L^1(\Omega \times V)\,.\end{align*}

In the above, M is the Maxwellian satisfying

\begin{align*}\displaystyle \int_{\Omega \times V} M(x,w) \, dx dw = \displaystyle \int_{\Omega} \bar{p}\,^0(x) dx\end{align*}

where we assume that $p^0$ has finite mass $\displaystyle \int_{\Omega} \bar{p}\,^0 \, dx$ and entropy $H_{\Phi}[p^0|M](0)$ . Hence, $\beta(x)$ is such that $\displaystyle \int_{\Omega} \beta(x) dx =\displaystyle \int_{\Omega} \bar{p}^0(x) \, dx $ and, in particular, if a stationary state $p^{\infty}(x,w)$ exists, then $\beta(x)=\bar{p}^{\infty}(x)$ , so that

\begin{align*}p^{\infty}(x,w)=\bar{p}^{\infty}(x) T(x,w).\end{align*}

Under non-absorbing boundary conditions, via exactly the same procedure followed in [Reference Pettersson41] it is possible to prove that

\begin{align*}\frac{dH_{\Phi}[p|M]}{dt}(t) \leq 0 \quad \mbox{for} \quad t\ge 0 \qquad\end{align*}

and

\begin{align*}\frac{dH_{\Phi}[p|M]}{dt} (t) =0, \, \qquad \mbox{iff} \qquad p(t,x,w)=M(x,w),\end{align*}

where p(t,x,w) is the unique $L^1$ solution to (2.1),(2.3) with initial condition $p^0$ and non-absorbing boundary conditions. Furthermore, continuing to argue as in [Reference Pettersson41], it can be proved that provided $\displaystyle \int_{\Omega}\int_{V}\big( 1+w^2 +|\log p_0|\big) \, dw dx < \infty$ , then

\begin{align*}\lim_{t\rightarrow \infty} \displaystyle \int_{\Omega} \int_V |p(t,x,w)-M(x,w)| \, dw \, dx =0.\end{align*}

3.1 Diffusion-dominated case

In the diffusion-dominated case, we assume that the macroscopic drift term ${\bf E}_T=0$ . We consider this case first to introduce the necessary technical notations.

Theorem 3.1 Let assumptions A1– A3 be satisfied and assume

(3.3) \begin{equation}{\bf E}_T=0.\end{equation}

Consider equation (3.1) with expansion (3.2). The leading order term, $p_0(\tau,X,w)$ , of equation (3.2) satisfies

\begin{align*}p_0(\tau,X,w)=\bar p_0(\tau,X) \, T(X,w),\end{align*}

where $\bar p_0(\tau,X)$ solves the macroscopic anisotropic diffusion equation:

(3.4) \begin{equation}\frac{\partial}{\partial \tau} \bar p_0(\tau,X) = \nabla \cdot \int_V \frac{1}{\mu(X,w)}\nabla\cdot\Bigl[\bar p_0(\tau,X) \mathbb{D}_{T}(X,w)\Bigr] dw\,, \end{equation}

with microscopic anisotropic diffusion tensor:

\begin{align*}\mathbb{D}_{T} (X,w)= w \otimes w\, T(X,w).\end{align*}

Proof. We have seen that

\begin{align*} \mbox{ker}(\mathcal{L}) = \langle T\rangle. \end{align*}

In order to invert $\mathcal{L}$ on the orthogonal complement of its kernel, we use a specific weight for the inner product in $L^2_\eta(V)$ , with

(3.5) \begin{equation}\eta(X,w) = \frac{\mu(X,w)}{T(X,w)},\end{equation}

that is, given $f,g \in L^2$ the weighted scalar product is defined by:

\begin{align*}(f,g)_{\eta}=\int_V f(w)g(w)\eta(X,w) \, dw.\end{align*}

We observe that $\phi(X,w)\in \langle T \rangle^\perp$ if and only if

\begin{align*} 0= \int_V \phi(X,w) T(X,w) \frac{\mu(X,w)}{T(X,w)}dw = \int_V \phi(X,w) \mu(X,w) dw = \bar\phi_\mu(X). \end{align*}

The range of $\mathcal{L}$ is given by all functions that integrate to zero, since

\begin{align*} \int_V \mathcal{L} \phi(X,w) dw =\overline{\mathcal{L}\phi} =0\end{align*}

from mass conservation. Then, this range, as subset of $L^2_\eta(V)$ , can be written as:

\begin{align*} \mbox{Range}(\mathcal{L}) = \langle \eta^{-1}\rangle^\perp,\end{align*}

since

\begin{align*} 0= \int_V \phi(X,w) dw = \int_V\phi(X,w) \eta^{-1}(X,w) (\eta(X,w) dw) = (\phi, \eta^{-1})_\eta.\end{align*}

Then

\begin{align*} \mathcal{L}^\perp: \langle T\rangle^\perp \to \langle \eta^{-1}\rangle^\perp, \qquad \mbox{with} \quad \mathcal{L}^\perp =\mathcal{L}|_{\langle T\rangle^\perp}\end{align*}

and this restricted operator $\mathcal{L}^\perp$ is the operator we would like to invert. Given $\psi\in\langle \eta^{-1}\rangle^\perp$ , we need to find $\phi \in \langle T\rangle^\perp$ such that $\mathcal{L}\phi = \psi$ . Applying $\mathcal{L}$ from (2.9), we obtain (ignoring arguments for clarity of presentation)

\begin{align*}\mathcal{L} \phi = \bar \phi_\mu q - \mu \phi = \psi.\end{align*}

Since on $\langle T\rangle^\perp$ we have $\bar \phi_\mu=0$ , we can solve for $\phi$ as:

\begin{align*}\phi(X,w) = -\frac{1}{\mu(X,w)} \psi(X,w) .\end{align*}

Hence, we can write

(3.6) \begin{equation}\left(\mathcal{L}^\perp \right)^{-1}:\langle\eta^{-1}\rangle^\perp \to \langle T \rangle^\perp, \quad \psi\mapsto -\frac{1}{\mu} \psi.\end{equation}

We may observe that the pseudo-inverse depends on the microscopic velocity and on the macroscopic space coordinate through $\mu(X,w)$ .

We now substitute expansion (3.2) into equation (3.1) and match orders of $\varepsilon$ .

• • For $\varepsilon^{0}$ , we find

  \begin{align*}\mathcal{L}p_{0}(\tau,X,w)=0,\end{align*}
  hence $p_0\in \mbox{ker}(\mathcal{L})$ and we can write
  (3.7) \begin{equation}p_0(\tau, X, w) = \bar{p}_0(\tau, X)\; T(X, w) .\end{equation}

• • At order $\varepsilon^{1}$ , we have

  (3.8) \begin{equation}\nabla\cdot\big( w\, p_{0}(\tau,X, w)\big)= \mathcal{L} p_1(\tau,X,w)\,.\end{equation}
  To solve this equation for the next order correction term $p_1$ , we need to invert $\mathcal{L}$ . We saw earlier in (3.6) that $\mathcal{L}$ is invertible on $\langle \eta^{-1}\rangle^\perp $ . Hence, we check the solvability condition:
  \begin{align*} \nabla\cdot(wp_0(\tau,X,w)) \in \langle\eta^{-1} \rangle^\perp.\end{align*}
  Using (3.7), this condition becomes
  \begin{eqnarray*}0 &=& (\nabla\cdot(wp_0) , \eta^{-1})_\eta \\&=& \int_V \nabla\cdot(w\bar p_0(\tau,X) T(X, w)) \eta(X,w)^{-1} \eta(X,w) dw \\&=& \nabla \cdot \left[ \int_V w T(X,w) dw \,\bar p_0(\tau, X)\right] \\&=& \nabla\cdot \left[ \mathbb{E}_T(X) \bar p_0(\tau,X)\right]\,.\end{eqnarray*}
  Hence, in this step, it is necessary to assume
  (3.9) \begin{equation}\mathbb{E}_T(X) = 0.\end{equation}
  In words, we assume the distribution of run times along the fibre network has no dominant direction. In this case, we can invert (3.8) and find
  (3.10) \begin{equation}p_1(\tau,X, w) = \frac{-1}{\mu(X, w)} \nabla\cdot(w\, \bar p_0(\tau,X) \, T(X,w))\,.\end{equation}

• • In $\varepsilon^{2}$ :

  \begin{align*}\dfrac{\partial }{\partial \tau}p_{0}(\tau,X, w)+\nabla\cdot\big(p_{1}(\tau,X, w)w \big)=\mathcal{L} p_2(\tau,X,w)\,.\end{align*}
  Integrating this equation over V, we obtain
  \begin{align*} \int_V \frac{\partial }{\partial\tau}p_0(\tau, X, w) dw +\nabla\cdot \int_V w p_1(\tau,X, w) dw =0.\end{align*}
  Using the expressions (3.7) for $p_0$ and (3.8) for $p_1$ , we obtain
  \begin{eqnarray*}\dfrac{\partial }{\partial \tau}\bar p_0(\tau,X) \underbrace{\int_V T(X,w) dw}_{=1} &=& \nabla \cdot \int_V \frac{1}{\mu(X,w)} w\otimes w \nabla \Bigl[\bar p_0(\tau, X) T(X, w) \Bigr] dw \\ &=& \nabla \cdot\int_V \frac{1}{\mu(X,w)}\nabla\cdot\Bigl[\bar p_0(\tau,X) w\otimes w T(X, w)\Bigr] dw \\ &=& \nabla \cdot \int_V \frac{1}{\mu(X,w)}\nabla\cdot\Bigl[\bar p_0(\tau,X) \mathbb{D}_{T}(X,w)\Bigr] dw\,, \end{eqnarray*}
  with an anisotropic diffusion tensor:
  \begin{align*}\mathbb{D}_{T} (X,w)= w \otimes w\, T(X,w)\,.\end{align*}
  $\mathbb{D}_T$ is a microscopic scale diffusion tensor, where
  \begin{align*}\mathbb{V}_T(X) = \int_V \mathbb{D}_T(X,w) dw \end{align*}
  describes the variance of the run-time distribution along the network fibres.

Lemma 1 The above parabolic limit equation (3.4) can be written as an anisotropic drift-diffusion model:

(3.11) \begin{equation}\dfrac{\partial \bar{p}_0}{\partial \tau}(\tau,X) +\nabla \cdot \big(\mathbf{a}(X)\bar{p}(\tau,X)\big)= \nabla \cdot \nabla \cdot \big(\mathbb{D}(X)\bar{p_0}(\tau,X)\big) ,\end{equation}

with macroscopic diffusion tensor, $\mathbb{D}$ , and advection speed, $\mathbf{a}$ , given by:

(3.12) \begin{eqnarray}\mathbb{D}(X) &=& \int_V w\otimes w\, \frac{T(X,w)}{\mu(X,w)} dw \end{eqnarray}

(3.13) \begin{eqnarray}\mathbf{a}(X) &=& - \int_V w\otimes w\, \frac{\nabla\mu(X,w)}{\mu(X,w)}\frac{T(X,w)}{\mu(X,w)} dw.\end{eqnarray}

Proof. The proof relies on straightforward application of the quotient rule. Omitting arguments for readability,

\begin{align*} \nabla\cdot\nabla\cdot \left(\bar p_0 \int_V w\otimes w \frac{T}{\mu} dw\right) = \nabla \cdot \int_V \frac{1}{\mu} \nabla\cdot (\bar p_0 w\otimes w T) dw - \nabla\cdot \left(\bar p_0 \int_V \frac{\nabla \mu}{\mu^2} w\otimes w T dw \right). \end{align*}

This representation allows some interesting physicobiological interpretations.

• • In (3.5), earlier we defined the weight function $\eta=\frac{\mu}{T}$ . We can write the above macroscopic quantities in terms of $\eta$ as:

  \begin{eqnarray*}\mathbb{D}(X) &=& \int_V w\otimes w\, \frac{dw}{\eta(X,w)}\,, \\\mathbf{a}(X) &=& - \int_V w\otimes w\, \nabla\ln(\mu(X,w))\frac{dw}{\eta(X,w)}\,.\end{eqnarray*}
  $\mathbb{D}$ then appears as an anisotropy matrix related to the measure $\eta^{-1}dw $ , while $\mathbf{a}$ measures the anisotropic logarithmic gradient with the same measure $\eta^{-1} dw$ .

• • We can also relate the terms back to the original network structure given by q(X,w). Recall that $T=\frac{q}{C\mu}$ , where C(X) was a normalisation constant. By considering the distance travelled in direction w, we can write

  \begin{eqnarray*}\mathbb{D}(X) &=&\frac{1}{C(X)} \int_V \chi(X,w) \otimes \chi(X,w) q(X,w)\, dw\,, \\\mathbf{a}(X) &=& - \frac{1}{C(X)} \int_V \chi(X,w) \otimes \chi(X,w)\, \nabla\ln(\mu(X,w)) \, q(X,w) dw\,,\end{eqnarray*}
  using the definition of $\chi$ from (2.12). $\mathbb{D}$ is then the scaled variance–covariance matrix of the directed mean run time along the fibre network and $\mathbf{a}$ is the scaled mean logarithmic derivative of the turning rate, weighted by the mean run times along the fibre network.

• • The equations simplify drastically when the turning rate $\mu$ does not depend on space. In this case, $\nabla\mu =0$ and hence $\mathbf{a}=0$ , that is, we get a pure anisotropic diffusion equation. Viewed this way, the advection velocity $\mathbf{a}$ clearly results from spatial variation in the turning rate, $\mu(X,w)$ , and relates to an anisotropic taxis term measuring the drift due to the gradient of $\mu$ . Spatial dependence in $\mu(X,w)$ generates an advection pushing cells towards decreasing values of the turning frequency.

Due to the fact that V and $\Omega$ are compact, we can directly apply a convergence result of [Reference Degond, Goudon and Poupaud13], giving us

Lemma 2 Suppose A1– A3 and M1– M2 hold and assume (3.3). Let us also suppose that $\exists \, C_1,C_2 \geqslant0$ such that

\begin{align*}|w\cdot\nabla q(x,w)| \leqslant C_1\mu(x,w), \quad a.e. \quad \textit{in} \quad \Omega\times V\end{align*}

and

\begin{align*}\displaystyle\int_V w^2 \dfrac{q(x,w)}{\mu(x,w)^2} dw \le C_2 \quad a.e. \quad in \quad \Omega.\end{align*}

Let the initial condition $p^{0}_{\epsilon}(x,w)$ satisfy $\displaystyle \int_{\Omega \times V} \dfrac{p_{{\epsilon}^0}(x,w)^2}{T(x,w)}\, dx dw < \infty$ and $\bar{p}_{\epsilon}^0 \rightharpoonup \bar{p}^0$ in $H^{-1}(\Omega)$ . Let $p_{\epsilon}$ be the solution to (3.1). Then there exists a subsequence $\bar{p}_{\epsilon} \rightharpoonup \bar{p}$ in $L^2((0,T)\times \Omega)$ where $\bar{p}$ satisfies (3.11) with initial condition $\bar{p}_0(x)$ .

3.2 Drift-diffusion case

The authors of [Reference Goudon and Mellet22] study the parabolic scaling also in the case where the macroscopic drift ${\bf E}_T \ne 0$ . The key lies in a transformation to moving spatial coordinates, shifting the solution in the direction of ${\bf E}_T$ as $Z=X-{\bf E}_T t$ . This method also works here and we introduce

\begin{align*} p(\tau,X,w) := u(\tau, X-{\bf E}_T \tau, w). \end{align*}

Then, the transport equation (2.1) transforms to

\begin{align*} \frac{\partial}{\partial \tau}u +\nabla\cdot [ (w-{\bf E}_T) u] = \mathcal{L} u.\end{align*}

Therefore, the scaled transport equation (3.1):

\begin{align*}\varepsilon^2 \frac{d}{d\tau} p +\varepsilon \nabla \cdot (w p) = \mathcal{L} p \end{align*}

transforms to

\begin{align*} \varepsilon^2 \frac{\partial}{\partial \tau}u +\varepsilon \nabla\cdot [ (w-{\bf E}_T) u] = \mathcal{L} u.\end{align*}

For this modified transport equation, we perform the same scaling analysis as before. We consider

\begin{align*} u(\tau,Z,w) = u_0(\tau, Z, w) +\varepsilon u_1(\tau,Z, w) +\varepsilon^2 u_2(\tau, Z, w) + \cdots. \end{align*}

Upon comparing orders of $\varepsilon$ we find, to leading order, that

\begin{align*} \mathcal{L} u_0 = 0.\end{align*}

Hence, $u_0$ is in the kernel of $\mathcal{L}$ and we can write

\begin{align*} u_0(\tau,Z, w) = \bar u_0(\tau,Z) \; T(X, w). \end{align*}

The order $\varepsilon$ terms are

\begin{align*} \nabla\cdot\big[(w-{\bf E}_T) \bar u_0 T\big] = \mathcal{L} u_1.\end{align*}

To solve this equation for $u_1$ , we require the solvability condition:

\begin{align*}\nabla\cdot \left[\int_V (w-{\bf E}_T) T dw \; \bar u_0 \right] =0\,,\end{align*}

which is true since

\begin{align*} \int_V(w-{\bf E}_T) T dw = {\bf E}_T - {\bf E}_T =0.\end{align*}

Then

\begin{align*} u_1 = -\frac{1}{\mu} \nabla\cdot\big[(w-{\bf E}_T) \bar u_0 T\big].\end{align*}

The terms of order $\varepsilon^2$ are

\begin{align*}\frac{\partial}{\partial \tau} (\bar u_0 T) + \nabla\cdot \big[(w-{\bf E}_T) u_1\big] = \mathcal{L} u_2.\end{align*}

We integrate this equation over V to obtain

(3.14) \begin{equation}0= \frac{\partial}{\partial \tau}\bar u_{0} - \nabla\cdot \underbrace{\int_V (w-{\bf E}_T) \left(\frac{1}{\mu}\right) \nabla\cdot \big[(w-{\bf E}_T) \bar u_0 T\big] dw}_{(I)}.\end{equation}

Generalising the previous definitions of the diffusion tensor (3.12) and the drift velocity (3.13), we define now a macroscopic diffusion tensor $\mathbb{D}_h$ and macroscopic drift velocity $\mathbf{a}_h$ as:

(3.15) \begin{eqnarray}\mathbb{D}_h(X) &=& \int_V (w-{\bf E}_T)\otimes (w-{\bf E}_T)\, \frac{T}{\mu} dw\,, \end{eqnarray}

(3.16) \begin{eqnarray}\mathbf{a}_h(X) &=& \int_V (w-{\bf E}_T) \nabla \cdot\left[\frac{w-{\bf E}_T}{\mu}\right] T\; dw\,.\end{eqnarray}

Note that for ${\bf E}_T=0$ , we return to the previous definitions (3.12) and (3.13) of the diffusion-dominated case.

Therefore, we have that

(3.17) \begin{eqnarray}\nabla \cdot (\mathbb{D}_h \bar u_0) &=& \int_V (-\nabla \cdot {\bf E}_T) (w-{\bf E}_T) \frac{T}{\mu} dw\; \bar u_0 +\underbrace{\int_V(w-{\bf E}_T) \frac{1}{\mu} \nabla \cdot\big[(w-{\bf E}_T) \bar u_0 T\big] \, dw}_{(I)}\nonumber \\&&+ \int_V (w-{\bf E}_T) \otimes (w-{\bf E}_T) T \nabla\left(\frac{1}{\mu}\right) dw\; \bar u_0\nonumber\\&=& (I) + \int_V(w-{\bf E}_T) \left[\nabla\cdot(w-{\bf E}_T) \frac{1}{\mu} + (w-{\bf E}_T) \cdot\nabla\left(\frac{1}{\mu}\right) \right] Tdw \;\bar u_0 \nonumber \\&=& (I) + {\bf a}_h \bar u_0\,. \end{eqnarray}

With these definitions, we obtain from (3.14) a fully anisotropic drift-diffusion model for $\bar u_0$ :

(3.18) \begin{equation} \frac{\partial}{\partial \tau}\bar u_{0} +\nabla\cdot ({\bf a}_h \bar u_0) = \nabla\cdot\nabla\cdot (\mathbb{D}_h \bar u_0) \,.\end{equation}

Finally, transforming back to $\bar p_0 (\tau, X) = \bar u_0(\tau, X-{\bf E}_T \tau)$ we find

(3.19) \begin{equation} \frac{\partial}{\partial \tau}\bar p_{0} +\nabla\cdot (({\bf a}_h+{\bf E}_T) \bar p_0) = \nabla\cdot\nabla\cdot (\mathbb{D}_h \bar p_0) \,.\end{equation}

If ${\bf E}_T=0$ , we get the same result as already shown in Theorem 3.1 and equation (3.11) in Lemma 1.

3.3 Hyperbolic limit

In drift-dominated phenomena, we expect that non-dimensionalisation leads to a scaling in which the macroscopic time and space scales are of the form $\tau=\varepsilon t, X=\varepsilon x$ , $\varepsilon \ll 1$ . Effectively, cells do not have a large turning frequency, drift dominates and we can perform a hyperbolic limit. The rescaled equation is

(3.20) \begin{equation}\varepsilon \dfrac{\partial p}{\partial \tau} + \varepsilon w\cdot \nabla p =\mathcal{L} p \,.\end{equation}

We consider the following expansion:

(3.21) \begin{equation}p= p_0 + \varepsilon g +\mathcal{O}(\varepsilon^2),\end{equation}

where we assume that the correction term g carries zero mass ( $\bar g =0$ ) and, therefore, $\bar p_0=\bar p$ . Substituting (3.21) into ( $3.20$ ), we find that the leading order terms ( $\varepsilon^{0}$ ) are $\mathcal{L} p_0= 0$ , hence

(3.22) \begin{equation}p_0(\tau,X, w) = \bar p(\tau,X)\, T(X,w)\,.\end{equation}

With this choice of $p_0$ , the remaining terms in (3.20) are

(3.23) \begin{equation}\dfrac{\partial p}{\partial \tau} + \varepsilon g_\tau + w\cdot \nabla p + \varepsilon w\cdot \nabla g = \mathcal{L} g + O(\varepsilon^2).\end{equation}

We integrate this equation over V, use the form of $p_0$ from above (3.22), and assume the $O(\varepsilon^2)$ -terms are negligible. Then

\begin{align*} \dfrac{\partial \bar p}{\partial \tau} +\varepsilon\int_V g_\tau dw + \nabla\cdot \int_V w \bar p T dw + \varepsilon \nabla\cdot \int_V w gdw = 0 .\end{align*}

Since g carries no mass, we have $\frac{\partial}{\partial \tau} \bar g =0$ . Using the definition of the expectation of T from (2.11), we obtain

(3.24) \begin{equation} \dfrac{\partial \bar p}{\partial \tau}+\nabla \cdot ( {\bf E}_T \bar{p}) + \varepsilon \nabla\cdot \int_{V} w g dw =0 \,,\end{equation}

which, to leading order, becomes the pure drift model:

(3.25) \begin{equation}\dfrac{\partial \bar p}{\partial \tau}(\tau, X) +\nabla\cdot \Bigl({\bf E}_T(X) \bar p(\tau,X)\Bigr) =0 .\end{equation}

The macroscopic drift velocity ${\bf E}_T$ is the expected movement direction based on the average time spent on the fibres. If ${\bf E}_T$ is of order 1, then (3.25) is the leading order model. However, for small or zero expectation ${\bf E}_T$ , we can compute the next-order correction term g. Specifically, we use $\displaystyle \int_V wg \, dw= \displaystyle \int_V (w-{\bf E}_T)g \, dw$ and rewrite equation (3.24) as:

(3.26) \begin{equation}\dfrac{\partial \bar{p}}{\partial \tau}+\nabla \cdot \big( {\bf E}_T\bar p \big)\\= -\varepsilon \nabla \cdot\Bigl[\int_V(w-\mathbf{E}_T) g dw\Bigr].\end{equation}

From the previous equation (3.23), we find to leading order that

(3.27) \begin{equation}\mathcal{L} g = \dfrac{\partial p}{\partial \tau} +w\cdot \nabla p = \dfrac{\partial \bar{p}}{\partial \tau} T + w\cdot \nabla (\bar p T) = -\nabla\cdot({\bf E}_T\bar p) T + w\cdot \nabla(\bar p T),\end{equation}

where we used the drift model (3.25) in the last step. To solve for g, we need to satisfy the solvability condition:

(3.28) \begin{eqnarray}0 &=& \int_V -\nabla \cdot (\mathbf{E}_T \bar p) T + w\cdot \nabla (\bar p T) dw \nonumber\\[3pt]&=& -\nabla\cdot({\bf E}_T \bar p) +\nabla \cdot\left(\int_V w T dw \bar p\right)\nonumber \\[3pt]&=& -\nabla\cdot({\bf E}_T \bar p) +\nabla\cdot({\bf E}_T \bar p) \,,\end{eqnarray}

which is true since $\int T dw = 1$ . Hence, the right-hand side of equation (3.27) is in the range of $\mathcal{L}$ . To solve an equation of the form $\mathcal{L} g = R$ with $R\in \mbox{Range}(\mathcal{L})$ , we can use the pseudo-inverse of $\mathcal{L}^\perp$ and add a term which lies in the kernel of $\mathcal{L}$ , that is, we write

\begin{align*} g(\tau,X,w) = -\frac{1}{\mu(X,w)} R(\tau,X,w) + \alpha(\tau,X) T(X,w), \end{align*}

where $\alpha(\tau,X)$ is independent of w. Here, we assumed the correction term has no mass, $\bar g =0$ , and hence we choose $\alpha(\tau, X)$ in such a way that $\bar g =0$ . We solve (3.27) for the correction term, g, to obtain

\begin{align*}g = -\frac{1}{\mu} \Bigl[-\nabla \cdot({\bf E}_T \bar p) T + w \cdot \nabla (\bar p T) \Bigr] + \alpha T\end{align*}

with

(3.29) \begin{equation}\alpha(\tau,X) = +\int_V \frac{1}{\mu} \Bigl[-\nabla \cdot({\bf E}_T \bar p) T + w \cdot \nabla (\bar p T) \Bigr]dw.\end{equation}

Note that if the turning rate is independent of velocity w (i.e. $\mu(X,w)=\mu(X)$ ), then from (3.28) we note $\alpha=0$ and return to the standard case discussed in the introduction (see [Reference Hillen and Painter27]). We write g in the following form, more convenient for later manipulation:

\begin{eqnarray*}&g=& -\frac{1}{\mu}\Big[ T(w-{\bf E}_T ) \cdot \nabla \bar p + \bar p(w\cdot \nabla T) - \bar p_0 ((\nabla\cdot {\bf E}_T) T ) \Big] + \alpha T \\[3pt]&=& -\frac{1}{\mu}\Big[ T(w-{\bf E}_T) \cdot \nabla \bar p + \nabla T \cdot (w-{\bf E}_T) \bar p + {\bf E}_T \cdot \nabla T\bar p + \big[\nabla\cdot(w-{\bf E}_T) \big]\bar p T\Big] + \alpha T \\[3pt]&=& -\frac{1}{\mu} \Big[\nabla\cdot\big[(w-{\bf E}_T) \bar p T \Big] + {\bf E}_T\cdot \nabla T \bar p \Big] + \alpha T.\end{eqnarray*}

Then, the term in the square brackets of (3.26) becomes

(3.30) \begin{eqnarray} \int_V(w-\mathbf{E}_T) g dw &=& -\underbrace{\int_{V}(w-{\bf E}_T)\frac{1}{\mu}\nabla\big[(w-{\bf E}_T) \bar p T\big] dw}_{(I)}\end{eqnarray}

(3.31) \begin{eqnarray} && - \int_V \frac{1}{\mu} (w-{\bf E}_T) {\bf E}_T \cdot \nabla T \, dw\; \bar p + \int_V(w-{\bf E}_T) \alpha T \, dw \end{eqnarray}

(3.32) \begin{eqnarray}&=& -\nabla \cdot(\mathbb{D}_h \bar p ) + {\bf a}_h \bar p - \int_V (w-{\bf E}_T) {\bf E}_T \cdot \frac{\nabla T}{\mu} \, dw\; \bar p\,.\end{eqnarray}

In the above, we used the relation (3.17) between the integral (I) and the diffusion and drift terms $\mathbb{D}_h$ and ${\bf a}_h$ , respectively, as well as the fact that

\begin{align*} \int_V(w-{\bf E}_T) \alpha T dw = \alpha \left(\int_V w T dw - {\bf E}_T \int_V T dw \right) =0. \end{align*}

Combining these calculations with the macroscopic limit equation (3.26), we obtain the hyperbolic limit equation with correction term as:

(3.33) \begin{equation}\dfrac{\partial \bar{p}}{\partial \tau}+\nabla \cdot \big( ({\bf E}_T +\varepsilon {\bf a}_h) \bar{p}\big)=\varepsilon \nabla \cdot \nabla \cdot \Big( \mathbb{D}_h\bar{p}\Big)+\varepsilon\nabla \cdot\left(\int_V (w-{\bf E}_T) {\bf E}_T \cdot \frac{\nabla T}{\mu} \, dw\; \bar p\right),\end{equation}

where $\mathbb{D}_h$ and $\mathbf{a}_h$ are given by (3.15) and (3.16), respectively.

A particularly pleasing aspect of the above hyperbolic limit lies in its generalisation of the earlier parabolic limit. Specifically, for the case ${\bf E}_T=0$ , we obtain the same macroscopic quantities: (3.15) coincides with (3.12) and (3.16) coincides with (3.13). Moreover, for ${\bf E}_T=0$ , the rather clumsy integral correction term vanishes and we obtain a rescaled version of the parabolic limit equation (3.11):

(3.34) \begin{equation}\dfrac{\partial \bar{p}}{\partial \tau}+\varepsilon \nabla \cdot \big({\bf a}_h \bar{p}\big)=\varepsilon \nabla \cdot \nabla \cdot \Big( \mathbb{D}_h\bar{p}\Big)\,.\end{equation}

A further scaling of time with $\varepsilon$ then reproduces the parabolic limit (3.11).

3.4 Time-dependent turning distribution

For time-varying tissues, q, $\mu$ and C will all depend on time. Rescaling (2.1) and (2.3) with $X=\varepsilon x$ , regardless of using timescale $\tau=\varepsilon^2 t $ or $\tau= \varepsilon t$ , and comparing equal orders of $\varepsilon$ allows us to obtain the leading order function of (3.2), that is,

\begin{align*}p_0(\tau,X,w)=\bar{p}(\tau,X)T(\tau,X,w),\end{align*}

where the equilibrium distribution is now also time-dependent, $T(\tau,X,w)=\dfrac{q(\tau,X,w)}{\mu(\tau,X,w)}$ .

The diffusive limit does not change, but the hyperbolic limit has different correction terms. Proceeding as in the previous section, we find that the correction g is a solution to

(3.35) \begin{equation}\mathcal{L}g=\dfrac{\partial}{\partial \tau} (\bar{p}T)+w\cdot\nabla \Big( \bar{p}T\Big) +\mathcal{O}(\varepsilon).\end{equation}

Let us suppose again that g is of the form:

\begin{align*}g=-\dfrac{R(\tau,X,w)}{\mu(\tau,X,w)}+\alpha(\tau,X)T(\tau,X,w),\end{align*}

where again $\alpha(\tau,X)$ does not depend on w. Inverting (3.35), we find

\begin{align*}g=-\dfrac{1}{\mu}\dfrac{\partial \bar{p}}{\partial \tau} T+\bar{p}\dfrac{\partial T}{\partial \tau} +w\cdot\nabla \Big( \bar{p}T\Big)+\alpha(\tau,X) T(\tau,X,w)\end{align*}

and then

\begin{align*}g=-\dfrac{1}{\mu}\Big[\Big(w T-{\bf E}_T T\Big)\cdot \nabla \bar{p}+\Big( w\cdot\nabla T-\nabla \cdot{\bf E}_T T+\dfrac{\partial T}{\partial \tau} \Big)\bar{p} \Big]+\alpha(\tau,X) T(\tau,X,w) .\end{align*}

Since we assumed g carries no mass and have $\dfrac{\partial}{\partial \tau}\bar{T}=0$ , we obtain the same $\alpha$ as in (3.29). Concluding, equation (3.33) becomes in this case:

(3.36) \begin{align}\dfrac{\partial \bar{p}}{\partial \tau}+\nabla \cdot \big( \bar{p}({\bf E}_T +\epsilon a_h) \big)&=\varepsilon \nabla \cdot \nabla \cdot \Big( \mathbb{D}_h\bar{p}\Big)+\varepsilon\nabla \cdot\left(\int_V (w-{\bf E}_T) {\bf E}_T \cdot \frac{\nabla T}{\mu} \, dw\; \bar p\right) \nonumber\\& \quad + \varepsilon \nabla \cdot \Big(\dfrac{\partial}{\partial \tau}{\bf E}_T \bar{p}\Big).\end{align}

4.1 Dependencies of turning rates and turning kernels

1. (1) Suppose $\mu$ does not depend on w. The run time distribution, T, then becomes

   \begin{align*} T(x,w) = \frac{q(x,w)}{C(x)\mu(x)}, \quad \mbox{with }\quad C(x) = \int_V \frac{q(x,w) }{\mu(x)} dw = \frac{1}{\mu(x)}\,. \end{align*}
   Hence, in this case,
   \begin{align*} T(x,w) = q(x,w)\end{align*}
   and the model reduces to previously studied cases (e.g. see [Reference Hillen24,Reference Painter and Hillen38]). The diffusion matrix here is given by the variance–covariance matrix of the fibre network distribution, q(x,w), divided by the turning rate:
   \begin{align*}\mathbb{D}(x) =\int_V w\otimes w\, \frac{T(x,w)}{\mu(x)} dw = \frac{1}{C(x)\mu(x)^2}\int_V w\otimes w q(x,w) dw = \dfrac{\mathbb{V}_q(x)}{\mu(x)}.\end{align*}
   The drift velocity (3.13), meanwhile, is
   (4.2) \begin{equation}{\bf a}(x)=- \dfrac{\nabla \mu(x)}{\mu^2(x)}\mathbb{V}_q(x).\end{equation}
   In this case, if q is even and then ${\bf E}_T$ vanishes, we obtain a parabolic limit equation (3.11) as:
   (4.3) \begin{equation}\dfrac{\partial u}{\partial t} - \nabla \cdot \left[\mathbb{V}_q \dfrac{\nabla\mu}{\mu^2} u \right] = \nabla \cdot \nabla \cdot \left(\dfrac{\mathbb{V}_q}{\mu} u\right).\end{equation}
   This can again be written in Fickian form, where we obtain
   (4.4) \begin{equation}\dfrac{\partial u}{\partial t}= \nabla \cdot \left(\dfrac{\mathbb{V}_q}{\mu}\nabla u \right)+\nabla \cdot \left[ \dfrac{\nabla \cdot \mathbb{V}_q}{\mu} u\right]\,.\end{equation}
   When $\mu$ is spatially heterogeneous, the diffusion process (4.3) is a fully anisotropic process with advection given by the taxis term (4.2). This leads the dynamics towards decreasing values of $\mu$ . The turning rate $\mu$ also scales the diffusivity, with large diffusivity for small turning rates and vice versa. Hence, frequently turning particles will not show large diffusive spread compared to those turning less often. The drift in the direction of the negative gradient of $\mu$ indicates a drift away from regions of low diffusion towards regions of high diffusion.

2. (2) If we further assume that $\mu$ is constant and consider $q=q(x,w)$ , the Fokker–Planck equation (3.11) leads to

   (4.5) \begin{equation}\dfrac{\partial u}{\partial t}= \frac{1}{\mu} \nabla \cdot \nabla \cdot(\mathbb{V}_q u).\end{equation}
   When the oriented habitat is not spatially homogeneous, (4.5) describes a fully anisotropic process where diffusion is again linked to an anisotropic environment (which determines the direction of anisotropy) and to the frequency of reorientations (which determines the intensity of diffusion).

3. (3) Suppose now $\mu$ is constant and further assume $q=q(w)$ , that is, the network is spatially homogeneous but can still describe an oriented environment. Then the diffusion equation (4.5) becomes

   (4.6) \begin{equation}\dfrac{\partial u}{\partial t} = \frac{1}{\mu} \nabla \cdot (\mathbb{V}_q \nabla u)\,.\end{equation}
   Equation (4.6) describes diffusion in a spatially homogeneous environment, where anisotropy is due to the presence of a dominant alignment in the environment.

4. (4) Suppose both turning rates and turning kernels depend on w but not x, that is, $\mu = \mu (w)$ and $q=q(w)$ . If we further assume ${\bf E}_T=0$ , as in this case $\nabla\mu(X,w)=0$ , then the macroscopic drift from (3.13) satisfies ${\bf a}=0$ and the macroscopic diffusion (3.12) is

   (4.7) \begin{equation}\mathbb{D} = \int_V w\otimes w \frac{T(X,w)}{\mu(w)} dw .\end{equation}
   We use this case later to analyse the movement patterns of cells migrating on fabricated anisotropic surfaces, such as the fibronectin stripe arrangements devised in [Reference Doyle, Wang, Matsumoto and Yamada16] (see the schematic in Figure 1 and the simulations in Figure 7).

5. (5) Finally, suppose a constant fibre distribution, that is,

   (4.8) \begin{equation}q=1/2\pi,\end{equation}
   but $\mu = \mu (w,x,t)$ , that is, the turning rate can depend on the orientation of the environment. If we suppose ${\bf E}_T=0$ , we are again in the previous case, but the eventual anisotropy will be the direction orthogonal to the dominant direction of $\mu$ , as it is the orientation of cells that tend to turn more frequently.

4.2 Isotropy and anisotropy

Generally, we observe an anisotropic diffusion process (3.12), a consequence of both cells orienting with respect to the network alignment via q(x,w) and the direction-dependent turning frequency $\mu(x,w)$ . Specifically, diffusive anisotropy is encapsulated through the second moment of the distribution:

\begin{align*} \frac{T(x,w)}{\mu(x,w)} = \frac{q(x,w)}{C(x) \mu^2(x,w)}. \end{align*}

The nominator/denominator localisations of q and $\mu$ naturally suggest that turning direction and turning rate might have opposing impact on macroscopic dynamics.

1. (1) For example, suppose that

   (4.9) \begin{equation}q(x,w) \sim \mu^2(x,w)\,.\end{equation}
   Here, we have $T(x,w)\sim \mu(x,w)$ and obtain a macroscopic drift term of the form ${\bf E}_T(x) \sim \int_V w\mu(x,w) dw.$ The macroscopic diffusion tensor (3.12) in this case will be isotropic:
   \begin{align*}\mathbb{D} \sim \displaystyle \int_V w\otimes w dw = c I.\end{align*}
   Hence, for $\int w\mu(x,w)dw \neq 0$ , the macroscopic regime will be isotropic and drift-dominated.

2. (2) Suppose instead

   (4.10) \begin{equation}q(x,w) \sim \mu(x,w)\,.\end{equation}
   Here, T(x,w) is constant and the macroscopic regime will always be diffusive with ${\bf E}_T=0$ . The diffusion might be anisotropic, since
   \begin{align*} \mathbb{D}\sim \displaystyle \int_V w\otimes w \dfrac{1}{q(x,w)} dw.\end{align*}

3. (3) It is also relevant to consider the case in which

   (4.11) \begin{equation}q(x,w) \sim \frac{1}{\mu(x,w)}\,,\end{equation}
   because in many biological cases, cells are likely to both choose the fibre direction and turn less frequently when moving in the direction of dominating alignment. In this case, we have both a drift-driven phenomenon, with ${\bf E}_T=\int w\frac{1}{\mu^2(x,w)}dw$ , and the diffusion can be anisotropic according to
   \begin{align*} \mathbb{D}\sim \displaystyle \int_V w\otimes w q^3(x,w) dw.\end{align*}

4.3 Comparison to diffusion of passive particles

The somewhat curious diffusive terms derived here arise due to the study of active movers: biological particles, such as cells and organisms, generate their own movement energy and are therefore unconstrained by energy and momentum conservation as demanded by classical physics. For passive movers (movers simply transported by a fluid environment or some other stream), the situation is different: in that context, an equation of the form (4.6) corresponds to the equation describing diffusion in an anisotropic but homogeneous environment, that is, without spatial variation; equation (4.5) corresponds to the equation derived in [Reference Chapman8], describing particle diffusion in a heterogeneous environment; equation (4.3) corresponds to the case in which $\mu$ depends on the position, and it is then the only case in which the macroscopic equation has that particular mixed structure. Within the diffusion theory of physical particles, this mixed structure arises when describing a so-called thermal effect, that is, when there is a temperature gradient [Reference Wereide55]. Here, the variable turning rate can be viewed somewhat analogously to temperature in a system of physical particles, ^{Footnote 1} and the two factors influencing the dynamics are the spatial heterogeneity and the adaptability to environmental heterogeneity.

An example of adaptability in a biological context is provided in [Reference Cho and Kim10,Reference Chung, Kim, Kwong and Yoon11], addressing starvation-driven diffusion. Here, starvation induces the organism to increase motility and find a better environment, even if it is not known where that may be, that is, the migration is unbiased. We further note that the three forms of diffusion equation (4.3),(4.5) and (4.6) can be derived from space-jump processes, where variation in jumping depends on environment assessment of the current location, between locations or the destination [Reference Chung, Kim, Kwong and Yoon11,Reference Lutscher and Hillen31,Reference Othmer and Stevens33], with spatially homogeneous jumping times. When both the jumping time and length are spatially heterogeneous, a mixed structure similar to (3.11) arises. We note that the coefficients appearing in (4.3),(4.5) and (4.6) may also depend on time, since both q and $\mu$ could be time-dependent. This would be relevant, for example, in the context of mesenchymal motion where there is significant ECM remodelling [Reference Hillen24,Reference Painter36].

5.1 Arrangements

We specify q and $\mu$ using bimodal von Mises distributions, a standard circular distribution with known analytical forms for the first and second moments (e.g. see [Reference Hillen, Murtha, Painter and Swan25]). For notational convenience, we introduce the following shorthand for a bimodal von Mises distribution, with a given concentration parameter $k>0$ and a given unit vector $\theta$ :

(5.1) \begin{equation}{\mathsf{bvM}}_{k, \theta} (w) = \frac{1}{4\pi I_0(k)} \left(e^{k w\cdot\theta} + e^{-k w\cdot \theta} \right).\end{equation}

To further simplify notation, we identify a unit vector by its angle, that is, writing $\theta= (\cos\theta, \sin\theta)^T$ .

Fibres are arranged in two principal patterns, schematised in Figure 2. We base q and $\mu$ on (5.1), stipulating functions $\theta_q$ and $k_q$ for q and $\theta_\mu$ and $k_\mu$ for $\mu$ . $k_q = 0$ , therefore corresponds to unaligned fibres and large $k_q$ generates highly aligned fibres along the axis $\theta_q, \theta_q+\pi$ . Tests are devised as follows.

• • Test 1. Here we set $\theta_q = \theta_{\mu} = \pi/4$ , $k_{q}=50e^{-0.25((x-5)^2+(y-5)^2)}$ and $\mu$ as one of $\mu \sim constant, \sim q, \sim q^2, \sim 1/q.$

• • Test 2. We set $\theta_q$ , $\theta_{\mu}$ and $k_q$ as in Test 1, but shift $k_{\mu}$ to $50e^{-0.25((x-4.5)^2+(y-4.5)^2)}$ .

• • Test 3. Setting $\Omega_{XY}=\Omega_X \cap \Omega_Y$ , where $\Omega_X=\lbrace (x,y): 4 \leq x \leq 6\rbrace$ , $\Omega_Y \lbrace (x,y): 4\leq y \leq 6 \rbrace$ , we consider

  (5.2) $$q(x,w) = \left\{ {\matrix{ {{1 \over {2\pi }}\quad {\rm{on}}\quad \Omega {\rm{ - }}{\Omega _{{\rm{XY}}}},} \hfill \cr {bv{M_{{k_q},{\theta _q}}}\quad {\rm{on}}\quad {\Omega _{\rm{X}}},} \hfill \cr {bv{M_{{k_q},\theta _q^ \bot }}\quad {\rm{on}}\quad {\Omega _{\rm{Y}}},} \hfill \cr {{1 \over 2}(bv{M_{{k_q},{\theta _q}}} + bv{M_{{k_q},\theta _q^ \bot }}){\mkern 1mu} {\rm{on}}{\mkern 1mu} {\Omega _{{\rm{XY}}}},} \hfill \cr } } \right.\mu (x,w) = \left\{ {\matrix{ {{1 \over {2\pi }}\quad {\rm{on}}\quad \Omega {\rm{ - }}{\Omega _{{\rm{XY}}}},} \hfill \cr {bv{M_{{k_\mu },{\theta _\mu }}}\quad {\rm{on}}\quad {\Omega _{\rm{X}}},} \hfill \cr {bv{M_{{k_\mu },\theta _\mu ^ \bot }}\quad {\rm{on}}\quad {\Omega _{\rm{Y}}},} \hfill \cr {{1 \over 2}(bv{M_{{k_\mu },{\theta _\mu }}} + bv{M_{{k_\mu },\theta _\mu ^ \bot }}){\mkern 1mu} {\rm{on}}{\mkern 1mu} {\Omega _{{\rm{XY}}}}.} \hfill \cr } } \right.$$
  Specifically, we will set $\theta_q = \pi$ and $k_q = 50$ to form the cross-configuration on the right of Figure 2, where away from the cross-fibres are isotropic and on the horizontal (vertical) arm fibres are aligned in the horizontal (vertical) direction. At the centre, fibres cross. For $\mu$ , we consider related but subtly distinct forms, allowing distinct turning rate to turning distribution behaviour.

Figure 2.

Left: Arrangement and intensity of fibre orientation for Tests 1 and 2. The maximum of $\mu$ is positioned at the red square for Test 1 and the green circle for Test 2. Right: Arrangement and intensity of fibre orientation for Test 3. According to the precise simulation, we again shift the position of the maximum of $\mu$ from the red square to the green circle. Note that our model assumes that the fibre orientation directly informs the turning kernel, q, as explained in Section 2.

5.2 Simulations

Test 1 was designed to explore the extent to which isotropic or anisotropic behaviour arises with the choices of q and $\mu$ , for example, related as in equations (4.9), (4.10) or (4.11). Initialising according to Figure 3(a), we plot the macroscopic density at $T=7.5$ in Figure 3(b–e) for (b) $\mu \sim constant = 1$ , (c) $\mu \sim q$ , (d) $\mu \sim q^2$ , (e) $\mu \sim 1/q$ . In (f), we set q constant but maintain anisotropy through $k_\mu = 50e^{-0.25((x-5)^2+(y-5)^2)}$ and $\theta_\mu=\pi/4$ .

Figure 3.

Test 1. Anisotropic/isotropic spread according to the choices of q and $\mu$ . See text for details. (a) Initial macroscopic density $\bar{p}(0,x)$ , (b)–(f) macroscopic density $\bar{p}(t,x)$ at $t=7.5$ . The white curve in (b)–(f) is the level set defined by $\bar{p}(t,x)=0.1$ . In (b–e) q is as described in the text and the turning rate is given by: (b) $\mu=1$ , (c) $\mu=q$ , (d) $\mu=\sqrt{q}$ and (e) $\mu=1/q$ . In (f), $q=1/2\pi$ while $\mu$ is given by (5.1), with $\theta_{\mu}=\pi/4$ and $k_{\mu}=50e^{-0.25((x-5)^2+(y-5)^2)}$ .

In particular, for Figure 3(b), we recover the original model of [Reference Hillen24] and observe anisotropic diffusion according to the dominating fibre alignment. In Figure 3(c), we set $\mu=q$ , that is, situation (4.10). There is no drift, but anisotropy arises through the directional bias in q. This generates a conflict in which cells preferentially choose the dominating fibre alignment but, when facing those directions, turn more frequently. Consequently, the anisotropy becomes orthogonal to the dominating fibre alignment. In Figure 3(d), we consider case (4.9) by choosing $\mu=\sqrt{q}$ . Again, preferential movement along an axis is countered by increased turning, but the weighting now generates isotropic diffusion. Further, there is a non-trivial drift ${\bf E}_T$ . The dynamics are similar to those in Figure 3(c), yet the pattern is more symmetric due to the isotropic diffusion. In Figure 3(e), we consider case (4.11), that is, where $\mu=1/q$ . Cells now turn less frequently along the directions of preferential alignment and, intuitively, we should expect enhanced anisotropic diffusion along certain axes. Simulations confirm this, with the cells remaining even more tightly aligned along the dominating fibre orientation as compared to Figure 3(b). Finally, in Figure 3(f), the fibre network does not impact on orientation, but cells oriented along the axis ( $\pi/4, 5\pi/4$ ) turn more frequently. Diffusion remains anisotropic as predicted by (4.8).

Test 2 was designed to show the taxis induced by the spatial variability of the turning rate $\mu$ and simulation results are reported in Figure 4. Specifically, we shift the peak of the turning rate away from the centre of the domain, to the point represented by the green dot. There is a subsequent tendency of cells to avoid this new location, coherent with equation (4.2) indicating greater diffusion with higher values of the turning frequency (Figure 4).

Figure 4.

Test 2. Taxis induced through variable turning rates. q is as described in the left of Figure 2 (see text for details), while the peak of the turning frequency $\mu$ is shifted to centre on the green dot. Simulations plot the macroscopic density at successive times $t = 2.5, 5, 7.5$ .

Test 3 extends these analyses to a more complicated environment, see Figure 2. Figure 5(a) plots the initial (macroscopic) cell distribution for the experiments (b) and (c), while Figure 5(d) shows the distribution used for (e) and (f). In Figure 5(b), we set $\mu \sim q$ so that cells orient and migrate along fibres but this is counterbalanced by turning more frequently when moving in those directions. Spread is subsequently inhibited by the cross-arrangement. In Figure 5(c), we consider instead $\theta_\mu=\theta_q^{\bot}$ , so that particles both follow the fibres and turn less frequently when moving in their direction. As expected, there is a clear tendency of cells to follow the cross-structure. The second row performs the same simulations, but relocating the initial cell distribution (now centred at (4,4)) and the peak of $k_\mu$ to (3,3). The latter generates an advection away from this point, a consequence of the decreasing gradient of $\mu$ experienced by cells. Figure 5 (e) shows even more clearly the inhibition resulting from fibres at the centre of the cross, while in (f) cells are seen to spread rapidly along the arms once it has been reached. The different times in Figure 5(c) and (f) were chosen as to best represent the dynamics in the two different settings.

Figure 5.

Test 3. Dynamics for q given by the cross-structure in Figure 2 (left). (a) Initial macroscopic density (centred at (5,5)) for the simulations represented in (b,c); (d) initial macroscopic density (centred at (4,4)) for simulations represented in (e,f). q is given by (5.2), see text for details, while $\mu$ has a similar structure, but for (b,e) $\theta_{\mu} = \theta_{q}$ and for (c,f) $\theta_\mu=\theta_q^{\bot}$ .

6.1 Control case

As a control, consider a constant turning rate $\mu$ and, in turn, a constant mean travel time $\tau$ . Here, the macroscopic diffusion tensor is computed from formula (4.6) as:

\begin{align*} \mathbb{D} = \frac{1}{\mu} \mathbb{V}_q.\end{align*}

Under Case A, the middle region is completely non-oriented, hence $q=\frac{1}{2\pi}$ (the uniform distribution). Then,

(6.1) \begin{equation}\mathbb{V}_q = \int_V w\otimes w \frac{1}{2\pi} dw = \frac{2\pi}{2} \mathbb{I}\frac{1}{2\pi} = \frac{1}{2} \mathbb{I},\end{equation}

where $\mathbb{I}$ denotes the identity matrix. The criss-cross configuration of Case B can be described by combining two bimodal von Mises distributions, in perpendicular directions $e_1 = (1,0)$ and $e_2 = (0,1)$ but with the same concentration parameter k:

(6.2) \begin{equation}q(w) = \frac{1}{2}\left( {\mathsf{bvM}}_{k, e_1} + {\mathsf{bvM}}_{k, e_2}\right).\end{equation}

Following standard calculations (e.g. see [Reference Hillen, Murtha, Painter and Swan25]), we obtain

\begin{align*} \mathbb{V}_q = \frac{1}{2}\left( 1 -\frac{I_2(k)}{I_0(k)}\right) \mathbb{I} + \frac{I_2(k)}{I_0(k)} \frac{1}{2}\left(e_1 e_1^T + e_2 e_2^T\right) .\end{align*}

Now,

\begin{align*} e_1 e_1^T = \left( \begin{array}{cc} 1 & 0 \\ 0 & 0 \end{array}\right), \qquad \mbox{and} \quad e_2 e_2^T = \left( \begin{array}{cc} 0 & 0 \\ 0 & 1 \end{array}\right),\end{align*}

so

(6.3) \begin{equation}\mathbb{V}_q = \frac{1}{2}\mathbb{I} -\frac{I_2(k)}{2I_0(k)} \mathbb{I} + \frac{I_2(k)}{2I_0(k)} \mathbb{I} = \frac{1}{2} \mathbb{I},\end{equation}

which coincides exactly with the calculation for Case A, that is, (6.1). Therefore, under a constant turning rate, there should be no macroscopic difference between Case A and Case B, even if cells bias their orientation along the criss-crossed fibres.

6.2 Direction-dependent turning

We now consider w-dependence in the turning rate, that is, $\mu(w)$ . Specifically, we choose $\mu$ such that the rate of turning is reduced if cells are migrating along the direction of dominating alignment. For the analysis, we focus only on the central region, so both Case A and Case B can be regarded as spatially homogeneous (i.e. not depending on x) and we can use equation (4.7):

\begin{align*} \mathbb{D} = \int w\otimes w \frac{T}{\mu} dw, \qquad T = \frac{q}{C \mu}, \qquad C(x) = \int_V \frac{q}{\mu} dw. \end{align*}

We again choose q to be a combination of bimodal von Mises distributions in the two perpendicular directions $e_1$ and $e_2$ , as given in (6.2). However, now we assume that $\mu\sim q^{-1}$ , so that

\begin{align*}\mu=\frac{1}{2\pi} \dfrac{2}{\left({{\mathsf{bvM}}_{k, e_2}}+{\mathsf{bvM}}_{k, e_1}\right)}\,,\end{align*}

where we chose the normalisation constant to be $({2\pi)}^{-1}$ such that in the isotropic limit of $k\to 0$ we obtain

\begin{align*} \lim_{k\to 0} \mu = \frac{1}{2\pi} \frac{2}{\frac{1}{2\pi}+\frac{1}{2\pi} } = 1.\end{align*}

Then,

(6.4) \begin{equation}\frac{T}{\mu} = \frac{\pi^2}{2 C } \left({{\mathsf{bvM}}_{k, e_2}}+{\mathsf{bvM}}_{k, e_1}\right)^3.\end{equation}

For this choice of $\mu$ and q, we can compute the normalisation constant C as:

(6.5) \begin{eqnarray} C(x) &=& \int_V\frac{q}{\mu} dw = \int_V\frac{\pi}{2} \left({\mathsf{bvM}}_{k, e_2}+{\mathsf{bvM}}_{k, e_1}\right)^2 dw \nonumber \\ &=& \frac{\pi}{2}\frac{1}{(4\pi I_0(k))^2} \int_V 4\pi I_0(2k) (\mathsf{bvM}_{2k, e_1} +\mathsf{bvM}_{2k, e_2}) + 4 \nonumber \\ && + 8\pi I_0(\sqrt{2}k) \left(\mathsf{bvM}_{\sqrt{2} k, \frac{e_1+e_2}{\sqrt{2}}} + \mathsf{bvM}_{\sqrt{2}k, \frac{e_1-e_2}{\sqrt{2}}}\right) dw \nonumber \\ &=& \frac{1}{4 I_0(k)^2 }\bigl(I_0(2k) + 1 + 2 I_0(\sqrt{2} k) \bigr). \end{eqnarray}

We note the isotropic limit $\lim_{k\to 0} C(x) = 1$ , since $I_0(0)=1$ .

Next, we compute the third power in (6.4)

\begin{eqnarray*}\frac{T}{\mu} &=& \frac{\pi^2}{2 C} \left( \frac{I_0(3k)}{I_0(k)} \Bigl[\mathsf{bvM}_{3k,e_1}+\mathsf{bvM}_{3k,e_2}\Bigr] + 3\Bigl[{\mathsf{bvM}}_{k, e_1}+{\mathsf{bvM}}_{k, e_2}\Bigr]\right.\\&& + 3\Bigl[ \mathsf{bvM}_{k, 2e_1+e_2} +\mathsf{bvM}_{k,2e_1-e_2} + \mathsf{bvM}_{k, e_1+2e_2} + \mathsf{bvM}_{k,e_1-2e_2} \\&& + 2({\mathsf{bvM}}_{k, e_1} + {\mathsf{bvM}}_{k, e_2})\Bigr]\Bigr)\end{eqnarray*}

Notably, vectors $2e_2+ e_1$ etc are not unit vectors so we rescale

\begin{align*} \zeta_1=\frac{1}{\sqrt{5}}(2,1)^T\,,\qquad \zeta_2=\frac{1}{\sqrt{5}}(1,-2)^T\,,\qquad \xi_1=\frac{1}{\sqrt{5}}(2,-1)^T\,,\qquad \xi_2=\frac{1}{\sqrt{5}}(1,2)^T.\end{align*}

For the example $2e_1+e_2$ , this gives

\begin{align*} \mathsf{bvM}_{k,2e_1+e_2} = \frac{I_0(\sqrt{5} k)}{I_0(k)} \mathsf{bvM}_{\sqrt{5}k, \zeta_1}, \end{align*}

and similar for the other terms. With unit vectors everywhere which, moreover, are pairwise perpendicular ( $\zeta_1\cdot \zeta_2 =0, \, \xi_1\cdot\xi_2=0$ ), we have

\begin{eqnarray*}\frac{T}{\mu} &=& \frac{\pi^2}{2 C} \left( \frac{I_0(3k)}{I_0(k)} \Bigl[\mathsf{bvM}_{3k,e_1}+\mathsf{bvM}_{3k,e_2}\Bigr] + 3\Bigl[{\mathsf{bvM}}_{k, e_1}+{\mathsf{bvM}}_{k, e_2}\Bigr]\right.\\&& + 9 ({\mathsf{bvM}}_{k, e_1}+{\mathsf{bvM}}_{k, e_2})+ 3\frac{I_0(\sqrt{5} k)}{I_0(k)} \left(\mathsf{bvM}_{\sqrt{5}k, \zeta_1} +\mathsf{bvM}_{\sqrt{5}k, \zeta_2}\right) \\&& + 3\frac{I_0(\sqrt{5} k)}{I_0(k)} \left(\mathsf{bvM}_{\sqrt{5}k, \xi_1} +\mathsf{bvM}_{\sqrt{5}k, \xi_2}\right)\Bigr).\end{eqnarray*}

Noting that the second moment of bimodal von Mises distributions with pairwise perpendicular unit vectors is $\frac{1}{2}\mathbb{I}$ (see (6.3)), we find a diffusion tensor:

\begin{align*} \mathbb{D} = d(k) \mathbb{I}, \qquad d(k) = \frac{\pi^2}{2 C (4\pi I_0(k))^2} \left( \frac{I_0(3k)}{I_0(k)} + 6 \frac{I_0(\sqrt 5 k)}{I_0(k)} + 9\right).\end{align*}

Substituting the normalisation constant C from (6.5), we obtain

\begin{align*}d(k) = \frac{I_0(3k) + 9 I_0(k) + 6 I_0(\sqrt{5}k) } {8 I_0(k) ( I_0(2k) + 2I_0(\sqrt{2} k) + 1)}. \end{align*}

Again we consider the isotropic limit:

\begin{align*} \lim_{k\to 0} d(k) = \frac{1 + 9 + 6}{ 8(1+2+1) } = \frac{1}{2}, \end{align*}

which has the correct scaling as for the isotropic case. The diffusion coefficient d(k) is plotted in Figure 6, where we observe that for small k there is negligible change but for $k>3$ we see clear and sustained increase for the diffusion coefficient.

Figure 6.

Diffusion coefficient d(k) as function of k for $k\in[0,10]$ on the left and $k\in[0,100]$ on the right. The diffusion coefficient grows very slowly for small k, but then grows more rapidly for $k \gtrsim 3$ . For large k, the increase is approximately $\propto \sqrt{k}$ .

6.3 Transport model simulations

The above analysis indicates that a direction-dependent turning rate coupled to a criss-cross fibre network substantially increases the diffusion, compared to either a constant turning rate or a completely isotropic network. To simulate this situation, we consider a rectangular domain $\Omega=[0,15]\times[0,5]$ with $V=\mathbb{S}^1$ ( $s=1$ ) and define local environments corresponding to those illustrated in Figure 1. Specifically, to describe the parallel alignment in the left and right regions, we set $q(x,w)=e^{k w\cdot \theta}/2\pi I_0(k)$ with $k=25$ and $\mu(x,w)=1$ when $x<5$ or $x>10$ . To replicate the completely isotropic scenario of Case A, we define the central region ( $5\le x\le 10$ ) by $q(x,w)=1/2\pi$ , while for the criss-cross network of Case B we choose $q(x,w)=\left({\mathsf{bvM}}_{k, e_1}+{\mathsf{bvM}}_{k, e_2}\right)/2$ in the central region and take the anisotropy constant k to be a variable model parameter. Consequently, for Case A $\mu(w)=1$ for the full domain, while in Case B we chose $\mu(x,w)=1/(2\pi q)$ if $5\le x\le 10$ . We initialise the density distribution as uniformly distributed in $\mathbb{S}^1$ with the initial macroscopic density ( $\bar{p}_0$ ) as defined in Figure 3(a), but centred at the coordinate $(3.5,2.5)$ , that is, inside the left-most region of highly aligned fibres. Simulations of the corresponding transport equations (2.1)–(2.9) are shown in Figure 7. The first row illustrates results for Case A, while subsequent rows are for Case B with increasing values $k = 0,3,10,25$ , respectively.

Figure 7.

Cell movement on fabricated anisotropic surfaces. For all simulations, we set $q(x,w)=e^{k w\cdot \theta}/(2\pi I_0(k))$ with $k=25$ and $\mu(x,w)=1$ for $x<5$ and $x>10$ . First row: Case A with $q(x,w)=1/(2\pi)$ if $5\le x\le 10$ and $\mu(w)=1$ on the full domain. Second to final row: $q(x,w)=\left({\mathsf{bvM}}_{k, e_1}+{\mathsf{bvM}}_{k, e_2}\right)/2$ and $\mu(x,w)=1/(2\pi q)$ if $5\le x\le 10$ , so that $s=1$ . In particular: (second row) $k=0$ , (third row) $k=3$ , (fourth row) $k=10$ , (final row) $k=25$ .

Under Case A (Figure 7 first row), as cells reach the isotropic region, we observe a gradual diffusive-like spread, consistent with the earlier analysis. Under Case B, but setting $k=0$ (Figure 7 second row), generates equivalent behaviour: in line with the prediction that isotropic criss-cross networks do not alter the macroscopic dynamics when the turning rate is constant. For $k>0$ , however, we observe distinct behaviour. This is minimal for small k (e.g. $k=3$ , Figure 7 third row) but becomes clear for large k, for example, $k=10$ (fourth row) and $k=25$ (fifth row), respectively. A noteworthy phenomenon lies in the ‘droplet’ detaching from the main swarm for high anisotropy parameter values ( $k=10, 25$ ). Here, a fraction of invaders maintain the left to right direction on reaching the central region, detaching from the main swarm. This observation is unexpected and could be of interest to confirm experimentally.

As a final remark, we note that while the simulations have been performed within a non-dimensional setting, direct comparison with the data of [Reference Doyle, Wang, Matsumoto and Yamada16] allows a more precise comparison of the spatial and temporal timescales of movements. Typical cell migration speeds reported in [Reference Doyle, Wang, Matsumoto and Yamada16] are of the order 1 $\mu\mbox{m}\,{min}^{-1}$ , while the central square region is of approximate side length $40\,\mu\mbox{m}$ . Based on these measurements, the frames at $t=7.5$ , 20 and 30 correspond to approximately 60, 160 and 240 min, respectively. Accordingly, for the Case B setting under higher k values, the majority of cells traverse the central region within a few hours. This appears in line with the videos of cell movements in [Reference Doyle, Wang, Matsumoto and Yamada16].