Persistent exclusion processes: inertia, drift, mixing and correlation

TL;DR

This paper develops and analyzes lattice-based models for motile agents with persistence and crowding effects, revealing nonlinear diffusion behavior and extending to multiple species and drift, with validation against simulations.

Contribution

It introduces a mean-field framework for persistent exclusion processes, deriving nonlinear PDEs and generalizing to multiple species and drift, advancing understanding of collective motion with crowding.

Findings

Models are well described by nonlinear diffusion equations.

Persistence and exclusion cause nonlinearity in diffusion behavior.

Mean-field predictions align with stochastic simulations.

Abstract

In many biological systems, motile agents exhibit random motion with short-term directional persistence, together with crowding effects arising from spatial exclusion. We formulate and study a class of lattice-based models for multiple walkers with motion persistence and spatial exclusion in one and two dimensions, and use a mean-field approximation to investigate relevant population-level partial differential equations in the continuum limit. We show that this model of a persistent exclusion process is in general well described by a nonlinear diffusion equation. With reference to results presented in the current literature, our results reveal that the nonlinearity arises from the combination of motion persistence and volume exclusion, with linearity in terms of the canonical diffusion or heat equation being recovered in either the case of persistence without spatial exclusion, or spatial exclusion without persistence. We generalise our results to include systems of multiple species of interacting, motion-persistent walkers, as well as to incorporate a global drift in addition to persistence. These models are shown to be governed approximately by systems of nonlinear advection-diffusion equations. By comparing the prediction of the mean-field approximation to stochastic simulation results, we assess the performance of our results. Finally, we also address the problem of inferring the presence of persistence from simulation results, with a view to application to experimental cell-imaging data.

Tables

Table 1. Table 1: Master equations (a) in discrete lattice coordinates (b) in continuous coordinates in the original basis (the equation for 𝑷 𝑷 \bm{P} ) and (c) the eigenvector basis (the equation for 𝚽 𝚽 \bm{\Phi} ). In (b) and (c), terms which are 𝒪 ​ ( τ 2 ) 𝒪 superscript 𝜏 2 \mathcal{O}(\tau^{2}) or 𝒪 ​ ( Δ 3 ) 𝒪 superscript Δ 3 \mathcal{O}(\Delta^{3}) are negligible in the continuum limit compared to terms exhibited here and are not shown. Here 𝑷 𝑷 \bm{P} and 𝚽 𝚽 \bm{\Phi} and their partial derivatives are all evaluated at ( x , y , t ) 𝑥 𝑦 𝑡 (x,y,t) .

	${𝑷}_{n+1}(i,j)$	$=\left(1-P\right){𝑷}_n(i,j)+\frac{P}{4}\left\{\mathrm{diag}\left[\begin{array}{c}1+\phi \\ 1-\phi \\ 1\\ 1\end{array}\right]C_n(i+1,j)+\mathrm{diag}\left[\begin{array}{c}1-\phi \\ 1+\phi \\ 1\\ 1\end{array}\right]C_n(i-1,j)+\mathrm{diag}\left[\begin{array}{c}1\\ 1\\ 1+\phi \\ 1-\phi \end{array}\right]C_n(i,j+1)+\mathrm{diag}\left[\begin{array}{c}1\\ 1\\ 1-\phi \\ 1+\phi \end{array}\right]C_n(i,j-1)\right\}{𝑷}_n(i,j)$

Table 2. Table 2: Master equation for the generalised two-parameter model.

Equations

$${\mathbf{P}}_n(i,j)={\left[\begin{array}{cccc}{R}_n(i,j)& L_n(i,j)& U_n(i,j)& D_n(i,j)\end{array}\right]}^{\top },$$

$$C_n(i,j)=R_n(i,j)+L_n(i,j)+U_n(i,j)+D_n(i,j).$$

$${\mathbf{P}}_{n+1}(\mathbf{k})={\mathbf{T}}_n(\mathbf{k}|\mathbf{k}){\mathbf{P}}_n(\mathbf{k})+\sum _{{\mathbf{k}}^{\prime }\in \mathcal{N}(\mathbf{k})}{\mathbf{T}}_n(\mathbf{k}|{\mathbf{k}}^{\prime }){\mathbf{P}}_n({\mathbf{k}}^{\prime }),$$

$$\tau \frac{\partial \mathbf{P}}{\partial t}$$

$$+\frac{P[1-C(x,y,t)]}{4}\{\left[\begin{array}{cccc}-3+\phi & 1-\phi & 1& 1\\ 1-\phi & -3+\phi & 1& 1\\ 1& 1& -3+\phi & 1-\phi \\ 1& 1& 1-\phi & -3+\phi \end{array}\right]\mathbf{P}(x,y,t)+\Delta \left(\left[\begin{array}{cccc}-1-\phi & -1+\phi & -1& -1\\ 1-\phi & 1+\phi & 1& 1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\frac{\partial \mathbf{P}}{\partial x}+\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ -1& -1& -1-\phi & -1+\phi \\ 1& 1& 1-\phi & 1+\phi \end{array}\right]\frac{\partial \mathbf{P}}{\partial y}\right)$$

$$+\frac{{\Delta }^2}{2}\left(\left[\begin{array}{cccc}1+\phi & 1-\phi & 1& 1\\ 1-\phi & 1+\phi & 1& 1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\frac{{\partial }^2\mathbf{P}}{\partial x^2}+\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 1& 1+\phi & 1-\phi \\ 1& 1& 1-\phi & 1+\phi \end{array}\right]\frac{{\partial }^2\mathbf{P}}{\partial y^2}\right)\}$$

$$\begin{array}{rl}\tau \frac{\partial \mathbf{\Phi }}{\partial t}& =\frac{P}{4}\left\{2\phi \Delta \left(\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 1& 0\\ 1/2& 1/2& 0& 0\\ 0& 0& 0& 0\end{array}\right]\frac{\partial C}{\partial x}+\left[\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 0& -1\\ 0& 0& 0& 0\\ 1/2& -1/2& 0& 0\end{array}\right]\frac{\partial C}{\partial y}\right)+{\Delta }^2\left({\nabla }^{2}C\right)\mathbf{I}\right\}\mathbf{\Phi }(x,y,t)\\ & +\frac{P[1-C(x,y,t)]}{4}\left\{\mathrm{diag}\left[\begin{array}{c}0\\ -4\\ -4+2\phi \\ -4+2\phi \end{array}\right]\mathbf{\Phi }(x,y,t)+\Delta \left(\left[\begin{array}{cccc}0& 0& -2\phi & 0\\ 0& 0& -2\phi & 0\\ -2& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\frac{\partial \mathbf{\Phi }}{\partial x}+\left[\begin{array}{cccc}0& 0& 0& -2\phi \\ 0& 0& 0& 2\phi \\ 0& 0& 0& 0\\ -2& 0& 0& 0\end{array}\right]\frac{\partial \mathbf{\Phi }}{\partial y}\right)+\frac{{\Delta }^2}{2}\left(\left[\begin{array}{cccc}2& 0& 0& 0\\ 2& 0& 0& 0\\ 0& 0& 2\phi & 0\\ 0& 0& 0& 0\end{array}\right]\frac{{\partial }^2\mathbf{\Phi }}{\partial x^2}+\left[\begin{array}{cccc}2& 0& 0& 0\\ -2& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 2\phi \end{array}\right]\frac{{\partial }^2\mathbf{\Phi }}{\partial y^2}\right)\right\}\end{array}$$

$$\mathbf{P}(x,y,t)={\left[\begin{array}{cccc}R(x,y,t)& L(x,y,t)& U(x,y,t)& D(x,y,t)\end{array}\right]}^{\top },$$

$$C(x,y,t)=R(x,y,t)+L(x,y,t)+U(x,y,t)+D(x,y,t).$$

$$\left[\begin{array}{cccc}-3+\phi & 1-\phi & 1& 1\\ 1-\phi & -3+\phi & 1& 1\\ 1& 1& -3+\phi & 1-\phi \\ 1& 1& 1-\phi & -3+\phi \end{array}\right]$$

$${\mathbf{M}}_{\mathrm{P}\to \Phi }=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 1& -1& -1\\ 1& -1& 0& 0\\ 0& 0& 1& -1\end{array}\right],$$

