Displacement of transport processes on networked topologies

TL;DR

This paper introduces the winding distance to measure particle displacement on networks, revealing a topology-induced transition from diffusive to ballistic behavior for Brownian particles and extending analysis to anomalous diffusion.

Contribution

It proposes the winding distance as a new measure and derives a macroscopic model showing topology-driven displacement transitions in networked transport.

Findings

Confinement can induce a diffusive to ballistic transition in displacement.

A topological condition predicts when the transition occurs.

Network topology influences long-time behavior of anomalous diffusion.

Abstract

Consider a particle whose position evolves along the edges of a network. One definition for the displacement of a particle is the length of the shortest path on the network between the current and initial positions of the particle. Such a definition fails to incorporate information of the actual path the particle traversed. In this work we consider another definition for the displacement of a particle on networked topologies. Using this definition, which we term the winding distance, we demonstrate that for Brownian particles, confinement to a network can induce a transition in the mean squared displacement from diffusive to ballistic behaviour, $\langle x^2(t) \rangle \propto t^2$ for long times. A multiple scales approach is used to derive a macroscopic evolution equation for the displacement of a particle and uncover a topological condition for whether this transition in the mean squared displacement will occur. Furthermore, for networks satisfying this topological condition, we identify a prediction of the timescale upon which the displacement transitions to long-time behaviour. Finally, we extend the investigation of displacement on networks to a class of anomalously diffusive transport processes, where we find that the mean squared displacement at long times is affected by both network topology and the character of the transport process.

Figures

Equations

$${\rho }_k^{+}=\sum _{\nu \in {\mathcal{V}}_L}\frac{d_{\nu }-1}{d_{\nu }}\times p_k(\nu ).$$

$${\rho }_{\infty }^{+}=\sum _{\nu \in {\mathcal{V}}_L}\frac{d_{\nu }-1}{d_{\nu }}\times \frac{d_{\nu }}{{\sum }_{\omega \in {\mathcal{V}}_L}{d}_{\omega }}=1-\frac{1}{\bar{d}},$$

$$\frac{\mathrm{d}{p}_n}{\mathrm{d}t}={\lambda }_{(n-1)\to n}{p}_{n-1}+{\lambda }_{(n+1)\to n}{p}_{n+1}-\left({\lambda }_{n\to (n-1)}+{\lambda }_{n\to (n+1)}\right)p_n.$$

$$\frac{\partial P_n(x,t)}{\partial t}$$

$$\frac{\partial P_i}{\partial t}(x,t)$$

$${ϵ}^2\frac{\partial \vec{P}}{\partial t}(x,t)=\mathbf{A}\vec{P}+ϵ\mathbf{B}\frac{\partial \vec{P}}{\partial x}+{ϵ}^2\mathbf{C}\frac{{\partial }^2\vec{P}}{\partial x^2}+\mathcal{O}\left({ϵ}^3\right),$$

$$\color[rgb]{0,0,0}{\mathbf{A}=\begin{pmatrix}-1&1/2&0&\cdots&0&0&\rho^{+}\\ 1/2&-1&1/2&\cdots&0&0&0\\ 0&1/2&-1&\cdots&0&0&0\\ 0&0&1/2&\cdots&0&0&0\\ &\vdots&&\ddots&&\vdots&\\ 0&0&0&\cdots&1/2&0&0\\ 0&0&0&\cdots&-1&1/2&0\\ 0&0&0&\cdots&1/2&-1&\rho^{-}\\ 1/2&0&0&\cdots&0&1/2&-1\end{pmatrix},}$$

$$\color[rgb]{0,0,0}{\mathbf{B}=\begin{pmatrix}0&1/2&0&\cdots&0&0&-\rho^{+}\\ -1/2&0&1/2&\cdots&0&0&0\\ 0&-1/2&0&\cdots&0&0&0\\ 0&0&-1/2&\cdots&0&0&0\\ &\vdots&&\ddots&&\vdots&\\ 0&0&0&\cdots&1/2&0&0\\ 0&0&0&\cdots&0&1/2&0\\ 0&0&0&\cdots&-1/2&0&\rho^{-}\\ 1/2&0&0&\cdots&0&-1/2&0\end{pmatrix},}$$

$$\color[rgb]{0,0,0}{\mathbf{C}=\begin{pmatrix}0&1/4&0&\cdots&0&0&\rho^{+}/2\\ 1/4&0&1/4&\cdots&0&0&0\\ 0&1/4&0&\cdots&0&0&0\\ 0&0&1/4&\cdots&0&0&0\\ &\vdots&&\ddots&&\vdots&\\ 0&0&0&\cdots&1/4&0&0\\ 0&0&0&\cdots&0&1/4&0\\ 0&0&0&\cdots&1/4&0&\rho^{-}/2\\ 1/4&0&0&\cdots&0&1/4&0\end{pmatrix}.}$$

$$ϵ\frac{\partial \vec{P}}{\partial \hat{t}}(x,\hat{t})=\mathbf{A}\vec{P}+ϵ\mathbf{B}\frac{\partial \vec{P}}{\partial x}+{ϵ}^2\mathbf{C}\frac{{\partial }^2\vec{P}}{\partial x^2},$$

$${\left[{\vec{u}}_0(ϵ)\right]}_n=2{\rho }^{+}ϵ-2n\left({\rho }^{+}-{\rho }^{-}\right){ϵ}^2,$$

$$\frac{\partial {\vec{P}}_0}{\partial \hat{t}}=\mathbf{A}{\vec{P}}_1+\mathbf{B}\frac{\partial {\vec{P}}_0}{\partial x}.$$

$$\mathbf{A}\vec{P}_{1}=\color[rgb]{0,0,0}{\vec{u}_{0}(\epsilon)}\color[rgb]{0,0,0}{\dfrac{\partial f}{\partial\hat{t}}-\mathbf{B}}\color[rgb]{0,0,0}{\vec{u}_{0}(\epsilon)}\color[rgb]{0,0,0}{\dfrac{\partial f}{\partial x}}.$$

$$\dfrac{\partial f}{\partial\hat{t}}=\dfrac{\vec{v}^{T}\mathbf{B}\color[rgb]{0,0,0}{\vec{u}_{0}(\epsilon)}}{\vec{v}^{T}\color[rgb]{0,0,0}{\vec{u}_{0}(\epsilon)}}\dfrac{\partial f}{\partial x}.$$

$${\left[{\vec{u}}_1(ϵ)\right]}_n={\beta }_0(ϵ)+{\beta }_1(ϵ)n+{\beta }_2(ϵ)n^2+{\beta }_3(ϵ)n^3,$$

$${\beta }_0(ϵ)$$

$${\beta }_1(ϵ)$$

$${\beta }_2(ϵ)$$

$${\beta }_3(ϵ)$$

$${\left[{\vec{u}}_1(ϵ)\right]}_N=ϵ+\frac{{\rho }^{+}-{\rho }^{-}}{6}\left(1-{ϵ}^2\right).$$

$$\frac{\partial {\vec{P}}_1}{\partial \hat{t}}=\mathbf{A}{\vec{P}}_2+\mathbf{B}\frac{\partial {\vec{P}}_1}{\partial x}+\mathbf{C}\frac{{\partial }^2{\vec{P}}_0}{\partial x^2}.$$

$$\mathbf{A}\vec{P}_{2}={\color[rgb]{0,0,0}\vec{u}_{0}(\epsilon)}\dfrac{\partial g}{\partial\hat{t}}-{\color[rgb]{0,0,0}V(\epsilon)\vec{u}_{1}(\epsilon)}\dfrac{\partial^{2}f}{\partial x^{2}}-\mathbf{B}{\color[rgb]{0,0,0}\vec{u}_{1}(\epsilon)}\dfrac{\partial^{2}f}{\partial x^{2}}-\mathbf{C}{\color[rgb]{0,0,0}\vec{u}_{0}(\epsilon)}\dfrac{\partial^{2}f}{\partial x^{2}}-\mathbf{B}{\color[rgb]{0,0,0}\vec{u}_{0}(\epsilon)}\dfrac{\partial g}{\partial x}.$$