$$\mathbf{\Phi }={\mathbf{M}}_{\mathrm{P}\to \Phi }\mathbf{P}=\left[\begin{array}{c}R+L+U+D\\ (R+L)-(U+D)\\ R-L\\ U-D\end{array}\right].$$

$$\mathbf{\Phi }={\left[\begin{array}{ccccccc}{\varphi }_1& & \Delta {\varphi }_2& & \Delta {\varphi }_3& & \Delta {\varphi }_4\end{array}\right]}^{\top },$$

$$\mathbf{P}=\frac{{\varphi }_1}{4}\left[\begin{array}{c}1\\ 1\\ 1\\ 1\end{array}\right]+\frac{\Delta {\varphi }_2}{4}\left[\begin{array}{c}1\\ 1\\ -1\\ -1\end{array}\right]+\frac{\Delta {\varphi }_3}{2}\left[\begin{array}{c}1\\ -1\\ 0\\ 0\end{array}\right]+\frac{\Delta {\varphi }_4}{2}\left[\begin{array}{c}0\\ 0\\ 1\\ -1\end{array}\right].$$

$$(1-C)P{\varphi }_2(x,y,t)$$

$$\frac{P}{4}\left\{\frac{\partial C}{\partial x}\phi {\varphi }_1-2(1-C)\left[(2-\phi ){\varphi }_3+\frac{\partial {\varphi }_1}{\partial x}\right]\right\}$$

$$\frac{P}{4}\left\{\frac{\partial C}{\partial y}\phi {\varphi }_1-2(1-C)\left[(2-\phi ){\varphi }_4+\frac{\partial {\varphi }_1}{\partial y}\right]\right\}$$

$$\tau \frac{\partial {\varphi }_1}{\partial t}+\mathcal{O}({\tau }^2)=\frac{P{\Delta }^2}{4}\left\{{\nabla }^{2}C{\varphi }_1+2\phi \nabla C\cdot \left[\begin{array}{c}{\varphi }_3\\ {\varphi }_4\end{array}\right]+(1-C)\left[-2\phi \left(\frac{\partial {\varphi }_4}{\partial y}+\frac{\partial {\varphi }_3}{\partial x}\right)+{\nabla }^2{\varphi }_1\right]\right\}+\mathcal{O}({\Delta }^3).$$

$$\frac{\partial u}{\partial t}=\nabla \cdot \left[\mathcal{D}(u)\nabla u\right],$$

$$\mathcal{D}(u)$$

$${\mathbf{P}}_n^{(k)}(i,j)={\left[\begin{array}{cccc}{R}_n^{(k)}(i,j)& L_n^{(k)}(i,j)& U_n^{(k)}(i,j)& D_n^{(k)}(i,j)\end{array}\right]}^{\top },$$

$$C_n^{(k)}(i,j)=R_n^{(k)}(i,j)+L_n^{(k)}(i,j)+U_n^{(k)}(i,j)+D_n^{(k)}(i,j),$$

$$C_n^{\text{tot}}(i,j)=\sum _{k=1}^M{C}_n^{(k)}(i,j).$$

$$\frac{\partial u^{(k)}}{\partial t}$$

$${\mathcal{D}}^{(k)}(\mathbf{u})$$

$${\mathbf{v}}^{(k)}(\mathbf{u})$$

$$D^{(k)}$$

$$\mathbf{\Phi }=\left[\begin{array}{c}{\varphi }_1\\ L_2+\Delta {\varphi }_2\\ L_3+\Delta {\varphi }_3\\ L_4+\Delta {\varphi }_4\end{array}\right],$$

$$\mathbf{A}=$$

$$+\frac{1-\lambda }{4}\left[\begin{array}{cccc}1+h& 1+h& 1+h& 1+h\\ 1-h& 1-h& 1-h& 1-h\\ 1+v& 1+v& 1+v& 1+v\\ 1-v& 1-v& 1-v& 1-v\end{array}\right].$$

$$\mathbf{\pi }=\left[\begin{array}{c}\frac{-2+2h(-1+\lambda )+\lambda \phi }{-8+4\lambda \phi }\\ \frac{1}{4}+\frac{h(1-\lambda )}{-4+2\lambda \phi }\\ \frac{-2+2v(-1+\lambda )+\lambda \phi }{-8+4\lambda \phi }\\ \frac{1}{4}+\frac{v(1-\lambda )}{-4+2\lambda \phi }\end{array}\right].$$

$$\begin{array}{rl}{L}_2& =({\pi }_1^{\ast }+{\pi }_2^{\ast })-({\pi }_3^{\ast }+{\pi }_4^{\ast })=0,\\ L_3& ={\pi }_1^{\ast }-{\pi }_2^{\ast }=\frac{h(1-\lambda )}{2-\lambda \phi }{\varphi }_1,\\ L_4& ={\pi }_3^{\ast }-{\pi }_4^{\ast }=\frac{v(1-\lambda )}{2-\lambda \phi }{\varphi }_1.\end{array}$$

Full text

I Introduction

Random walks are a broad class of stochastic processes relevant for modelling a wealth of phenomena Hughes (1995); Weiss (1994). An area of application which has received significant interest is the use of such models for the study of cellular motility Codling et al. (2008); Mort et al. (2016); Krummel et al. (2016); Baker et al. (2014). Much of the existing literature on random walks considers a single walker. However, in some contexts, it is useful to model a population of walkers on a lattice, for which it is necessary to account for interactions Simpson et al. (2009a); Codling et al. (2008). As a first model of a system of interacting walkers (sometimes referred to as agents), occupancy of a given lattice site by multiple agents is forbidden, with attempted moves by agents into already occupied sites being aborted. If the process evolves in continuous time with each individual agent waiting an exponentially distributed time between move attempts, any two agents will attempt simultaneous moves with probability zero and so clashing move attempts do not occur. Alternatively, we may allow the process to evolve in discrete time with a random sequential update procedure: if there are $N$ agents present, we make $N$ sequential random choices of agent with replacement, offering each chosen agent an opportunity to move. Once $N$ such offers have been made, we increment time by a constant amount. When a simple random walk stepping rule for individual agents is combined with either of these two time evolution protocols and the abortion of moves onto occupied sites, we obtain a model called the simple exclusion process (SEP). A number of rigorous results are available for this particular model Liggett (1985, 1999); Komorowski et al. (2007), the most important of these being that in an appropriate continuum limit, the probability of occupancy of a location evolves under the classical linear diffusion (or heat) equation Liggett (1999).

For most systems of interacting agents apart from the SEP, exact continuum limit partial differential equations (PDEs) have been not been established and indeed may not exist due to subtle effects of inter-agent correlation (which happen to conveniently remove themselves from some questions that can be answered for the SEP). However, useful approximate continuum limit PDEs can be obtained by mean-field arguments in which correlation effects are either neglected (see, e.g., Simpson et al. (2009a), although the literature in this area is quite extensive), or modelled in a simple approximate manner Markham et al. (2013). For the SEP, the simplest mean-field approach actually produces the correct PDE Simpson et al. (2009a, b). Further results have been developed for variations on the simple exclusion process, including incorporation of additional bias effects and accounting for interactions between multiple species of walkers Simpson et al. (2009b), as well as consideration of alternative interaction rules Penington et al. (2011); Almet et al. (2015); Nan et al. (2018); Fernando et al. (2010). In many cases the predictions of the mean-field PDE agree well with simulation averages, especially at low densities of agents, but the quality of agreement degrades when inter-agent interactions have a significant attractive component Fernando et al. (2010), or when agents reproduce at random as well as moving at random.

In many scenarios including modelling of vehicle movement Michelini and Coyle (2008), chemical reaction kinetics Vilensky et al. (1994), cell motility Baker et al. (2014) and animal movement Wu et al. (2000), it is clear that agents do not move completely at random but exhibit short-term persistence in the direction of motion, thus giving rise to the notion of a persistent random walk (PRW) Codling et al. (2008); Gorelik and Gautreau (2014); Wu et al. (2000). The case of non-interacting walkers has received substantial attention Patlak (1953); Masoliver and Lindenberg (2017); Codling et al. (2008); Renshaw and Henderson (1981); Henderson et al. (1984). In particular, it is known that the non-interacting persistent random walk is governed by a linear diffusion equation Renshaw and Henderson (1981); Henderson et al. (1984). However, current work considering multiple interacting persistent random walks remains sparse Gavagnin and Yates (2018); Teomy and Metzler (2019a, b).

In the persistent exclusion process, we consider swarms of interacting, motion-persistent agents. Mean-field approximations may be used to derive an approximate continuum description of such a system at the macroscopic level. Gavagnin and Yates Gavagnin and Yates (2018) presented such an analysis for one model formulation in which individual agents follow an explicit velocity-jump process. Agents possess one of four possible polarisations on a two-dimensional lattice, and either move following a polarisation-dependent rule for step directions, or repolarise uniformly. Movement and repolarisation are driven by independent Poisson processes, and in the continuum limit this model yields a system of four nonlinear PDEs describing the dynamics of agents in each of the four possible orientations. Even more recently, Teomy and Metzler Teomy and Metzler (2019a, b) reported results for transport in a general model of exclusion processes in which agents possess a finite memory of attempted moves and move in a memory-dependent manner. The authors showed that in the continuum limit such a system may be described by nonlinear advection-diffusion equations, verifying their results for the case of steady-state transport in one dimension.

In this paper, we present and analyse in one and two dimensions a related but distinct model for a persistent exclusion process where the motion persistence of individual agents is incorporated as a single-step memory of the last successful move taken Renshaw and Henderson (1981); Patlak (1953); Jones et al. (2015); Galanti et al. (2013); Hughes (1995); Weiss (1994), and hence the sequence of successful moves made by each agent follows a strict Markov rule. We begin by formulating our agent-based model for a number of cases, from which we then derive population-level descriptions in the continuum limit. For the basic case of a single species of motion-persistent agents, we find that such systems are approximately governed by nonlinear diffusion equations. Although our analysis yields results resembling those presented by Teomy and Metzler Teomy and Metzler (2019a), key differences are revealed regarding the descriptions of population-level behaviour. We extend our analysis to consider generalisations of the basic model involving multiple species of interacting agents as well as a superimposed global drift effect, showing that in general such systems are approximately governed by systems of nonlinear advection-diffusion equations. By comparing these continuum results to simulation data, we verify that our derived continuum approximations hold for evolving, inhomogeneous systems. Finally, we investigate methods for inferring the presence of persistence from simulation results, with a view to application to experimental cell-imaging data Mort et al. (2016); Agnew et al. (2014).

II Agent-based model

II.1 One-dimensional model

We consider $N\gg 1$ agents on $\mathbb{Z}$, the usual one-dimensional lattice of sites with integer coordinates. Agents on the lattice are able to move right (R) or left (L), and to each agent we ascribe the direction of its last successful move to be its orientation. We thus partition the agent population by orientation into two groups denoted R and L. At each timestep, we sample $N$ agents with replacement following the random sequential update procedure Simpson et al. (2009a) previously described. Each sampled agent, with probability $P\in(0,1]$, will attempt to move. If the destination site is vacant the attempted move is successful, otherwise it is aborted. Move attempts are made with the following assigned probabilities: in the direction of orientation with probability $(1+\varphi)/2$, and against the direction of orientation with probability $(1-\varphi)/2$. We illustrate this in Fig. 1(a) for a right-oriented agent. The parameter $\varphi\in[-1,1]$ controls the extent of motion persistence. We make note that the case $\varphi<0$ corresponds to a ‘negative’ persistence effect where agents have a propensity to undo their previous move.

II.2 Two-dimensional model

Consider now the square lattice $\mathbb{Z}^{2}$. We may extend the one-dimensional model to now consider four possible orientations, right (R), left (L), up (U) and down (D). An attempted move is taken in the direction of orientation with probability $(1+\varphi)/4$, against the direction of orientation with probability $(1-\varphi)/4$, and in one of the two choices orthogonal to the direction of orientation each with probability $1/4$. We illustrate this in Fig. 1(b) for an agent whose most recent move was to the right.

Note that in this two-dimensional case, the parameter $\varphi$ controls only the propensity of agents to move in the direction parallel to their orientation, and that we enforce that agents move orthogonal to their orientation with a fixed probability. We have obtained corresponding results for a more general two-parameter model from which the simpler model can be derived (see Appendix A), but for simplicity we limit our discussion to the single-parameter model.

II.3 Two-dimensional model with global drift