$$\dfrac{\partial g}{\partial\hat{t}}+{\color[rgb]{0,0,0}V(\epsilon)}\dfrac{\partial g}{\partial x}=\dfrac{\vec{v}^{T}\mathbf{B}{\color[rgb]{0,0,0}\vec{u}_{1}(\epsilon)}+\vec{v}^{T}\mathbf{C}{\color[rgb]{0,0,0}\vec{u}_{0}(\epsilon)}+{\color[rgb]{0,0,0}V(\epsilon)}\vec{v}^{T}{\color[rgb]{0,0,0}\vec{u}_{1}(\epsilon)}}{\vec{v}^{T}{\color[rgb]{0,0,0}\vec{u}_{0}(\epsilon)}}~{}\dfrac{\partial^{2}f}{\partial x^{2}}.$$

$$\vec{P}=\vec{P}_{0}+\epsilon\vec{P}_{1}+\ldots=(f(x,\hat{t})+\epsilon g(x,\hat{t})){\color[rgb]{0,0,0}\vec{u}_{0}(\epsilon)+\epsilon f_{x}(x,\hat{t})\vec{u}_{1}(\epsilon)}+\ldots.$$

$${\color[rgb]{0,0,0}\epsilon\dfrac{\partial h_{\epsilon}}{\partial t}+V(\epsilon)\dfrac{\partial h_{\epsilon}}{\partial x}=\epsilon D(\epsilon)\left(\dfrac{\partial^{2}h_{\epsilon}}{\partial x^{2}}-\epsilon\dfrac{\partial^{2}g}{\partial x^{2}}\right),}$$

$$\frac{\partial h}{\partial t}+\hat{V}\frac{\partial h}{\partial x}=\hat{D}\frac{{\partial }^{2}h}{\partial x^2},$$

$$\hat{V}$$

$$\hat{D}$$

$$⟨x(t)⟩$$

$$⟨x^2(t)⟩$$

$$h(x,t)=\frac{1}{\sqrt{\pi \hat{D}t}}\text{exp}\left(-\frac{(x-\hat{V}t{)}^2}{4\hat{D}t}\right)+\frac{\hat{V}}{2\hat{D}}\text{exp}\left(-\frac{\hat{V}x}{\hat{D}}\right)\text{erfc}\left(\frac{x+\hat{V}t}{2\sqrt{\hat{D}t}}\right),$$

$$u(x)=2{\rho }^{+}(1-x)+2{\rho }^{-}x,$$

$$⟨x(t)⟩\sim \frac{2\sqrt{t}}{\sqrt{\pi }},⟨x^2(t)⟩\sim t.$$

$$t_{\text{sw}}\sim {\left({\rho }^{+}-{\rho }^{-}\right)}^{-2}={\left(1-\frac{2}{\bar{d}}\right)}^{-2},$$

Full text

1 Introduction

The migration of stem cells during embryogenesis [1]; the transport of proteins along microtubules within eukaryotic cells [2]; and the rotational dynamics of polymers [3] are all examples of biological transport processes that take place in complex environments. Transport processes in such environments are not just of interest within the biological sciences, but are ubiquitous in the study of traffic flow [4, 5] and human crowd management [6]. Many of these transport processes are inherently random, and they can be investigated through studying statistics of stochastic transport models. A critical transport statistic used to classify the character of a transport process is the mean squared displacement (MSD), which at time $t$ is denoted by $\langle x^{2}(t)\rangle$. A particle whose position undergoes Brownian motion and lies on the real line, has MSD that evolves linearly in time, $\langle x^{2}(t)\rangle\propto t$. This is a defining feature of what is referred to as classical diffusion. For more general transport processes the MSD follows a power law, $\langle x^{2}(t)\rangle\propto t^{\alpha}$, and the transport process is referred to as: sub-diffusive if $0<\alpha<1$; diffusive if $\alpha=1$; super-diffusive or sub-ballistic if $1<\alpha<2$; and ballistic if $\alpha=2$.

For simple geometries such as the real line it is unambiguous how to define displacement: setting the origin of the real line to be at the initial position of the particle the displacement is defined as the current position of the particle, $x(t)$. However, for more complex geometries there are several ways to measure displacement. Consider a Brownian particle on a comb lattice [7, 8, 9, 10, 11], a topological structure consisting of an infinite backbone (the real line) with teeth (dendrites) of infinite length that protrude periodically out of the backbone. The MSD, when displacement is measured as the distance traversed along the backbone only, is sub-diffusive with $\langle x^{2}(t)\rangle\propto t^{1/2}$ for long times [10]. The particle spends large amounts of time diffusing along the protruding teeth, during which no transport along the backbone occurs. These long periods of time, where the displacement along the backbone does not change, are responsible for the sub-diffusive nature of the MSD. However, if we define displacement to also include the distance travelled down the teeth then we do not retain the sub-diffusive MSD. The choice of definition for displacement can therefore have a large impact on what we can learn about the transport process.

For topological structures that are embedded within a two- or three-dimensional domain a common definition of displacement is the Euclidean distance between the current and initial positions. With this definition, the potential for geometry to induce a qualitative change in the displacement of a particle has been widely investigated for a class of complex structures known as fractals [12, 13, 14, 15, 16, 17]. In fact, for Brownian particles on these structures, a random walk dimension, $d_{w}$, can be defined such that $\langle x^{2}(t)\rangle\propto t^{2/d_{w}}$ [18, 19]. However, if the topological structure cannot be embedded within a Euclidean domain then this measure of displacement is no longer appropriate. For random walks on general networked topologies (not necessarily embedded in a Euclidean domain), a measure of displacement sometimes used in the literature is the shortest path between the current and initial positions [20]. However, the shortest path does not include information about the actual path the particle took. Statistics of the path of a particle are of interest across research fields such as polymer physics [21], flux lines within superconductors [22], and more broadly in the winding statistics of stochastic processes. Previous studies of winding statistics have focussed on the winding angle of a Brownian particle in a planar domain [23, 24]. However, there has been little research on winding statistics for particles constrained to more complex environments. Kundu et. al [25] have studied winding statistics of a Brownian particle constrained to the unit circle. They calculate the exact distribution of the net number of clockwise or anticlockwise rotations as a function of time. We investigate the winding statistics of particles constrained to more complex networks. We introduce a new definition of displacement, which we term the winding distance, and investigate how it evolves over time as a function of network topology.

The remainder of this paper is organised as follows. In Section 2 we motivate and define our notion of displacement of a particle on a finite network by introducing a topological structure we term the displacement tree. In order to analyse the displacement of a Brownian particle on a finite network at long times, we identify a relationship between displacement along the displacement tree and displacement on the half line with periodic bias. We calculate the strength of this bias as a function of the average degree of the network, a very simple topological property. In Section 3 we study the displacement of a diffusive particle on the half line subject to periodic bias via a multiple scales approach introduced in [26]. We derive a macroscopic advection-diffusion partial differential equation (PDE) to describe the evolution of the distribution of the winding distance of a particle. The diffusive and advective transport coefficients are calculated in terms of topological properties of the network. We uncover a topological condition that enables us to distinguish whether the long-time behaviour of the displacement of a Brownian particle is characteristically diffusive or ballistic, and a prediction of the timescale on which the transition between the two behaviours occurs. In Section 4 we generalise the transport process of the particle from Brownian motion to a continuous time random walk (CTRW) model of anomalous diffusion. Repeating the multiple scales approach in Section 3 we derive a macroscopic fractional advection-diffusion PDE for the evolution of the winding distance of an anomaloulsy diffusive particle on a network. As in Section 3, we obtain a prediction of the timescale on which a particle transitions from the short-time to the long-time behaviour and discuss its asymptotic properties for small values of $\alpha$. Finally, in Section 5 we summarise our results and discuss the scope for future work.

2 Displacement on networks

Consider a network $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$, where $\mathcal{V}$ is the set of vertices and $\mathcal{E}$ the set of edges. Each edge has an associated edge length $L_{e}>0$ for all $e\in\mathcal{E}$. A point particle traverses the edges of the network, and on any given edge its position at time $t$ is denoted by $X(t)$. In order for the position $X(t)$ on the edge $e$ to uniquely identify the location of the particle on the network, each edge is oriented such that the vertex located at one end of the edge has position [math] and the other vertex has position $L_{e}$ (the choice of which vertex has which position is arbitrary). The position of a particle on each edge evolves according to a standard Brownian motion. When a particle reaches a vertex, located at positions [math] or $L_{e}$, an adjacent edge $e_{*}$ is selected uniformly at random111The standard numerical scheme for simulating a Brownian motion is the Euler-Maruyama method. For this scheme a particle performs jumps at discrete time steps of length $\Delta t$ and, as such, a particle will never reach a vertex but jump past it. For a discussion of how to implement the selection of adjacent edges see Appendix A.222The sampling procedure for the adjacent edge is formally incorporated in the definition of a Walsh Brownian motion, for more details see [27].. The new position of the particle is then either $X(t)=0$ or $X(t)=L_{e_{*}}$ depending on how the edge $e_{*}$ is oriented. The initial condition for the particle can be either a position $x_{0}\in(0,L_{e})$ on a known edge $e$ of the network, or the particle can reside initially on a vertex of the network. In the latter case an initial edge is sampled uniformly from all adjacent edges to the initial vertex.