We generalise the two-dimensional model to incorporate both a local persistence and global drift. We introduce an additional parameter $\lambda\in[0,1]$ that specifies the relative contribution of local persistence and global drift, and parameters $h,v\in[-1,1]$ controlling the direction and magnitude of global drift in the horizontal and vertical directions respectively. An agent offered the opportunity to move chooses its move as follows. With probability $\lambda$, choose how to move as in Section II.2. Alternatively with probability $1-\lambda$, attempt a move that respects the imposed global drift: right with probability $(1+h)/4$, left with probability $(1-h)/4$, up with probability $(1+v)/4$, or down with probability $(1-v)/4$. This movement rule is illustrated in Fig. 2 for an agent whose most recent move was to the right.

It is important to note here that in our agent-based model, agents only remember successful moves, with aborted move attempts being ‘forgotten’. This is in contrast with the model presented by Teomy and Metzler Teomy and Metzler (2019a), in which the memory is updated regardless of the outcome of a move attempt. The resulting differences for bulk system behaviour are immediately appreciable. For instance, in one dimension with extreme negative persistence ($\varphi=-1$), Teomy and Metzler report that mean square displacement (MSD) exhibits $\sqrt{t}$ behaviour for starting densities above $1/2$. For our model, however, such behaviour is never possible in the equivalent setting.

III Continuum-limit analysis

III.1 Single species

For brevity, we show the derivation of the continuum-limit PDE in the two-dimensional case, as the analysis for the one-dimensional model is analogous. We model the agent-based model as a random walk with internal states Hughes (1995); Weiss (1994). As mentioned earlier, we partition agents by orientation into subpopulations henceforth labelled R, L, U and D. We introduce $\bm{P}_{n}(i,j)$ as the vector of occupancy probabilities for the lattice site $(i,j)$ at time step $n$, with

[TABLE]

where $R_{n},L_{n},U_{n}$ and $D_{n}$ are the respective occupancy probabilities for right, left, up and down-oriented agents. We also define the total occupancy $C_{n}(i,j)$ to be

[TABLE]

Under a mean-field approximation, we assume that average lattice occupancies of distinct sites are independent, and thus the evolution of $\bm{P}_{n}$ with time across the lattice is approximately governed by the equation

[TABLE]

where $\mathcal{N}(\bm{k})=\left\{\bm{k}+\bm{e_{1}},\bm{k}-\bm{e_{1}},\bm{k}+\bm{e_{2}},\bm{k}-\bm{e_{2}}\right\}$ denotes the set of sites immediately adjacent to $\bm{k}$, and $\bm{e_{1,2}}$ denote the basis vectors in the $x$ and $y$ directions respectively. $\bm{T}_{n}(\bm{k}|\bm{k^{\prime}})$ is a matrix of transition probabilities for the transition from site $\bm{k^{\prime}}$ to $\bm{k}$, which is in general not constant but dependent on $\bm{k},\bm{k^{\prime}}$ and $n$ as well as $\bm{P}_{n}$. From this, we obtain the master equation in Table 1(a).

Following the typical procedure Renshaw and Henderson (1981); Codling et al. (2008); Simpson et al. (2009a, b); Gavagnin and Yates (2018), we introduce continuous coordinates $(x,y,t)$, a lattice spacing $\Delta$, and time step $\tau$. We take $x=\Delta i,y=\Delta j,t=\tau n$. By making a formal substitution of variables, our vector of occupancy probabilities becomes

[TABLE]

and the total occupancy becomes

[TABLE]

We then substitute formally into Table 1(a) and take Taylor expansions of $\bm{P}$ and $C$ to first order in $t$ and second order in $(x,y)$. Cancelling a term in $\bm{P}(x,y,t)$ from both sides and rearranging, we arrive at the system in continuous variables shown in Table 1(b).

We now perform a judicious change of basis for the system into the eigenvectors of the matrix

[TABLE]

found to feature prominently within the continuous-variable master equation in Table 1(b). We use the change of basis matrix

[TABLE]

to form from $\bm{P}$ the transformed quantity $\bm{\Phi}$ given by

[TABLE]

Note that all components of $\bm{\Phi}$ except the first correspond to differences between subpopulation occupancy probabilities. By symmetry of the problem, we reason that as $\Delta\to 0$ in the continuum limit, these components should vanish since no individual direction is preferred above others. We thus formulate an ansatz for the form of $\bm{\Phi}$ in which the vanishing components are scaled by $\Delta$, that is,

[TABLE]

where $\phi_{i}=\phi_{i}(x,y,t)$ for $1\leq i\leq 4$. From this, we obtain

[TABLE]

Transforming the system in Table 1(b), we obtain the simplified system in Table 1(c). We first apply the continuum limit $\Delta,\tau\to 0$ for components $\phi_{2},\phi_{3},\phi_{4}$, which respectively yields

[TABLE]

Further, the equation for the first component $\phi_{1}$ (without taking any limits) transpires to be

[TABLE]

Solving (11), (12), (13) for $\phi_{2},\phi_{3},\phi_{4}$, we differentiate to obtain $\partial_{y}\phi_{4}$ and $\partial_{x}\phi_{3}$. Substituting into (14), we take the continuum limit in the usual manner holding $\Delta^{2}/\tau$ constant Codling et al. (2008); Simpson et al. (2009a, b); Renshaw and Henderson (1981); Gavagnin and Yates (2018). Noting that $C=\phi_{1}$, we obtain a nonlinear diffusion equation of the general form

[TABLE]

where we take $u=\phi_{1}$ and the nonlinear diffusivity is given by

[TABLE]

In the absence of persistence, we set $\varphi=0$ and our diffusivity reduces to $\mathcal{D}(u)=D$, that is, we recover the linear diffusion equation as a description of the simple exclusion process Simpson et al. (2009a, b).

In the case of the one-dimensional model, similar arguments also produce the nonlinear diffusion equation given in (15) where we interpret $\bm{\nabla}$ as acting in one dimension, and

[TABLE]

We note in both cases that for $\varphi>0$ ($\varphi<0$), the diffusivity $\mathcal{D}$ is linear in the density $u$ and decreases (increases) with increasing density, whereas for $\varphi=0$ in the non-persistent case, the diffusivity remains density-independent. Holding $u$ constant, the diffusivity is also an increasing function of the persistence parameter $\varphi$. Introduction of motion-persistence has thus led to the nontrivial result of a diffusivity intriguingly dependent on both the density $u$ and extent of persistence $\varphi$.

In the one-dimensional case, the diffusivity becomes negative for $u>(2\varphi)^{-1}$. Together with the observation that on a one-dimensional lattice the relative ordering of agents must remain fixed, we conclude that at large $\varphi$ and $u$, spatial correlation of site occupancies on the one-dimensional lattice become significant and so mean-field arguments employed in this analysis are invalidated. We also observe that $\mathcal{D}$ becomes unbounded as $\varphi\to 1$. This corresponds qualitatively to a degenerate case for the one-dimensional discrete process in which agents are unable to change orientation and so are forced to move in a fixed direction. Noting that non-interacting walkers subject to this movement rule would exhibit ballistic rather than diffusive motion, we suggest that this is in some way reflected in the blowing up of $\mathcal{D}$. Thus, we expect the continuum results in one dimension to be useful mostly for relatively dilute systems, or for low persistence. For the two-dimensional case the diffusivity is always non-negative.

Whilst our analysis yields a diffusivity linear in the density $u$ for both one and two dimensions, the analogous result for diffusivity in one dimension obtained by Teomy and Metzler Teomy and Metzler (2019a) varies quadratically with density. In particular, the authors reported a critical density ($1/2$ in the mean-field approximation), below which the diffusivity increases with increasing persistence, and above which the diffusivity instead decreases with increasing persistence. In contrast, our model displays no such behaviour, suggesting that the reported behaviour is characteristic of the specific model formulation used where agents keep track of the last attempted move.

The model presented by Gavagnin and Yates Gavagnin and Yates (2018) yielded systems of four coupled PDEs which required explicit solutions for each subpopulation in the continuum limit. On the other hand, analysis of memory-based systems yields only a single equation in the total population density, e.g. (15) which naturally lends itself to interpretation as a nonlinear diffusion equation. This is a result of the scaling argument (9), where we take the limit in which the occupancy probability at each position is uniformly distributed among the possible orientations. In Section IV, we find that the match between the resulting continuum description and simulation is very good. Nonetheless, we hypothesise that in principle it may be possible to derive a system of PDEs describing the individual subpopulations and thus account for cases such as initial conditions with non-uniform distribution of agent orientation. From simulations of the agent-based model however, we find that for mild persistence such initial anisotropies generally dissipate quickly to recover the isotropic case considered in our continuum-limit analysis.

Curiously, it is known that both the non-interacting persistent random walk Renshaw and Henderson (1981); Henderson et al. (1984) and the non-persistent exclusion process Simpson et al. (2009a) can be described by a canonical linear diffusion equation in the continuum limit. Based on this knowledge, our findings indicate that only upon combining both persistence and exclusion does the system yield behaviour governed by a nonlinear diffusion equation.

III.2 Multiple species

Consider now $M$ distinct species of random walkers on the two-dimensional lattice, each of which possess species-specific parameter values for motility $P^{(k)}$ and persistence $\varphi^{(k)}$, $1\leq k\leq M$. Subject to these parameters, agents of each species are governed by the movement rules described in Section II.

For the $k$th species, we introduce the vector of occupancy probabilities

[TABLE]

where $R^{(k)}_{n}(i,j)$ is the occupancy probability for right-oriented agents of species $k$, and so on. Further, we define the total occupancy probability for the $k$th species (across all four orientations) to be

[TABLE]

and the total occupancy probability across all species to be

[TABLE]

We perform the analysis for the multi-species case parallel to the previous analyses, except we now consider $\bm{P}^{(k)}_{n}(i,j),P^{(k)},$ and $\varphi{(k)}$ for each individual species $k$. In the continuum limit, we obtain a system of coupled advection-diffusion PDEs for $1\leq k\leq M$ described generally by

[TABLE]

where we have written $\bm{u}=\left(u^{(1)},\dots,u^{(M)}\right)$ and the diffusivity and velocity field are given respectively by

[TABLE]

with

[TABLE]

where for brevity we have written $C^{\mathrm{tot}}$ as defined in (20).

The diffusivity term $\mathcal{D}^{(k)}$ is similar to that obtained for a single species, except that it is now always a decreasing linear function in the total density $C^{\text{tot}}$ regardless of the value of $\varphi$. In contrast to the case of a single species, each individual species is now subject to a velocity field $\bm{v}^{(k)}$ oriented parallel to $-\bm{\nabla}C^{\text{tot}}$, i.e. along the direction of steepest decrease in total agent density. This is consistent with the exclusion effect exerted by the total population on agents of each individual species. As expected, by setting $M=1$ we recover the single-species case (15, 16), and we note that in the non-persistent case $\varphi=0$, we recover a previously known result derived by Simpson et al. Simpson et al. (2009b).

III.3 Incorporating global drift

For the two-dimensional model with global drift described in Section II.3, the analysis of Sections III.1 and III.2 require some key modifications due to the drift. A discrete master equation is formulated which, under the typical analysis, we convert into continuous variables $(x,y,t)$ and take Taylor expansions. As previously, a transformation (7) into the natural coordinates of the system is then performed. With the presence of drift, we can no longer apply the scaling arguments presented in (9) since the presence of drift pushes the system away from isotropy. Instead, we put forward an ansatz for the form of $\mathbf{\Phi}$ for the current system to be

[TABLE]

where $L_{2},L_{3},L_{4}$ are non-negative functions of $\phi_{1}$ only, and $\phi_{k}=\phi_{k}(x,y,t),1\leq k\leq 4$ are functions in $(x,y,t)$. We next consider the case of a single isolated agent, in which all attempted moves are successful, and formulate the transitions between orientations as a discrete-time Markov chain, for which the state space comprises of the four available orientations $\{R,L,U,D\}$. This is described by the transition matrix

[TABLE]

Solving for the stationary distribution $\bm{\pi}$ of the system described by $\bm{A}$ yields the following, where entries correspond respectively to the proportions of time spent by the walker in right, left, up and down orientations under an ergodic interpretation:

[TABLE]

We now suppose that the stationary distribution $\bm{\pi}$ applies also to the case of interacting agents. Noting that $\phi_{1}$ corresponds to the total occupancy probability across all orientations, we scale $\bm{\pi}$ by $\phi_{1}$ to obtain in absolute terms the limiting occupancy probabilities for each orientation, $\bm{\pi}^{*}=\phi_{1}\bm{\pi}$, and so we have that

[TABLE]

We thus arrive at an ansatz for the form of $\bm{\Phi}$,

[TABLE]

In order to take the continuum limit, we must further enforce a scaling for weak drift Simpson et al. (2009a, b) i.e. $h,v\in\mathcal{O}(\Delta)$ to ensure that drift-related terms remain well-behaved in the limit,

[TABLE]

Carrying out the remaining analysis using the ansatz (29) yields the resulting advection-diffusion PDE in the continuum-limit for the single species case,

[TABLE]

with diffusivity and velocity field given by

[TABLE]

where $D$ is as defined in Section III.1. We note that upon setting $\lambda=1$, we recover the previously presented result (16), corresponding to zero drift. For nonzero $\lambda$, we obtain also a velocity field $\bm{v}$ representing the drift velocity parallel to $(H,V)$. Further, the magnitude of each component of $\bm{v}$ is an increasing function of $\varphi$, indicating that motion-persistence reinforces the presence drift.

In the multispecies case, applying the arguments in Section III.2 yields a system of advection-diffusion PDEs described by (21) with diffusivity and velocity field respectively given by

[TABLE]

where $1\leq k\leq M$ for $M$ interacting species and $D^{(k)}$ is as defined previously in Section III.2). For brevity, we have written $C^{\mathrm{tot}}$ as defined in (20). We note that the velocity field $\bm{v}^{(k)}$ now contains a term in the direction of $-\bm{\nabla}C^{\text{tot}}$ arising from population exclusion as well as the drift term parallel to $(H,V)$. Setting $\lambda=1$ to eliminate drift, the diffusivity and velocity field revert to the previously obtained results (22, 23) as expected.

IV Results