Rather than the position of a particle at time $t$, it is often of interest to consider some measure of how far the particle has travelled up to and including time $t$. Such measures can be considered to be different definitions for the displacement of the particle. To aid discussion of different definitions of displacement, we introduce a simple class of networks. Let $\mathcal{G}_{n}$ be the network comprising two vertices connected by $n$ edges of the same length, $L$. The shortest path definition for displacement on networks sometimes used in the literature [20] is the length of the shortest path between the position at time $t$ and the initial position (at time $t=0$). In certain situations, for example when considering the rotation of polymers around cylindrical fibres [21], we are interested in the winding statistics of stochastic processes in constrained topologies [25]. Therefore, using the shortest path definition for displacement can be an undesirable transport statistic as it does not include any information about the winding of a particle between the initial and current times. Instead, consider a particle attached to an “extendable leash”. In Figure 1(a) we present two trajectories on the network $\mathcal{G}_{3}$. The initial position of the particle is the left-most vertex (Fig. 1(a); unfilled circle), and one end of the leash remains stationary at this vertex at all times. The position of the particle evolves according to a Brownian motion. After selecting an initial edge uniformly at random, the particle moves around the network and the leash extends and contracts accordingly. To emphasise the distinction between the two definitions of displacement, in Figure 1(a) we consider again the two possible trajectories. The first has the final position of the particle on the upmost edge with the leash hooked around the right-most vertex (Figure 1(a); solid line with label one). The second trajectory has the same final position as the first. However, in this case, the particle travelled directly along the upmost edge (Figure 1(a); dashed line with label two) and consequently the leash corresponding to the second trajectory is shorter. The length of the leash we refer to as the winding distance of the particle and is the definition of displacement we consider in this work.

For the circle network, $\mathcal{G}_{2}$, the winding distance of a particle is equivalent to the displacement of a particle whose position evolves according to a Brownian motion on the infinite one-dimensional line [25, 28]. Therefore, the winding distance, $x(t)$, of a Brownian particle at time $t$ on the unit circle is well studied, and the MSD is $\langle x^{2}(t)\rangle=2Dt$ for all times $t$, where $D$ is the diffusion coefficient. To investigate the evolution of the winding distance of a particle on the network $\mathcal{G}_{3}$ (Figure 1(a)), we perform numerical realisations of a Brownian motion. From $10^{4}$ sample paths an estimate for the MSD is calculated and Figure 1(b) shows the evolution of $\langle x^{2}(t)\rangle/t$ on a logarithmic scale. Initially $\langle x^{2}(t)\rangle=2Dt$, as on short timescales the particle has not left the initial edge and so is classically diffusive. However as time increases, the solid (red) curve in Figure 1(b) deviates away from the line $\langle x^{2}(t)\rangle=2Dt$ (Figure 1(b); dashed line) and in the long-time limit settles to $\langle x^{2}(t)\rangle\propto t^{2}$. This characteristic change in the MSD from diffusive to ballistic behaviour can be explained by considering the behaviour of a particle at the vertices. The first trajectory in Figure 1(a) (solid line) reaches the right-most vertex via the bottom edge. At this vertex the particle then selects a new edge uniformly at random. If the new edge is either the top or middle edge (which occurs with total probability $2/3$), then moving along the new edge towards the left-most vertex will increase the displacement. However, if the new edge is the bottom edge (which occurs with probability $1/3$), then the displacement will decrease if the particle travels (back) along the bottom edge. Consider temporarily, a particle which hops directly between the vertices in $\mathcal{G}_{3}$ over an edge, selected uniformly at random. The winding distance is now a stochastic process, which increases by $L$ if the new edge that the particle hops over is distinct from the previous edge, or decreases by $L$ if the new edge is the same as the previous edge. Thus the winding distance evolves as a biased random walk on $\{0,L,2L,\ldots\}$ with probability $2/3$ to jump in the positive direction. This is a well-studied random walk that is known to have a ballisitic MSD for all times $t$. Thus, the introduction of continuous stochastic motion along the edges of the network is responsible for the transition from diffusive to ballistic behaviour.

For diffusive transport processes, can we identify which topologies will cause the displacement of a particle to become characteristically ballistic? Given a topology that will induce ballistic behaviour, can we predict the timescale on which this transition will occur and how it depends on properties of the network? How does topology affect transport statistics for other transport processes? In order to investigate these questions we first introduce a useful topological structure.

2.1 The displacement tree

Recall that, for the circle $\mathcal{G}_{2}$, the winding distance of a particle is equivalent to the displacement of a particle on the real line. A similar relationship can be identified between the winding distance of a particle on a network and the displacement of a particle on a new topological structure we term the displacement tree. To construct a displacement tree, $\mathcal{T}_{3}$, for the network $\mathcal{G}_{3}$ we need to select an initial position for the particle, which we take as the left-most vertex (see Figure 2; unfilled circles). The three choices for the initial edge result in three potential initial trajectories for the particle. Consequently the tree $\mathcal{T}_{3}$ starts with its root, corresponding to the initial vertex, from which three parallel trajectories protrude (see Figure 2). Now consider a particle that has just reached the right-most vertex of $\mathcal{G}_{3}$; the particle must choose a new edge. If it chooses the edge it just traversed and travels back along that edge the displacement of the particle will decrease. Travelling along the other two edges will result in an increase in displacement. To represent these possibilities, each of the three initial trajectories of the displacement tree branch into two new trajectories for the two edge choices that if travelled along will increase displacement. Repeating this procedure produces an infinite tree that spans the entire ensemble of possible trajectories for a particle on $\mathcal{G}_{3}$, oriented such that trajectories that increase displacement lead to the right. The winding distance of a Brownian particle on $\mathcal{G}_{3}$ is therefore equivalent to how far a Brownian particle has travelled along $\mathcal{T}_{3}$ (Figure 2).

The displacement tree $\mathcal{T}_{3}$ has a very simple structure: as the edges of $\mathcal{G}_{3}$ are all of equal length, $L$, the branching points appear on the tree periodically along all trajectories at positions $kL$ for integers $k\geq 1$. We define a protruding edge on the displacement tree to be an edge that emerges from a branching point and leads to the right. Both the vertices in $\mathcal{G}_{3}$ have degree three and so the number of protruding edges from any of the branching points in the displacement tree is the same, and equal to two. We refer to these displacement trees as regular: every possible trajectory of a particle on $\mathcal{T}_{3}$ of some fixed displacement will have passed the same number of branching points, all of which had the same number of protruding edges. As a result, we can equate the displacement of a particle on such a regular displacement tree with the displacement of a particle undergoing a new transport process on the non-negative half line.

Consider a particle whose position on the regular displacement tree, $X(t)$, lies within the interval $[(k-1)L,kL]$ for some positive integer $k$. The position of the particle evolves according to a Brownian motion unless the particle is positioned on one of the boundaries of the interval, $X(t)=kL$ or $X(t)=(k-1)L$. These boundary positions correspond to the branching points on the displacement tree and, as such, a particle moves to the right of these points with probability $2/3$ and to the left with probability $1/3$. It is equivalent to just consider a particle on the non-negative half line, where the particle moves between the intervals $[(k-1)L,kL]$ for integers $k\geq 1$. For a particle with position $X(t)=kL$, the probability $2/3$ to move to the right is incorporated as a Robin boundary condition. The particle is absorbed into the new interval $[kL,(k+1)L]$ (a move to the right along the displacement tree) with probability $2/3$ and reflected back into the interval $[(k-1)L,kL]$ (a move to the left along the displacement tree) with probability $1/3$. For a particle positioned at $X(t)=(k-1)L$, the particle is absorbed into the new interval $[(k-2)L,(k-1)L]$ (a move to the left along the displacement tree) with probability $1/3$ and reflected back into the interval $[(k-1)L,kL]$ (a move to the right along the displacement tree) with probability $2/3$. This transport process is a Brownian motion on the non-negative half line with singular periodic bias of strength $2/3$.