The agent-based model was implemented for both the one and two-dimensional cases. To obtain qualitative insight on the effect of introducing motion persistence for the case of a single isolated walker, the agent-based simulation was run up to $t=2000$ on a $200\times 200$ lattice, starting at position $(100,100)$ with parameters $P=1,\varphi=0$ (non-persistent), $\varphi=1$ (persistent). Results displayed in Fig. 3 show that the non-persistent walker remains relatively close to its starting point, moving in both a locally and globally random fashion. In contrast, the persistent walker covers much more ground, moving in locally directed bursts that are interspersed with changes in direction, reflecting the local persistence in the presence of the globally unbiased motion characteristic of the persistent random walk Gavagnin and Yates (2018); Treloar et al. (2011); Codling et al. (2008).

IV.1 Single species

In this subsection, we examine the results of our continuum-limit analysis for homogeneous populations of motion-persistent walkers. All simulations were performed up to $t=1000$ and sampled at $t=0,100,500,1000$, for which corresponding PDE solutions were found. For the one-dimensional case, we conducted simulations on a lattice of size $200$, starting with an initially uniform agent distribution at density 0.8 on $[80,120]$ for both non-persistent ($\varphi=0$) and persistent ($\varphi=\pm 0.5$) walkers, and $P=1$. For all simulations in this work, initial agent orientations were randomly assigned with uniform probability from two and four possible orientations for one and two dimensions respectively, in keeping with the symmetry assumption of Section III.1. Simulations were averaged over 500 realisations, and the continuum-limit PDE (15, 17) was solved using the pdepe routine offered by Matlab with default settings, subject to the corresponding initial conditions and no-flux boundary conditions. The results shown in Fig. 4(a) confirm a close agreement between the agent-based discrete model and its continuum-limit approximation. From these results, we note that a positive persistence parameter $\varphi=0.5$ results in faster dispersal of agents from the initial configuration, compared with the non-persistent case where $\varphi=0$. This is concordant with our observations made earlier on the motion of isolated agents. Conversely, $\varphi=-0.5$ enforces a negative persistence and we observe that agents disperse more slowly as a result — this agrees with the fact that agents are now averse to moving in the same direction for multiple steps. As expected from Section III.1, the match between simulation and continuum results deteriorated for large $\varphi$ and high densities (not shown) and eventually solution of the PDE (15, 17) failed for sufficiently large $\varphi$ since a negative diffusivity resulted. Solutions for higher values of $\varphi$ can be accommodated by starting with a sufficiently low initial density.

We now turn our attention to the two-dimensional case. For ease of visualisation, simulations were conducted on a $200\times 20$ lattice with $y$-invariant initial conditions and lattice occupancies were averaged along columns Gavagnin and Yates (2018); Simpson et al. (2009b). For the PDE model, computing a column-averaged density amounts to integrating out $y$-dependence. We find that the column-averaged density satisfies the same PDE as the original density when we drop $y$-terms in the gradient operator, i.e. we take

[TABLE]

for column-averaged solutions. Agents were initially uniformly distributed on $[80,120]\times[1,20]$ at density 0.8. Parameters corresponding to non-persistent ($\varphi=0$) and persistent ($\varphi=\pm 1$) walkers were used, and $P=1$ held fixed. Simulations were averaged across 20 realisations, and solutions to the column-averaged PDE (15, 16) were found as previously. Snapshots of the lattice agent distribution taken at each time for a single simulation realisation for each value of $\varphi$ are shown in Fig. 4(b), and a comparison of column-averaged simulation results to PDE solutions are presented in Fig. 4(c). As previously noted, the column-averaged results provide an excellent match to simulation. A positive persistence parameter ($\varphi=1$) again results in faster dispersal of agents compared to the non-persistent ($\varphi=0$) case, whilst a negative persistence parameter ($\varphi=-1$) has the opposite effect. In particular, we note that whilst in the one-dimensional case the PDE solution deviates from simulated results for large values of the persistence parameter and eventually becomes badly behaved, in two dimensions neither of these issues are apparent. This suggests that the additional degree of freedom has non-negligible consequences for the behaviour of the system, and that a level of caution may be necessary when attempting to generalise one-dimensional results into higher dimensions.

Subsequently we employed a $200\times 200$ lattice and an initial agent distribution on $[80,120]\times[80,120]$ at density 0.8 to investigate the correspondence between discrete and continuum results in a fully two-dimensional case. Parameters $\varphi=1,P=1$ were used. Results at $t=1000$ were averaged across 20 simulation realisations and smoothed using a $10\times 10$ convolution kernel $\bm{J}_{10}$. Solutions to the PDE (15, 16) for $t=1000$ were found numerically using the solvepde routine from the Matlab PDE Toolbox. Fig. 5 shows contour plots of both the simulated and continuum density at 10%, 25%, 50% and 75% of the maximum simulated occupancy and again a close match confirms that the derived continuum description accurately captures population-level behaviour in a general two dimensional scenario.

For the two-dimensional model with global drift described in Section II.3, we performed simulations using a setup identical that used previously for Fig. 4(b, c), except parameters $P=1,\varphi=-1,0,1$ were used and $\lambda=0.9,h=0.5,v=0$ in order to enforce a weak rightwards drift. Column-averaged simulation results were averaged over 20 realisations, and corresponding numerical solutions to the PDE (31, 32, 33) were found at each time. Representative simulation snapshots for each value of $\varphi$ are shown in Fig. 6(a), and a comparison between simulation and continuum results is shown in Fig. 6(b). From this, the addition of a global drift effect is clearly visible. At long times, agents accumulate near the right boundary due to the imposed no-flux condition. We note that the presence of persistence in the $\varphi=1$ case has a synergistic effect with the rightwards drift. Thus, from both our agent-based and continuum results, we conclude that presence of a positive persistence effect enhances both the rate of agent dispersal as well as the extent of agent migration in the presence of an external drift.

IV.2 Multiple interacting species

For the case of two interacting species, we investigated the correspondence between agent-based simulations and solutions to the continuum PDE (21, 22, 23) for a range of scenarios. As in Section IV.1, unless otherwise mentioned, all simulations were performed up to $t=1000$ and sampled at $t=0,100,500,1000$. First, we consider a $y$-invariant problem with two persistent species dispersing across a $200\times 20$ lattice. Initially, agents were uniformly distributed at density 0.8 with species 1 on $[70,95]\times[1,20]$ and species 2 on $[105,130]\times[1,20]$. Both species were subject to identical parameters $P^{(1,2)}=1,\varphi^{(1,2)}=1$. Simulations results were averaged across $20$ realisations, and column-averaged continuum PDE solutions were found as previously in Section IV.1. In Fig. 7(a), we show a snapshot of a single realisation of the discrete simulation and then compare simulation to continuum results. From this, we observe that the continuum approximation agrees well with simulation data for each species. In particular, we note that the species distributions become skewed, reflecting the large-scale effect of spatial exclusion between agents of different species at the interface near $x=100$.