2.2 Extending to general networked topologies

Networks that represent complex environments will not necessarily have regular displacement trees. For a general network, $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$, the edge lengths, $L_{e}$ for $e\in\mathcal{E}$, are not always of equal length and the vertices will not always have equal degree. An example of such a network is shown in Figure 3(a) where two vertices are connected by three multiple edges, similarly to $\mathcal{G}_{3}$, but now the top, middle and bottom edges have edge lengths $3L$, $L$ and $2L$, respectively. Consider two possible trajectories for a particle starting at the left-most node; trajectory one is highlighted with a solid line in Figure 3(b), has a length of $4L$, and finishes at a vertex on the network $\mathcal{G}$, i.e. a branching point on the displacement tree (see Figure 3(c); solid line with label one). Trajectory two is highlighted with a dashed line in Figure 3, also has a length of $4L$ but finishes along an edge, not at a branching point. Therefore knowledge of the displacement of a particle is not sufficient to determine whether the particle is at a branching point on the displacement tree. In order to circumvent this problem we introduce a topologically equivalent network $\mathcal{G}_{L}=\{\mathcal{V}_{L},\mathcal{E}_{L}\}$ which is formed by introducing additional vertices into $\mathcal{G}$ that occur with equal-spacing of length $L$ (see Figure 3(a); black squares). The corresponding displacement tree for $\mathcal{G}_{L}$ has branching points that occur periodically at positions $kL$ for integers $k\geq 1$ along every trajectory (see Figure 3(c); black squares). Note that all the additional branching points introduced will have only one protruding edge as the additional vertices in $\mathcal{G}_{L}$ have degree two. As a result, they will not introduce a bias. For general networks $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$, as long as the edge lengths, $L_{e}$, are rational, or can be approximated as such, there will always exist $L>0$ such that $L_{e}=m_{e}L$ for some positive integer $m_{e}$ for all $e\in\mathcal{E}$, and so the transformation from $\mathcal{G}$ to $\mathcal{G}_{L}$ can be applied.

As for $\mathcal{G}_{3}$, for these more general networks $\mathcal{G}_{L}$, it would be convenient to equate the displacement of a particle diffusing on the displacement tree to the displacement of a particle on the half line. The displacement tree for $\mathcal{G}_{L}$ has branching points along all trajectories at positions $kL$ for integers $k\geq 1$. However, the number of protruding edges is not constant across trajectories with the same length. Consider again the highlighted trajectories of length $4L$ in Figure 3(c). The branching point at the end of trajectory one has two protruding edges. A particle at this branching point travels along an edge that will increase displacement with probability $2/3$. However, trajectory two ends at a branching point with only one protruding edge. The probability of a particle travelling along an edge that will increase displacement is $1/2$. To make progress, we consider a particle on the displacement tree whose position lies in the interval $[(k-1)L,kL]$. For a particle at the boundary position $X(t)=kL$ we seek an averaged probability, $\rho^{+}_{k}$, to select an edge on the displacement tree which lies to the right of the branching points at $kL$. The probability of moving further along the displacement tree if the particle is at vertex $\nu\in\mathcal{V}_{L}$ is $(d_{\nu}-1)/d_{\nu}$ where $d_{\nu}$ is the degree of the vertex $\nu$. Weighting this by $p_{k}(\nu)$, the probability that the particle occupies vertex $\nu$ given it has displacement $kL$, and summing over $\nu\in\mathcal{V}_{L}$, yields

[TABLE]

The probabilities, $\rho^{+}_{k}$, are used to define the Robin boundary conditions for the transport process on the half line, just as for symmetric displacement trees. Consider a particle in the interval $[(k-1)L,kL]$ with position at the boundary $X(t)=kL$. The particle is absorbed into the adjacent interval $[kL,(k+1)L]$ (a move to the right along the asymmetric displacement tree) with probability $\rho^{+}_{k}$. The particle is reflected back into the interval $[(k-1)L,kL]$ (a move to the left along the asymmetric displacement tree) with probability $1-\rho^{+}_{k}$. Therefore, Equation 1 provides the sequence of probabilities, $\{\rho^{+}_{1},\rho^{+}_{2},\ldots\}$, to be absorbed at position $kL$ into the interval $[kL,(k+1)L]$ (see Figure 3(e)). The probabilities $p_{k}(\nu)$ are difficult to calculate as they depend heavily on the topology of the network and the initial position of the particle. However, the displacement of a particle transitions to ballistic behaviour when the periodic bias felt at positions $kL$ dominates over the unbiased Brownian motion. As such, we consider the displacement of a particle in the long-time limit, which is when the particle has traversed sufficiently many edges for the bias to dominate. The branching points in the displacement tree (corresponding to vertices in the network with degree greater than two) introduce an asymmetry, biasing particles further along the tree (to the right). As a particle moves further along the tree we are interested in the limit to which the sequence of effective biases $\{\rho^{+}_{1},\rho^{+}_{2},\ldots\}$ converges. Assume temporarily that the networks we consider are both finite and aperiodic. Then, as $k\rightarrow\infty$ the probability $p_{k}(\nu)$ converges to the equilibrium distribution for a discrete random walk on the vertices of $\mathcal{G}_{L}$ which is $p_{\infty}(\nu)=d_{\nu}/\sum_{\omega\in\mathcal{V}_{L}}d_{\omega}$ 333For general finite Markov Chains an equilibrium distribution for the occupany probability of each state in the chain exists if and only if the chain (here the network) is irreducible and aperiodic. All connected undirected networks however are necessarily irreducible.. The sequence $\{\rho^{+}_{1},\rho^{+}_{2},\ldots\}$ therefore converges to the limit

[TABLE]

where $\bar{d}=2|\mathcal{E}_{L}|/|\mathcal{V}_{L}|$ is the average degree of $\mathcal{G}_{L}$. Equation 2 provides an approximation for the probability of a particle to be absorbed at $kL$ into the interval $[kL,(k+1)L]$ once the particle is sufficiently far along the half line, that is, for large integers $k$. Note that the assumption that the network is finite is not a necessary one. Some infinite networks, such as a honeycomb lattice, have displacement trees where the equilibrium distribution $p_{\infty}(\nu)$ does exist. We can also relax the assumption that the network is aperiodic; see Appendix B.

In summary, in the long-time limit we have identified a link between the displacement of a diffusing particle on a general network and displacement on the half line with a constant periodic bias. The focus of the next section is to use multiple scales analysis to explore the transport properties of the process on the half line.

3 Periodic bias on the half line: multiple scales analysis

In this section we consider the transport properties of a Brownian motion on the half line with periodic bias using an approach from a recent paper by Chapman and Shabala [26]. The paper introduces a method of multiple scales to derive the macroscopic transport equations of a random walk on a periodic lattice with spatially dependent transition rates. The slow scale is continuous and evolves as the particle walks over many periods of the lattice. The fast scale is discrete and is defined on a unit periodic interval made up of $N$ lattice sites. In contrast to [26], in this work we are interested in the transport of a Brownian particle and so we take the limit $N\rightarrow\infty$ to obtain a continuous expression for the fast scale as well as the slow scale.

Brownian motion is defined in continuous space, however to make progress we formulate a continuous time Markov Chain (CTMC) in discrete space that has the same macroscopic properties. We non-dimensionalise in space by setting the length of each interval to be one ($L=1$), the intervals $[k,k+1]$ are then discretised into $N+1$ lattice sites each separated by distance $\epsilon=1/N$ (see Fig. 4). Consider the position of a particle on the discretised lattice that evolves according to a CTMC, where the transition rates to leave each site are symmetric and chosen to be $\lambda_{n\rightarrow n-1}=\lambda_{n\rightarrow n+1}=D/\epsilon^{2}$. In the limit $\epsilon\rightarrow 0$ and $N\rightarrow\infty$, the distribution of the position of a particle which evolves according to this CTMC is equivalent to the same distribution arising from a Brownian motion on the half line with diffusion coefficient $D$. For the remainder of this work we take $D=1/2$. To include the periodic bias at branching points, which occur at lattice sites with indices $kN$ for all integers $k>0$, we choose asymmetric transition rates at these points. For a particle to exit on the right of the branching point with probability $\rho^{+}_{\infty}$ we choose transition rates $\lambda_{kN\rightarrow kN+1}=\rho^{+}/\epsilon^{2}$ and $\lambda_{kN\rightarrow kN-1}=\rho^{-}/\epsilon^{2}$, where $\rho^{-}=1-\rho^{+}$ and we have dropped the $\infty$ subscripts from Section 2. For a diagrammatic representation of the CTMC see Figure 4.

Let $p_{n}(t)$ be the probability a particle occupies lattice site $n\geq 0$ at time $t$. The master equation for $p_{n}(t)$ is

[TABLE]

Initial conditions imply $p_{0}(0)=1$ and $p_{i}(0)=0$ for all $i\geq 1$ as the initial winding distance is always zero. The winding distance is defined such that it is always non-negative, therefore we impose a reflective boundary condition at the origin by setting $\lambda_{-1\rightarrow 0}=\lambda_{0\rightarrow-1}=0$ and $\lambda_{0\rightarrow 1}=1/\epsilon^{2}$. We define the spatial coordinate $x=\epsilon n$, and note that, for $\epsilon\ll 1$, $x$ evolves on a slower time scale than the lattice index $n$. Following [26], let $p_{n}(t)=P_{n}(x,t)$ so the probability density for the position of a particle at time $t$ depends on both the fast scale, $n$, and the slow scale $x$. Re-writing Eq. 3 in terms of $P_{n}(x,t)$, we have

[TABLE]

for $n\geq 0$ and $x\in[0,\infty)$. To derive a macroscopic PDE in the spatial variable $x$, we assume that the variables $x$ and $n$ are independent. This assumption reduces the infinite set of equations (3) to a finite set of $N$ distinct equations. Introducing circular notation $\overline{i\pm 1}=\left({i\pm 1}-1\text{ mod }N\right)+1$, we write the master equation for the probabilities $P_{i}(x,t)$ as

[TABLE]

for $i\in\{1,\ldots,N\}$ and $x\in[0,\infty)$. The finite set of equations defines the periodic unit interval. On the periodic unit interval there are symmetric transition rates of $1/\left(2\epsilon^{2}\right)$ for all adjacent sites, other than the rates $\lambda_{N\rightarrow 1}=\rho^{+}/\epsilon^{2}$ and $\lambda_{N\rightarrow N-1}=\rho^{-}/\epsilon^{2}$.

Following the derivation in [26], we introduce $\vec{P}(x,t)=\left[P_{1}(x,t),\ldots,P_{N}(x,t)\right]^{T}$ and expand about $x$ to obtain

[TABLE]

where

[TABLE]

Note that when $\rho^{+}=\rho^{-}=1/2$ we have symmetric transition rates of $1/(2\epsilon^{2})$ at all sites on the unit interval. These rates were chosen such that the time derivative on the right-hand side of Eq. 6 appears at the same order in $\epsilon$ as the second order spatial derivatives on the left-hand side. Therefore, as $\epsilon\rightarrow 0$, the probability distribution for the displacement of a particle follows a macroscopic diffusion equation. However, when $\rho^{+}\neq\rho^{-}$, recall that the displacement of a particle on a networked topology may transition from diffusive to ballistic behaviour. As such, we anticipate both advective and diffusive contributions to the macroscopic PDE describing the displacement of a particle and we rescale time, $\hat{t}=\epsilon^{-1}t$, to give

[TABLE]

where the time-dependent derivative now appears at the same order of magnitude as the first order spatial derivative. We make the series expansion $\vec{P}(x,\hat{t})=\vec{P}_{0}(x,\hat{t})+\epsilon\vec{P}_{1}(x,\hat{t})+\epsilon^{2}\vec{P}_{2}(x,\hat{t})+\ldots$ and substitute into Eq. 8. Collecting terms of $\mathcal{O}(\epsilon^{0})$ yields $\mathbf{A}\vec{P}_{0}=\vec{0}$, which has solution \vec{P}_{0}=f(x,\hat{t})\color[rgb]{0,0,0}{\vec{u}_{0}(\epsilon)} where $f(x,\hat{t})$ is some function to be determined and \color[rgb]{0,0,0}{\vec{u}_{0}}(\epsilon) is in the kernel of $\mathbf{A}$444The matrix $\mathbf{A}$ is the rate matrix for the CTMC on the periodic interval, therefore as $\vec{u}_{0}(\epsilon)$ is in the kernel of $\mathbf{A}$ it is proportional to the equilibrium distribution of the CTMC. which is given by

[TABLE]

for $1\leq n\leq N-1$ and $\left[\vec{u}_{0}(\epsilon)\right]_{N}=\epsilon$, such that $\sum_{n=1}^{N}\left[\vec{u}_{0}(\epsilon)\right]_{n}=1$. Collecting terms of $\mathcal{O}\left(\epsilon^{1}\right)$ we find

[TABLE]

Substituting in the expression for $\vec{P}_{0}$ and rearranging yields

[TABLE]

Applying the Fredholm alternative theorem555The Fredholm alternative theorem states that $\mathbf{M}\vec{x}=\vec{b}$ has a solution if and only if for all $\vec{y}$ such that $\mathbf{M}^{T}\vec{y}=\vec{0}$ we also have $\vec{y}^{T}\vec{b}=\vec{0}$. we see that, in order for Eq. 11 to have a solution, we need the right-hand side to be orthogonal to the null space of $\mathbf{A}^{T}$. The null space is spanned by the vector $\vec{v}=[1,\ldots,1]^{T}$, and the Fredholm alternative theorem therefore implies

[TABLE]

Introducing $V(\epsilon)=-\vec{v}^{T}\mathbf{B}\vec{u}_{0}(\epsilon)/\vec{v}^{T}\vec{u}_{0}(\epsilon)$ we can solve Eq. 11 to give $\vec{P}_{1}=f_{x}\vec{u}_{1}(\epsilon)+g\vec{u}_{0}(\epsilon)$, where $g(x,\hat{t})$ is some function to be determined, $\vec{u}_{1}(\epsilon)$ is such that $\mathbf{A}\vec{u}_{1}(\epsilon)=-\left(V(\epsilon)\mathbf{I}+\mathbf{B}\right)\vec{u}_{0}(\epsilon)$ and $\mathbf{I}$ is the identity matrix. Noting that $V(\epsilon)=\epsilon\left(\rho^{+}-\rho^{-}\right)$ we can solve for $\vec{u}_{1}(\epsilon)$ to obtain

[TABLE]

for $1\leq n\leq N-1$ where the coefficients are given by

[TABLE]

and the final term in $\vec{u}_{1}(\epsilon)$ is given by

[TABLE]

We now proceed to collect terms of $\mathcal{O}\left(\epsilon^{2}\right)$ which gives

[TABLE]

As before, we rearrange Eq. 16 to give

[TABLE]

Applying the Fredholm alternative theorem once again implies

[TABLE]

We can now rewrite the power series as

[TABLE]

Introducing $h_{\epsilon}(x,\hat{t})=f(x,\hat{t})+\epsilon g(x,\hat{t})$, rescaling into time coordinates $t=\epsilon\hat{t}$ and using Eq. 12 and Eq. 18 we have

[TABLE]

where {\color[rgb]{0,0,0}D(\epsilon)=\left(\vec{v}^{T}\mathbf{B}\vec{u}_{1}(\epsilon)+\vec{v}^{T}\mathbf{C}\vec{u}_{0}(\epsilon)+V(\epsilon)\vec{v}^{T}\vec{u}_{1}(\epsilon)\right)/\vec{v}^{T}\vec{u}_{0}(\epsilon)}. Upon dividing Eq. 20 by $\epsilon$, taking the limit $\epsilon\rightarrow 0$ and $N\rightarrow\infty$ where $\epsilon=1/N$, and introducing $h(x,t)=\text{lim}_{\epsilon\rightarrow 0}\left\{h_{\epsilon}(x,t)\right\}$, $\hat{V}=\text{lim}_{\epsilon\rightarrow 0}\left\{V(\epsilon)/\epsilon\right\}$, and $\hat{D}=\text{lim}_{\epsilon\rightarrow 0}\left\{D(\epsilon)\right\}$ we have the following linear advection-diffusion PDE

[TABLE]