Following Simpson et al. Simpson et al. (2009b), we consider a mixing scenario with a non-persistent ‘filler’ species 1 distributed on $([1,79]\cup[121,200])\times[1,20]$ at density 0.3, and a persistent ‘invading’ species 2 distributed on $[80,120]\times[1,20]$ at density 0.6. We take $P^{(1,2)}=1$ and $\varphi^{(1)}=0,\varphi^{(2)}=1$. Simulations were run as for the previous example but results were averaged over 40 realisations. Column-averaged continuum PDE solutions were found as previously in Section IV.1, and results are shown in Fig. 7(b). From these results, we again confirm that the continuum model provides an accurate description of the agent-based process. As observed by Simpson et al. Simpson et al. (2009b) for non-persistent agents, we also note that continuum density profiles obtained for the mixing problem display nonmonotone behaviour as a result of interactions, as shown in Fig. 7(b, inset).

Extending our investigation more generally to a two-dimensional problem, we consider again the $200\times 200$ lattice as previously, now with a fast-moving, persistent species 1 ($P^{(1)}=1,\varphi^{(1)}=1$) interacting with a slow-moving, non-persistent species 2 ($P^{(2)}=0.3,\varphi^{(2)}=0$), initially uniformly distributed at density 0.8 on $[75,95]\times[60,140]$ and $[105,125]\times[80,120]$ respectively. Simulations were performed up to time $t=300$. Results were averaged over 20 realisations and smoothed by convolution as previously. PDE solutions for $t=300$ were found as in Section IV.1. Fig. 8(a) shows snapshots of of a single simulation at $t=0,100,200,300$, and Fig. 8(b) shows contours from the PDE and simulated results at $t=300$ at levels of 10%, 25%, 50% and 75% of the maximum simulated occupancy of species 1. We find that the continuum model largely captures the simulated behaviour, again confirming that the continuum approximation holds well for two-dimensional systems with interacting species.

Finally we consider combined drift and multi-species interactions with persistence. We found that in general the continuum approximation incorporating all three effects is more fragile and often breaks down in situations where species distributions have differing drift velocities and thus ‘collide’ with each other. We deal with a $400\times 20$ lattice with a persistent, drifting species 1 ($\varphi^{(1)}=1$) and a non-persistent drifting species 2 ($\varphi^{(2)}=0$) initially distributed on $[110,159]\times[1,20]$ and $[160,200]\times[1,20]$ respectively with density 0.5. Parameters used were $P^{(1,2)}=1,\lambda^{(1,2)}=0.75,h^{(1,2)}=1,v^{(1,2)}=0$. Representative simulation snapshots at each time are shown in Fig. 9(a). Simulations were then averaged over 40 realisations and Fig. 9(b) shows a comparison of the simulation and continuum results. We note that the species 1 distribution is more sharp and pointed compared to species 2, a result of the interplay of persistence, drift, and spatial exclusion. This can also be seen from the ‘bunching up’ of black agents in the snapshot data. In particular, we note that at $t=100$ a small density front is formed (see inset, indicated by red arrow), an effect qualitatively captured by the continuum description. This effect occurs because motion persistence amplifies drift as observed previously, and disappears when persistence for species 1 is switched off.

Whilst the continuum description has held mostly accurate for the problems considered, we observe that the continuum PDE results tend to deviate from the true simulated results near the interface between species. This is especially evident in Fig. 8 and Fig. 9, suggesting that non-negligible occupancy correlations at such interfaces may render mean-field assumptions inaccurate. To correct for such effects in these regions, further analytical or computational work may be necessary Simpson and Baker (2011); Markham et al. (2013); Teomy and Metzler (2019a).

V Inferring persistence from data

In the following discussion, we investigate the use of pair correlation functions (PCFs) for inferring the presence or absence of persistence in agent movement from simulation results. PCFs are a common tool from statistical mechanics used to characterise the average interactions in a spatial particle distribution, and have recently been applied to the context of examining cellular motility Johnston and Crampin (2019); Agnew et al. (2014); Gavagnin et al. (2018). We hypothesise that PCFs could be applied to examine motion persistence in experimentally derived data, appropriately discretised to a lattice Mort et al. (2016); Agnew et al. (2014); Deroulers et al. (2009). We first consider the case of a swarm of walkers uniformly distributed on a two-dimensional lattice, and then extend discussion to the non-uniform case of a swarm of walkers diffusing from an initial condition.

Gavagnin, Owen and Yates Gavagnin et al. (2018) recently provided relevant formulae for computing a on-lattice pair correlation function based on the $\ell_{1}$ (Manhattan) norm. Following this, we define the pair correlation function for a given distance metric $d$ as

[TABLE]

where $n_{d}(m)$ is the number of unique particle pairs observed to be separated by distance $m$, and $\bar{n}_{d}(m)$ is the number of such pairs expected to be observed if the particle distribution was truly random, referred to as complete spatial randomness (CSR) Agnew et al. (2014); Gavagnin et al. (2018). We consider a rectangular lattice of dimensions $L_{x}\times L_{y}$ and enforce periodic boundaries. The Manhattan distance between $(x,y)$ and $(x^{\prime},y^{\prime})$ is thus

[TABLE]

As in Gavagnin et al. (2018), we consider

[TABLE]

and employ the normalisation factor for periodic boundaries,

[TABLE]

where $N$ is the total number of agents present on the lattice. In the following discussion, we take $f(m)=f_{d}(m)$ to be the PCF calculated using the Manhattan norm.

V.1 Pair correlation function for a uniform swarm

Simulations were performed on a $100\times 100$ lattice, initially populated at CSR and density 0.5 by agents with parameters $P=1$ and $\varphi=-1,0,1$ respectively. For each simulation the PCF was calculated at $t=20$, at which time the system was judged to have relaxed to its steady state. Results were averaged over $1000$ realisations for each value of $\varphi$. In Fig. 10(a), we overlay the averaged PCF with the PCF calculated for a single simulation instance, which we also display as a snapshot at $t=20$. The PCF profiles reveal a striking contrast between persistent and non-persistent agents. In the non-persistent case ($\varphi=0$), we observe that the PCF remains very close to unity, indicating negligible deviation from CSR as expected. However, in the persistent case ($\varphi=1$), we observe that the PCF increases significantly above unity at relatively small separation $m<5$. We also note the presence of a minor dip in the PCF at $m\approx 6$, as highlighted in the inset. With reference to the snapshot shown, this is consistent with visible spontaneous formation of small, localised agent aggregates, as has been previously reported to result from persistence Soto and Golestanian (2014). Indeed, for higher levels of persistence achieved using the more general two-parameter model, we observe the formation of much larger aggregates (see Appendix B). Considering the case $\varphi=-1$, we note the converse observation that the PCF decreases significantly below unity for small separation $m<3$, and increases to reach a slight peak again near $m\approx 4$ (see inset).

In experimentally relevant contexts, it is useful to consider whether this method is robust given limited observation data. The PCFs computed for a single simulation overlaid in Fig. 10(a) display behaviour at small distances consistent with our observations for the averaged PCF. From this, it is clear that it is possible to distinguish persistent agent behaviour using the PCF, even from a single lattice snapshot. This observation provides strong support for the experimental utility of the PCF as a means of detecting persistence.