for $x\in[0,\infty)$, with boundary conditions $\hat{D}h_{x}(0,t)+\hat{V}h(0,t)=\lim_{x\rightarrow\infty}\{\hat{D}h_{x}(x,t)+\hat{V}h(x,t)\}=0$ and initial conditions $h(x,0)=\delta(x)$. Equation 21 is a macroscopic PDE for the evolution of the displacement (defined as the winding distance) of a particle on a network $\mathcal{G}$, valid in the long-time limit. Note that the scalar coefficient of $\vec{u}_{1}(\epsilon)$ in Eq. 19 is $\epsilon f_{x}(x,\hat{t})$. Upon rescaling time $t=\epsilon\hat{t}$, it is clear from Eq. 12 that $f_{x}(x,t)$ is independent of $\epsilon$, and we find that $\lim_{\epsilon\rightarrow 0}\left\{\epsilon f_{x}(x,t)\right\}=0$. The entries of $\vec{u}_{1}(0)$ given in Eq. 15 are finite. Therefore, as $\epsilon\rightarrow 0$, there is no contribution to $\vec{P}(x,t)$ from the term $\epsilon f_{x}(x,t)\vec{u}_{1}(\epsilon)$. The terms $\hat{V}h_{x}$ and $\hat{D}h_{xx}$ in Eq. 21 represent the effective drift and diffusion contributions, respectively, to the distribution of the displacement of a particle. The transport coefficients $\hat{V}$ and $\hat{D}$ determine the strengths of these contributions as a function of the networked topology. We use the expressions for $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ in Eq. 7, as well as $\vec{u}_{0}(\epsilon)$ and $\vec{u}_{1}(\epsilon)$ in Eq. 9 and Eq. 15, respectively, to calculate the effective transport coefficients as

[TABLE]

From Eq. 21 we calculate expressions for the first two moments of the displacement as

[TABLE]

where $\bar{d}$ is the average degree of the network. The ballistic term $\hat{V}^{2}t^{2}$ in (23b) is nonzero when $\bar{d}>2$, that is, the average degree of $\mathcal{G}_{L}$ must be strictly greater than two for the network to induce ballistic behaviour in the winding distance of a diffusive particle.

3.1 Full distribution of displacement

We are not in fact limited to approximating the first two moments of the displacement but can also obtain an analytical approximation to the full distribution. Recall that the distribution of the displacement of a particle whose position evolves according to the CTMC is equivalent to the same distribution as that of a Brownian particle in the limit $\epsilon\rightarrow 0$ and $N\rightarrow\infty$ such that $\epsilon=1/N$. Note that $\vec{P}(x,t)\approx h_{\epsilon}(x,t)\vec{u}_{0}(\epsilon)$, where $h(x,t)=\lim_{\epsilon\rightarrow 0}\left\{h_{\epsilon}(x,t)\right\}$ satisfies Eq. 21 which has solution

[TABLE]

where $\text{erfc}(\cdot)$ is the complimentary error function [29], and $\vec{u}_{0}(\epsilon)$ is given in Eq. 9. We introduce the function $u(x)=\lim_{\epsilon\rightarrow 0}\{\epsilon^{-1}\left[\vec{u}_{0}(\epsilon)\right]_{n}\}$, where $x=\epsilon n$ and find that

[TABLE]

for $x\in[0,1)$ and $u(1)=1$. Extending the function $u(x)$ to be periodic on the half line with unit period, $u(x+1)=u(x)$ for all $x\in(0,\infty)$, provides us with a function that accounts for corrections on the microscopic scale to the displacement predicted from the macroscopic equation for $h(x,t)$. Then the full distribution is approximated by $P(x,t)\approx h(x,t)u(x)$. For $\mathcal{G}_{10}$, we simulate $10^{6}$ realisations of a Brownian particle terminating at time $t=50$. We use these realisations to empirically plot the distribution of displacement (Fig. 5; red solid curve) and compare with the analytical approximation using the solution $h(x,t)$ in Eq. 24 and the periodically extended function $u(x)$ in Eq. 25 (Fig. 5; black dashed curve). The two curves match well which suggests that the temination time $t=50$ is sufficiently large for the transport process to have reached its equilibrium behaviour for the network $\mathcal{G}_{10}$. However, in general how does the topology of the network affect the time taken to reach equilibrium?

3.2 Timescale for transition to equilibrium

From Eq. 23 we already have the dominant contributions to the first two moments of displacement for long times. Now consider short times, $t$, sufficiently short such that the particle has a displacement less than one, or equivalently has not reached the first branching point on the displacement tree. On this timescale, the particle undergoes unbiased Brownian motion on the half line, for which the first two moments of displacement are

[TABLE]

The dominant terms of the long-time and short-time analytical asymptotics in Eq. 23 and Eq. 26, respectively, are plotted in Figure 6 (dashed and dot-dashed lines). Equating the dominant contributions to the MSD on both the long and short timescales, $\hat{V}^{2}t^{2}$ and $t$, respectively, we obtain a prediction of the timescale, $t_{\text{sw}}$, upon which the displacement of a particle transitions from characteristically diffusive to ballistic:

[TABLE]

where $\bar{d}=2|\mathcal{E}_{L}|/|\mathcal{V}_{L}|$ is the average degree of the network. The timescale, $t_{\text{sw}}$, is given by the intersection of the dashed and dot-dashed lines in Figure 6(b) (vertical line).

To investigate the predictive capacity of the timescale $t_{\text{sw}}$, we consider the MSD over an ensemble of randomly generated aperiodic networks. The MSD of a particle on a general network has dominant contribution $\langle x^{2}(t)\rangle\sim\hat{V}^{2}t^{2}$ for long times. The rescaling $\bar{t}=t/t_{\text{sw}}$ allows for the MSD of a particle over an ensemble of topologies to be compared on the same temporal axis, as $\langle x^{2}(t)\rangle/t\sim\bar{t}$. In Figure 7(a) the MSD of a particle over an ensemble of $100$ networks is presented. The $100$ networks each have $10$ vertices and a range of $11$ to $110$ edges666The sampling procedure is as follows; select a desired number of edges and vertices; draw two distinct vertices uniformly at random and join them via an edge; repeat the previous step until the required number of edges is reached; check if the network is periodic or disconnected, if so reject and start over, if not take the network as a realisation. This procedure allows for multiple edges but avoids self loops.. For each network, the particle starts at the same vertex for all realisations, but that fixed vertex is randomly selected and highlighted in Figure 7(c) with a star777The initial vertex is fixed among realisations of the transport process, so that the displacement tree has a fixed root and the effects of the initial conditions on transition times are not averaged out..

On visual inspection of the data in Figure 7 we see that the grey curves rarely intersect as they collapse onto the curve $\langle x^{2}(t)\rangle/t=\bar{t}$. As such, the curves that intersect the line $t=t_{\text{sw}}$ (see Figure 7(a) inset) with a lower $y$ intercept will transition to ballistic behaviour before those curves with a greater $y$ intercept. Thus, we plot the $y$ intercept, $\langle x^{2}(t_{\text{sw}})\rangle/t_{\text{sw}}$, for all $100$ networks against the average degree of the network in Figure 7(b). We see a sharper correlation between average degree and the $y$ intercept for networks with a higher average degree. This suggests that $t_{\text{sw}}$ is a good prediction of the crossover time to ballistic motion when the average degree is large (Figure 7: the network highlighted in blue). However, for networks with a lower average degree the predictive capacity of $t_{\text{sw}}$ is less clear; we highlight networks in magenta888The magenta network has an average degree $\bar{d}=2.2$ however the fixed initial condition is the vertex of degree two highlighted with a star in Figure 7(c). The environment local to the initial vertex has a significantly higher average degree than the entire network, which results in a curve that deviates from the short-time behaviour earlier (see the magenta curve in Figure 7(a)). The dip in the gradient of the magenta curve occurs later, when the particle has explored the global environment, which has a lower average degree than the local environment. and yellow in Figure 7 that both have low average degrees yet the values of $\langle x^{2}(t_{\text{sw}})\rangle/t_{\text{sw}}$ are very different (Figure 7(b)). For topologies with a low average degree the timescale to transition from diffusive to ballistic behaviour is highly sensitive to the particular network and, as such, $t_{\text{sw}}$ may be less reliable as a predictor of the crossover to ballistic motion.

4 Anomalous diffusion on networks

Thus far, the only transport process we have considered is Brownian motion. An alternative modelling framework for particle transport is a position jump process, where the position of a particle is updated instantaneously after random intervals in time. A position jump process can be formally described as a continuous time random walk (CTRW) where the jump distance is sampled from $\lambda(x)$, the jump distribution, and the time between jump events is sampled from $\omega(t)$, the waiting time distribution. For a review of transport properties of these CTRWs see [30]. The freedom to prescribe any jump and waiting time distributions allows a CTRW to describe a broad range of physical processes. Consider a waiting time distribution that for long times has asymptotic power law behaviour, $\omega(t)\sim t^{-\left(1+\alpha\right)}$ for some $\alpha$. If $\lambda(x)$ has zero mean and finite variance then for $\alpha\geq 1$, the position of a particle evolving according to the CTRW on the real line in the long-time limit is equivalent to the position of a Brownian particle and hence has a MSD of $\langle x^{2}(t)\rangle\propto t$. However, for $0<\alpha<1$, $\omega(t)$ is heavy tailed and the mean waiting time between jump events diverges. Such transport processes for large times have a MSD on the real line of $\langle x^{2}(t)\rangle\propto t^{\alpha}$, and are known as anomalously diffusive.

Using the displacement tree (Section 2), the winding distance of anomalously diffusive transport on networks can be identified with displacement of the same transport process on the half line subject to periodic bias. As before, we consider a multiple scales approach and discretise the half line into intervals, $[k,k+1]$ for integers $k\geq 0$, each of $N+1$ lattice sites separated by a distance $\epsilon=1/N$. The position of a particle on the lattice evolves according to a CTRW rather than a CTMC because the waiting times for jumps between adjacent sites are no longer exponentially distributed. Similar to Section 3 we now select the fractional transition rates999The mean time for a jump to occur on the lattice diverges, as such we do not consider transition rates as we would for a CTMC. Instead we introduce a fractional transition rate [31]. to represent the periodic bias felt at the ends of the intervals $[k,k+1]$. Thus, for a particle to exit a lattice site with the index $kN$ to the right with probability $\rho^{+}$ we select the transition rates $\mu_{kN\rightarrow kN+1}=\rho^{+}/\epsilon^{2}$ and $\mu_{kN\rightarrow kN-1}=\rho^{-}/\epsilon^{2}$, and all other lattice sites $n$ have symmetric transition rates $\mu_{n\rightarrow n+1}=\mu_{n\rightarrow n-1}=1/\left(2\epsilon^{2}\right)$.

Let $p_{n}(t)$ be the probability a particle occupies lattice site $n$ at time $t$. The temporal evolution of the probability $p_{n}(t)$ evolves according to the fractional master equation [31], a generalisation of the master equation seen in Section 3. We define the operator $\partial^{\alpha}/\partial t^{\alpha}$ as the Riemann-Liouville fractional derivative

[TABLE]

The fractional master equation for the probabilities $p_{n}(t)$ is

[TABLE]

where $n\geq 0$. As before, the particle is initially at the origin, therefore $p_{0}(0)=1$ and $p_{i}(0)=0$ for all $i\geq 1$. To account for the reflective boundary condition at the origin, we set $\mu_{-1\rightarrow 0}=\mu_{0\rightarrow-1}=0$ and $\mu_{0\rightarrow 1}=1/\epsilon^{2}$. Introducing the spatial coordinate $x=\epsilon N$ where $\epsilon\ll 1$, we write $p_{n}(t)=P_{n}(x,t)$ so that the probability density for the position of a particle at time $t$ depends on both the fast scale $n$ and the slow scale $x$. Treating the two variables $x$ and $n$ independently and rewriting the fractional master equation Eq. 29 in terms of $P_{n}(x,t)$ gives the finite set of equations

[TABLE]

for $i\in\{1,\ldots,N\}$ and $x\in[0,\infty)$, where $\overline{i\pm 1}=\left({i\pm 1}-1\text{ mod }N\right)+1$. For the periodic interval, as before the transition rates out of all sites are symmetric and equal to $1/\left(2\epsilon^{2}\right)$, other than the rates $\mu_{N\rightarrow 1}=\rho^{+}/\left(\epsilon^{2}\right)$ and $\mu_{N\rightarrow N-1}=\rho^{-}/\left(\epsilon^{2}\right)$. Introducing $\vec{P}(x,t)=\left[P_{1}(x,t),\ldots,P_{N}(x,t)\right]^{T}$ and expanding the right-hand side of Eq. 30 about $x$, we find

[TABLE]

where the matrices are as given in Eq. 7. Let $h(x,t)$ denote the macroscopic probability density function for a particle to have a winding distance of $x$ at time $t$. The analysis proceeds identically to Section 3 and we derive the following fractional advection-diffusion PDE

[TABLE]

for $x\in[0,\infty)$, with boundary conditions $\hat{D}h_{x}(0,t)+\hat{V}h(0,t)=\lim_{x\rightarrow\infty}\{\hat{D}h_{x}(x,t)+\hat{V}h(x,t)\}=0$ and initial conditions $h(x,0)=\delta(x)$. The transport coefficients $\hat{V}$ and $\hat{D}$ are given in Eq. 22. This macroscopic fractional PDE is valid for large times and the first two moments of displacement are calculated to be

[TABLE]

Equation 33 demonstrates that anomalous diffusion on a network can induce a variety of transitional behaviours. Consider the analytical expression in the long-time limit for the MSD in Eq. 33b. At large times the dominant contribution is proportional to $t^{2\alpha}$ when $\bar{d}>2$. If $0<\alpha<1/2$, the winding distance of the particle in the long-time limit remains sub-diffusive but with an increased exponent of $2\alpha$. If $\alpha=1/2$, the winding distance becomes classically diffusive. For $1/2<\alpha<1$ the winding distance transitions to sub-ballistic or super-diffusive. Generalising to anomalously diffusive transport processes demonstrates that in order to capture the effect of network topology on the winding distance of a particle one must also incorporate details of the transport process itself.

As in Section 3 we shall compare the first two moments of displacement on both long and short timescales. On short timescales, where the particle has not left the first edge, the displacement of a particle is equivalent to the displacement of a particle whose position is described by a fractional diffusion equation on the half line. The first two moments of displacement [32] on short timescales are

[TABLE]

We compare the dominant short-time and long-time behaviours for the MSD to estimate the transition timescale

[TABLE]

where for $\alpha=1$ we recover Eq. 27.

Consider the limit where $\alpha\rightarrow 0$, which corresponds to a transport process where the waiting time distribution between jump events becomes increasingly heavy tailed. Expanding Eq. 35 for small $\alpha$ yields

[TABLE]

where $\gamma$ is the Euler-Mascheroni constant. In the limit $\alpha\rightarrow 0$ there are two cases for the limiting behaviour of the timescale $t_{\text{sw}}$. If the average degree of the network satisfies $\bar{d}<2\left(2+\sqrt{2}\right)$ then $t_{\text{sw}}\rightarrow\infty$, and if $\bar{d}\geq 2\left(2+\sqrt{2}\right)$ then $t_{\text{sw}}\rightarrow 0$. To explain this, in Figure 8 we estimate the rescaled MSD, $m(t)=\Gamma(1+\alpha)\langle x^{2}(t)\rangle/t^{\alpha}$, from simulations of the CTRW on the networks $\mathcal{G}_{3}$ and $\mathcal{G}_{10}$ over an ensemble of values for $\alpha$. The waiting times are Pareto distributed, $w(t)=\alpha t_{m}^{\alpha}/t^{\alpha+1}$ for $t\geq t_{m}$, and the jump lengths are normally distributed with zero mean and variance $\sigma^{2}$. We choose the minimum value of the waiting time, $t_{m}=\left(\sigma^{2}/\Gamma\left(1-\alpha\right)\right)^{1/\alpha}$, to ensure the simulations of the CTRW agree with the multiscale analysis.101010The Laplace transform of $w(t)$ for small values of $s$ is given by $\tilde{w}(s)\approx 1-\Gamma(1-\alpha)t_{m}^{\alpha}s^{\alpha}$, thus the cofficient of $s^{\alpha}$ provides a timescale, $\tau$, upon which the Pareto distribution decays, where $\tau^{\alpha}=\Gamma(1-\alpha)t_{m}^{\alpha}$. Following the derivation of the fractional diffusion equation seen in [30] the fractional diffusion coefficient is equal to $\sigma^{2}/\left(2\tau^{\alpha}\right)$. The hopping rates in Eq. 29 are selected such that as $\epsilon\rightarrow 0$, the resulting fractional diffusion coefficient is equal to $1/2$. Thus, we select $\tau^{\alpha}=\sigma^{2}$, or equivalently $t_{m}=\left(\sigma^{2}/\Gamma\left(1-\alpha\right)\right)^{1/\alpha}$. For $\mathcal{G}_{3}$ the average degree is less than $2(2+\sqrt{2})$ and so $t_{\text{sw}}\rightarrow\infty$ as $\alpha\rightarrow 0$ (see intersection points with $m(t)=1$ in Figure 8(a)). On inspection we see that the time to converge to equilibrium behaviour increases as expected when $\alpha\rightarrow 0$. In contrast, the average degree of $\mathcal{G}_{10}$ is greater than $2(2+\sqrt{2})$, meaning $t_{\text{sw}}\rightarrow 0$ as $\alpha\rightarrow 0$ (see intersection points with $m(t)=1$ in Figure 8(b)). But if we inspect the curve in Figure 8(b) corresponding to $\alpha=0.1$, we see it has transitioned away from the short-time behaviour significantly for times as small as $t=10^{-4}$ and is only just converging to equilibrium behaviour at $t=10^{10}$, a difference of $14$ orders of magnitude. Thus, for small values of $\alpha$, it is difficult to consider a single value (such as $t_{\text{sw}}$) that represents the time taken to transition to equilibrium behaviour as the change occurs gradually over a large window of time. That said, for values of $\alpha$ closer to one the transition occurs more rapidly, and $t_{\text{sw}}$ provides a good predictive timescale.