The effect of the persistence parameter $\varphi$ was further investigated by computing the PCF (averaged over 1000 realisations) for simulations with varying values of $\varphi=0.2,0.4,0.6,0.8,1.0$, performed up to $t=20$ with $P=1$ and agents initially distributed at CSR with density 0.5. The results shown in Fig. 11 (top) clearly show that increasing values of $\varphi$ result in roughly proportional increments in the PCF at $m=1$, corresponding to the probability of agents having an immediate ‘neighbour’. Similarly, to investigate the effect of the motility parameter $P$ on the PCF, we performed simulations for $P=0.2,0.4,0.6,0.8,1.0$, with $\varphi=1$ up to $t=50$ and the same initial configuration as previously. Resulting PCFs were averaged over 1000 realisations and shown also in Fig. 11 (bottom). These results show that the PCF is independent of $P$, and so the distribution pattern of agents is solely influenced by persistence parameter $\varphi$. In light of this, further work may investigate the PCF as a macroscopic statistic from which information about the microscopic persistence may be extracted.

V.2 Pair correlation function for a diffusing swarm

We now consider a non-uniform swarm of agents diffusing away from an initial condition. A $200\times 200$ lattice was used on which agents were initially distributed on $[80,120]\times[80,120]$ with unit density. Simulations were performed up to $t=500$ with parameters $P=1$ and $\varphi=-1,0,1$ respectively. For each value of $\varphi$, the PCF was averaged over 1000 realisations. For clarity, Fig. 10(b) displays the differences $\Delta f=f(m+1)-f(m)$ against $m$ for the averaged PCF, overlaid with the same result calculated from a single simulation. The original averaged PCF $f$ is shown in the inset, and a snapshot of a single simulation at $t=500$ is also displayed. As in the case of a uniform swarm, a sharp increase of the PCF for $m<5$ is characteristic of positive persistence ($\varphi=1$) and similarly a sharp decrease for negative persistence ($\varphi=-1$). In contrast to the case of a uniform swarm, varying values of $\varphi$ result in differing agent distribution profiles for diffusing swarms and so the overall PCF varies with both $\varphi$ and time. Despite this, the behaviour of the PCF at small $m$ remains characteristic of motion persistence (see Appendix C). These findings suggest that the PCF may be applicable for extracting information about persistence in both uniform and non-uniform swarms of motile agents.

VI Conclusion

We have presented a memory-based formulation for a lattice-based persistent exclusion process in one and two dimensions. In the basic case of a single species of motion-persistent agents, we have shown that such systems are approximately governed in the continuum limit by nonlinear diffusion equations similar to – but distinct from – a recent related result Teomy and Metzler (2019a). We note that key macroscopic differences arise from subtle differences in rules for agent memory in the respective agent-based model formulations. In two dimensions, our analysis has been extended to allow for the superposition of a global drift effect as well as interactions between multiple species of persistent walkers, resulting in systems of nonlinear advection-diffusion equations. To our knowledge, equivalent continuum-limit descriptions of the persistent exclusion process incorporating such effects have not been previously presented. We have also explored the utility of the pair correlation function (PCF) as a possible means of inferring the presence of motion persistence from simulated data. Our conclusions indicate that the PCF is a robust statistic from which persistent agent behaviour can be inferred in various settings.

Possible directions for extension of this work include application of the derived continuum models to problems of biological or other physical relevance, as well as investigation of the usefulness of the PCF for detecting motion persistence in experimental imaging data. Additionally, consideration of quantitative or empirical methods for corrected mean-field approximations Simpson and Baker (2011); Markham et al. (2013); Teomy and Metzler (2019a) could be useful for developing improved continuum models that better account for spatial correlations in both single-species and multispecies distributions of motion-persistent agents. Although we have considered only simple exclusion as our model for interactions between different agents, our approach should be able to be extended to a broad class of interactions between agents, in which the occupancies of a finite set of neighbours influence the choice of step direction in competition with the motion persistence Penington et al. (2011); Johnston et al. (2017).

Acknowledgements.

This research was supported by an Australian Research Council Discovery Project (DP140100339) and by Australian Mathematical Sciences Institute Vacation Research Scholarships. AC and BDH thank Kerry Landman for helpful discussions.

Appendix A Generalised two-parameter model

The model as presented in the main paper uses a single parameter $\varphi\in[-1,1]$ to capture the persistence in agent motion. Whilst the use of only a single parameter simplifies analysis, this may be restrictive since the probability of an agent moving orthogonal to its orientation is forced to be $0.5$. We have considered a two-parameter model which allows more control over the probabilistic rule followed by individual agents.

We now introduce two parameters $\alpha,\beta\in[0,1]$. Under this model, agents attempt moves in the direction of orientation with probability $\beta$, against the direction of orientation with probability $(1-\alpha)(1-\beta)$, and in either of the choices orthogonal to the direction of orientation with probability $\alpha(1-\beta)/2$. Thus, setting $\alpha=2/3,\beta=1/4$ we may recover the case of simple exclusion process. From this we may recover the single-parameter model using the relations

[TABLE]

A similar analysis to that presented for the one-parameter model can be performed for the two-parameter model. We show here results obtained for the model incorporating persistence and global drift, as the governing equations for the case of persistence without drift can be recovered by setting $\lambda=1$. Employing mean-field arguments, we obtain the approximate discrete master equation given in Table 2. This can be taken to its continuum limit in a similar fashion as demonstrated in the main paper to yield the PDE

[TABLE]

with diffusivity and velocity field given by

[TABLE]

where $K$ is given by

[TABLE]

and the constant $D$ is defined by

[TABLE]

For the case of multiple interacting species, we obtain the system of PDEs

[TABLE]

where the diffusivity and velocity fields are

[TABLE]

and

[TABLE]

For this general model, comparison of discrete simulation results to continuum approximations yielded similar conclusions to those presented in the main paper.

Appendix B Spontaneous aggregation for strong persistence

Using the two-parameter model, agents with a strong persistence were modelled by setting $\beta=0.9,\alpha=2/3$. Agents were initially distributed at complete spatial randomness at density 0.5, and a single simulation was performed up to $t=1000$ on a $100\times 100$ square lattice with periodic boundary conditions. Fig. 12 shows lattice snapshots at several times, revealing clear formation of slowly-evolving aggregates.

Appendix C PCF behaviour at small $m$ remains characteristic of persistence

In order to further examine the behaviour of the PCF for a diffusing swarm of persistent agents, simulations were performed as in Section V.2. in the main paper for $\varphi=0$ (non-persistent) and $\varphi=1$ (persistent) agents. Simulations were advanced up to time $t=20000$, by which time the agents were judged to have fully and uniformly dispersed. PCFs were calculated at $t=1000,5000,10000,20000$ and averaged over 100 realisations. As shown in Fig. 13, the characteristic behaviour of the PCF at small values of $m$ as discussed in the main paper is clearly evident across the full timescale of agent dispersal, confirming that it is not obscured by the bulk dispersal of agent population.