5 Summary and discussion

In this work, we have investigated the influence of network topology on the evolution of the displacement of a particle whose position evolves according to either a diffusive or an anomalously diffusive transport process. In Section 2 we define our notion of displacement of a particle as the winding distance, which is a measure of how far a particle has travelled as it winds around the network. To investigate the evolution of the winding distance we introduced a topological structure we term the displacement tree and noted that the winding distance of a particle is equivalent to the displacement of a particle along the tree. Furthermore, in the long-time limit, a link between the displacement of a particle on the tree and the displacement of a particle on the half line with periodic bias was highlighted. The simpler transport process on the half line has a periodic bias with strength $\rho^{+}_{\infty}$, a quantity that can be calculated in terms of the average degree of the network.

In Section 3 a multiple scales approach from [26] was used to derive an advection-diffusion PDE for the displacement of a particle whose position evolves according to a Brownian motion on the half line with periodic bias. The advective and diffusive transport coefficients are calculated in terms of the periodic bias $\rho^{+}_{\infty}$. From these coefficients a topological condition was derived for whether the long-time behaviour of the winding distance switches to ballistic motion or remains diffusive. By comparing the MSD for the short-time and long-time limits, a prediction of the timescale upon which a network induces the qualitative switch from diffusive to ballistic behaviour was obtained. For an ensemble of $100$ randomly generated networks we found that the predictive timescale performs better for networks with a higher average degree.

Finally, in Section 4, we discussed an extension to a class of anomalously diffusive transport processes using the CTRW framework. Through the use of the fractional master equation we derived a time fractional advection-diffusion PDE for the evolution of the displacement of a particle on the half line with periodic bias. Our results showed that the characteristic nature of the long-time behaviour depends upon the transport process itself. We found that, depending on the value of the exponent $\alpha$, the displacement of a particle in the equilibrium limit can remain sub-diffusive, become diffusive, or become super-diffusive. Similarly to Section 3 we derived a predictive timescale upon which these transitions occur and explored its validity as a function of the parameter $\alpha$, the exponent in the heavy-tailed waiting time distribution of the CTRW.

This work introduced the winding distance, a new measure of displacement for particles on networked topologies. Using the winding distance to investigate both diffusive and anomalously diffusive transport processes has highlighted that both the topology of the network and the nature of the transport process itself play critical roles in the evolution of particle displacement. Active transport through networked environments is a critical feature of many biological systems, such as the cytoplasmic streaming seen within Drosophilia oocytes [33, 34], the transport of bronchial mucus [35], and the transport of a broad range of biomolecules along microtubular networks [36, 37]. Explicitly incorporating network topology within these biological transport systems requires the theoretical study of models for stochastic transport in networked environments [38, 39]. Therefore, future work will extend our investigation of the winding distance to transport processes with a directional bias along each of the edges in a network.

Appendix A Simulating the transport processes

Here we present the algorithm used to compute realisations of a Brownian particle on a network $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$, using the Euler-Maruyama method.

1:Assign orientations to each edge in the network, select an initial vertex for the particle, set $t=0$, the termination time $T$ and the time step $\Delta t$.

2:Sample an adjacent edge, $e$, to the initial vertex uniformly at random. Let $X(0)=0$ or $X(0)=L_{e}$ depending on the orientation of edge $e$.

3:while $t<T$ do

4:     if $t=0$ then update the position as follows:

[TABLE]

where $\xi\sim\mathcal{N}(0,1)$.

5:     else  update the position as follows:

[TABLE]

where $\xi\sim\mathcal{N}(0,1)$.

6:     end if

7:     if  $X(t+\Delta t)<0$ then sample a new edge $e^{*}$, uniformly at random from the adjacent edges to edge $e$ (not including $e$ itself) at the end with position [math].

8:         if  the particle enters the new edge $e^{*}$ at the end with position [math] then update the new position as

[TABLE]

9:         else if  the particle enters the new edge $e^{*}$ at the end with position $L_{e^{*}}$ then update the new position as

[TABLE]

10:         end if

11:         Update $e=e^{*}$.

12:     else if $X(t+\Delta t)>L_{e}$ then sample a new edge $e^{*}$, uniformly at random from the adjacent edges to edge $e$ (not including $e$ itself) at the end with position $L_{e}$.

13:         if  the particle enters the new edge $e^{*}$ at the end with position [math] then update the new position as

[TABLE]

14:         else if  the particle enters the new edge $e^{*}$ at the end with position $L_{e^{*}}$ then update the new position as

[TABLE]

15:         end if

16:         Update $e=e^{*}$.

17:     end if

18:     Update the time $t=t+\Delta t$.

19:end while

For the CTRWs seen in Section 4 we select the jump distribution to be normal with zero mean and standard deviation $\sigma$. The above algorithm is easily adapted to simulate the CTRW on a network. Instead of fixed time steps, $\Delta t$, update the time $t$ with random waiting times from the Pareto distribution, and replace $\sqrt{2D\Delta t}$ in equations Eq. 37 and Eq. 38 with $\sigma$. The CTRW simulations must agree with the multiscale analysis in Section 4, on both long and short timescales. In particular, on timescales sufficiently short that a particle has not left the initial edge. The multiscale analysis discretises the unit interval into $N+1$ lattice sites. As $N\rightarrow\infty$ the width of each site, $\epsilon=1/N$, tends to zero, and a particle will undergo a large number of jumps before traversing an entire interval. Therefore, for the CTRW to agree with the multiscale analysis on short timescales, we take $\sigma^{2}\ll 1$, such that many jumps must occur before a particle traverses an edge.

Appendix B Periodic networks

Recall the sequence of biased probabilities $\{\rho^{+}_{1},\rho^{+}_{2},\ldots\}$ in Section 2. The effective bias for a particle with displacement $kL$ is given by

[TABLE]

where $p_{k}(\nu)$ is the probability of being at vertex $\nu$ given a displacement of $kL$. We previously assumed that the network was aperiodic when viewed as a CTMC, so that this sequence was guaranteed to converge to the unique value $\rho^{+}_{\infty}$ given in Eq. 2. However, if instead the network is periodic then the probability $p_{k}(\nu)$ converges to a periodic sequence $\{p_{\infty,1}(\nu),\ldots,p_{\infty,m}(\nu)\}$ where $m$ is the period. Weighting these probabilities by $(d_{\nu}-1)/d_{\nu}$ we calculate the effective probabilities to move to the right along the displacement tree as

[TABLE]

for $j\in\{1,\ldots,m\}$. We are now able to use the multiple scales approach of Section 3 for the case of periodic networks. The periodic unit interval must be extended to have $mN+1$ lattice sites, and there will now be internal lattice sites that have asymmetric rates. The analysis proceeds identically to before, the only changes are in the size and entries of the matrices in Eq. 7.