Modelling persistence of motion in a crowded environment: the diffusive limit of excluding velocity-jump processes

Understanding the properties of cell movement is of fundamental interest in many biological contexts such as embryogenesis (Mort et al., 2016), epidermal wound healing (Maini et al., 2004) and tumour growth (Sherratt and Chaplain, 2001). Mathematical models are now considered essential tools in cell biology for testing theoretical hypotheses, interpreting experimental data and extracting biological parameters (Simpson et al., 2007, 2009; Johnston et al., 2014; Ross et al., 2017). There are typically two approaches to modelling cell motion, either micro-scale discrete (Anderson and Chaplain, 1998; Simpson et al., 2007; Peirce et al., 2004; Deutsch and Dormann, 2007; Alber et al., 2003; Segovia-Juarez et al., 2004; Wang et al., 2015) or macro-scale continuum (Simpson et al., 2006; Murray, 2007; Anderson and Chaplain, 1998; Wise et al., 2008). The discrete approach, using agent-based models (ABMs), accounts for properties at the cell-scale, while the continuum approach, often presented as a system of partial differential equations (PDEs) or stochastic partial differential equations (SPDEs), gives a global description of the migration at the population-level. Continuum models have the advantage that they are generally more amenable to mathematical analysis and can lead to significant insights for situations in which the system comprises a large number of agents, at which point simulating the ABM becomes computationally expensive. Nevertheless, finding the appropriate continuum model to describe the collective behaviour of a system of moving agents can be a difficult task and continuum models are often specified on a phenomenological basis, which may reduce their predictive power. It is essential, therefore, to establish a connection between micro-scale properties, which can be inferred directly from experimental data, and macro-scale dynamics (Noble, 2002; Baker et al., 2010; Yates et al., 2012).

Many analyses of cell migration are based on the hypothesis that the movement of a single cell can be described as a simple random walk on a lattice Simpson et al. (2009); Mort et al. (2016); Baker et al. (2010); Deutsch and Dormann (2007). In many models, the behaviour of a single cell is assumed to be independent of the other cells’ positions and multiple cells can occupy the same lattice site simultaneously Codling et al. (2008); Berg and Brown (1972). In many applications, however, crowding effects play an important role that can not be neglected Simpson et al. (2007); Maini et al. (2004). Crowding is incorporated into such models via volume exclusion: each lattice site is allowed to be occupied by at most one cell Treloar et al. (2011, 2012); Slowman et al. (2016); Othmer and Stevens (1997); Stevens (2000). A macroscopic continuum description of this type of model can be obtained by considering an average mass conservation law for each lattice site and taking an appropriate limit as the spatial and temporal discretisation steps go to zero simultaneously Deutsch and Dormann (2007).

One of the key aspects of a simple random walk is that the direction of motion undergoes a series of uncorrelated jumps in space. In reality, experimental observations indicate that many types of cell tend to preserve their direction of motion for a certain time before reorienting Berg and Brown (1972); Hall (1977); Gail and Boone (1970); Wright et al. (2008), even though the movement is globally unbiased. This tendency is normally called persistence of motion or of direction Patlak (1953). There is vast literature about modelling persistence at multiple-scales for non-interacting agents Othmer et al. (1988); Othmer and Hillen (2000, 2002); Codling et al. (2008); Campos et al. (2010), but it is only in recent years that there has been an increasing interest in studying the role of persistence of motion for systems of self-interacting agents Treloar et al. (2011, 2012); Slowman et al. (2016); Thompson et al. (2011); Sepúlveda and Soto (2016); Soto and Golestanian (2014); Peruani et al. (2006); Redner et al. (2013) and persistence induced by crowding Grima et al. (2010); Smith et al. (2017).

Typically, to incorporate persistence in the ABM, cell movement is represented as a correlated random walk (CRW) which is also known as a velocity jump process Codling et al. (2008); Othmer et al. (1988). In this model, the cell has an assigned direction of motion (left or right in one dimension) and it moves in this direction with constant velocity, $v$, until the assigned direction is changed, which occurs according to a Poisson process with a given rate, $\lambda$. Notice that the temporary preferential direction induces bias in the motion for short time-scales, which represents persistence, but the resulting motion remains globally unbiased for longer time scales.

In the case of non-interacting agents, the macro-scale behaviour of the velocity-jump process in one dimension is well known to evolve according to the hyperbolic telegraph equation for the cell density $C(x,t)$ Kac (1974); Goldstein (1951); Othmer et al. (1988):

Notice that Eq. (1) was originally developed to describe the propagation of waves which travel and reflect trough a telegraph transmission line (Metzger, 2012). The same type of equation can be derived from a system of non-interacting agents performing a velocity-jump process in one dimension. In particular, Othmer et al. (1988) derive the telegraph equation for cells undergoing velocity-jump processes without interactions. Othmer and Hillen (2000) demonstrated that it is possible to obtain a parabolic limit as $v$ and $\lambda$ tend to infinity simultaneously, such that $v^{2}/\lambda$ remains constant. In this limit the canonical diffusion equation is recovered Othmer and Hillen (2000, 2002). This is not a surprise, since the short-term correlation effects become less evident at large time scales and so the limit process is effectively equivalent to a simple random walk.

When direct agent-agent interactions are introduced, however, the derivation of an exact closed form PDE for the total mean agent density is not possible (Codling et al., 2008). Recently, Treloar et al. have derived a system of macroscopic advective equations from a velocity-jump process with three different forms of direct interaction Treloar et al. (2011, 2012). Although their continuum model is successful in replicating the population-level behaviour of the ABM for a limited range of model parameters, the first order approximation considered in Treloar et al. (2011) enforces restrictions in the initial condition (which must be sufficiently smooth) and on the choice of the parameters.

The aim of our work is to ease these restrictions in order to provide a better connection between discrete and continuum models of volume excluding persistent agents. We consider a generalisation of the ABM of Treloar et al. (2011) in which we modulate the influence of persistence through an additional parameter $\varphi$. This allows us to take a diffusive limit if the new parameter $\varphi$ scales with the lattice size. The resulting PDE description includes a non-linear diffusive term, which encapsulates the long-term diffusive behaviour of cells, and an advective part, which is consistent with the findings of Treloar et al. (2011). Our new diffusive model represents an extension of the previous advective model and can be applied to study a wider range of scenarios. In particular, we can consider situations with steep gradient in cell density, which have not previously been investigated. Moreover, a diffusive limit is appropriate for the study of the long term behaviour of the system, especially if we are interested in statistical tools which are related to the diffusion coefficient, such as the mean squared displacement and the mean dispersal distance Codling et al. (2008).

In this paper we study the two-dimensional version of a model which incorporates persistence and volume exclusion. We explain the derivation of the diffusive continuum description and finally we test the agreement between our discrete and continuum models using some illustrative examples. Our new diffusive PDEs correctly represent the population-level behaviour of our ABMs, in particular in scenarios that could not be investigated with the previous advective models. Our investigation highlights some peculiar aspects of the excluding velocity-jump processes, which we discuss in the light of the new macroscopic description. In particular, a spontaneous form of aggregation, similar to that observed by Thompson et al. (2011) and Sepúlveda and Soto (2016), appears in both our agent- and population-level models. We believe this is the first reported example in which such a phenomenon appears at both micro- and macro-scales.

The paper is organised as follows. In Section 2 we define the ABM and introduce three forms of cell-cell interaction. In Section 3, we derive the continuum diffusive description of the ABM from the occupancy master equations. Our numerical results on the comparison between the ABMs and the corresponding PDEs are shown in Section 4, together with our observations on some of the interesting model behaviours. We conclude with a short discussion of our results and possible avenues for future research in Section 5.

In this section we describe the basic ABM. The models presented in the following sections are all adaptations of this basic model. Cells are represented by agents on a square lattice of size $L_{x}\times L_{y}$ sites and lattice step $\Delta$ with periodic boundary conditions in the $y$-direction and zero-flux boundary conditions in the $x$-direction²²2We implement periodic boundary conditions in the vertical direction to avoid edge-effects. In the horizontal direction we employ zero-flux boundary conditions, although in reality agents rarely, if ever, reach these boundaries. So other boundary conditions may be employed with little consequence.. Each site of the lattice can be occupied by at most one cell, in which case we will say the site is occupied, otherwise the site is said to be empty.

We assign to each agent a polarisation in one of the four directions of the lattice. We denote such polarisation with the corresponding initial capital letter: Right (R), Left (L), Up (U) and Down (D). Let $v\in\mathbb{N}^{+}$ be a positive integer which denotes the number of lattice sites that an agent can move during a single movement event. We can interpret this as a non-dimensional measure of the agent’s velocity. Agents can move or reorient their polarisation in continuous time. Both of these events occur at random as independent Poisson processes with rates $P_{m}$ and $P_{r}$, respectively. The role of the polarisation is to induce a temporary bias in the stochastic motion so that the polarised agent is more likely to move in the corresponding direction. Let $\varphi\in[0,1]$ be a parameter which characterises the intensity of the bias. Consider, for example, an R-polarised agent in two dimensions, located at site $(i,j)$. If the agent is chosen to move, one of the four sites $(i\pm v,j),\,(i,j\pm v)$ is selected, at random, as a target site. The right-hand site $(i+v,j)$, corresponding to the R-polarisation of that cell, is chosen with probability given by $\frac{1+\varphi}{4}$. In the opposite direction to the polarisation, the left-hand site $(i-v,j)$ is chosen with probability $\frac{1-\varphi}{4}$ and each of the sites in the vertical direction $(i,j\pm v)$ (orthogonal to the polarisation direction) are chosen with probability $\frac{1}{4}$ (see the schematic in Fig. 1)³³3Alternative transition probabilities could be considered. For example one could choose to reduce the movement probabilities in the direction orthogonal to the polarisation (Up and Down in the example of an R-polarised agent) as well as in the direction opposite to the current polarisation. A similar derivation can be applied, which will lead to a similar but slightly altered macroscopic model.. The transition probabilities for the other polarisations are obtained analogously. Finally, if a reorientation event occurs, with rate $P_{r}$, the agent changes its polarisation uniformly at random to one of the four possible polarisations (including the possibility of maintaining its current polarisation), see Fig. 1.

Notice that if $\varphi=0$, then the target site is chosen uniformly and the movement corresponds to a classic uncorrelated random walk (with non-local jumps for $v>1)$. On the contrary, if we let $\varphi=1$ we achieve the strongest bias where agents cannot move in the opposite direction to their polarisation⁴⁴4Note that if the model is specified in one dimension and $\varphi=1$, the model corresponds to the velocity-jump process described in Treloar et al. (2011). However, in the two-dimensional case the choice of maximum bias leads to a different model. In particular, the target site is still chosen at random (although not uniformly) between three of the nearest neighbours while in Treloar et al. (2011), upon an agent being selected to move, its target site is chosen deterministically..

When the rates of reorientation and movement are chosen such that $P_{r}\ll P_{m}$ together with a large value of the parameter $\varphi$, agents persist in their direction of motion. Fig. 2 shows two trajectories of a single agent for parameters $P_{m}=1$, $P_{r}=0.05$, $v=1$ and for panel 2 $\varphi=0$ (no persistence), for panel 2 $\varphi=0.8$ (strong persistence). Persistence of motion is clearly evident in the shape of the track in panel 2. The long term behaviour is unbiased, since none of the four directions is preferred in the long term. As there is only one agent in the domain, exclusion (specified in the next paragraph) does not play any role in the dynamics.

Once the target site is selected, the agent moves according to the exclusion property specified for the process. For consistency with Treloar et al. (2011), we consider four different exclusion properties, one without agent interaction and three with a variety of interactions. Fig. 3 shows two typical scenarios (for $v=3$) in which an agent (red) at position $i$ attempts to move to the target site at $i+3$. In Scenario A the target site is occupied by another agent (blue), while in Scenario B the target site is empty and the site $i+2$ is occupied. We use these two examples to explain the four exclusion properties as follows.

In this section we derive a family of PDEs which describe the behaviour of the ABMs introduced in Section 2 at the population-level. The general technique consists of writing down the continuous-time occupancy master equation for the average occupancy of a general site of the ABM, Taylor expanding and finally taking a limit as the lattice step, $\Delta$, and the time step, $\tau$, go to zero jointly while holding $\Delta^{2}/\tau$ constant (Deutsch and Dormann, 2007, Chapter 4). The four types of agent interactions considered lead to different PDEs, therefore, in what follows, we will consider those cases in different subsections.

We denote by $C^{t}_{i,j}(m)$ the occupancy of the site $(i,j)$ at time $t$ for the $m$-th simulation of the ABM, i.e. $C^{t}_{i,j}(m)=1$ if the site $(i,j)$ at time $t$ in the $m$-th simulation is occupied and $C^{t}_{i,j}(m)=0$ if it is empty. We denote with $\hat{C}^{t}_{i,j}$ the occupancy of the site $(i,j)$ at time $t$ averaged over all $M$ realisations, i.e.:

In addition, we define occupancy variables for the four subpopulations of agents classified according to their polarisation (Right, Left, Up, Down) and we denote these occupancies with the capital corresponding to their first letter. For example, $R^{t}_{i,j}(m)$ represents the occupancy of right-polarised agents at the site $(i,j)$ at time $t$ for the $m$-th simulation and $\hat{R}^{t}_{i,j}$ the corresponding occupancy averaged over the total number of realisations. The total occupancy can be obtained by adding together the occupancies of the four subpopulations as $\hat{C}^{t}_{i,j}=\hat{R}^{t}_{i,j}+\hat{L}^{t}_{i,j}+\hat{U}^{t}_{i,j}+\hat{D}^{t}_{i,j}$. From now on, we omit the hats in the averaged occupancies for simplicity. All the occupancies are defined for $t\in\mathbb{R}^{+}$ and $(i,j)\in\{1,\dots,L_{x}\}\times\{1,\dots,L_{y}\}$.

In this section we compare the discrete simulations of the ABM with the continuous approximation. Then we investigate how the model behaves under particular choices of the parameters. We reveal three previously unobserved aspects of the model that appear when a high level of persistence is enforced: spike formation, anisotropy and aggregation. Although such phenomena are interesting from a mathematical perspective, they represent potential obstacles for the application of such models to experimental data. We will discuss the implications of such issues and future challenges in Section 5.

All the ABMs and the corresponding PDEs are simulated on a two-dimensional domain. For the purpose of visualisation, in most examples, we show the results for column-averaged cell density profiles (i.e. averaged over the $y$ coordinates). For these examples, we define the total column-averaged density as:

and $\bar{R}$, $\bar{L}$, $\bar{U}$, $\bar{D}$ correspondingly. In these simulations we choose translationally invariant initial conditions in the vertical direction. In other words, $C(x,y,0)=C(x,0)$ for every $(x,y)\in[1,\Delta L_{x}]\times[1,\Delta L_{y}]$. Similarly, the polarised species, $R$, $L$, $U$ and $D$, are also initialised according to a translationally invariant condition. The periodic boundary conditions on the horizontal boundaries imply that translational invariance in the vertical direction is conserved as time evolves, namely:

for every $t\in\mathbb{R}^{+}$ and for every $(x,y)\in[1,\Delta L_{x}]\times[1,\Delta L_{y}]$ and similarly for the four subpopulations $R$, $L$, $U$ and $D$.

With this in mind, we can now derive the one-dimensional PDEs for the column-averaged densities from the corresponding two-dimensional equations. Formally, this is equivalent to dropping the dependence on $y$ in all the density functions. As an example, we write down the averaged PDEs for the model without interaction (type 0 interactions). We omit the expressions for the other three cases which can be obtained in a similar way. By column averaging the system (LABEL:eq:PDE_0), we obtain the following system of equations

where $\mu$ and $\nu$ are defined as in equations (6). The boundary conditions, for every $t\in\mathbb{R}^{+}$, are given by

For Figs. 4, 5, 6 and 7 we use the same computational set-up which can be described as follows. The domain is a $400\times 400$ lattice with $\Delta=1$. We impose periodic boundary conditions on the horizontal boundaries and zero-flux boundary conditions on the vertical boundaries. The ABM is simulated using the Gillespie algorithm Gillespie (1977) for $M=10$ identically prepared repeats. Few repeats are sufficient to compare the mean ABM behaviour to the solutions of the PDEs for the mean occupancy because column-averaging over 400 rows significantly reduces the noise in the ABM solutions. The numerical solutions of the PDEs are obtained through an implicit Euler method with spatial step $\delta x=0.1$ and time step $\delta t=0.1$, and using Picard iteration with tolerance $\epsilon=10^{-3}$ to solve the non-linear equations.

Fig. 4 shows the comparison between the total column-averaged density of the ABM and the PDE, for the four types of interactions at different times. The system is initialised such that the all sites with $x$ coordinate between $161$ and $240$ are populated uniformly at at random, with density $d=0.5$. The polarisation of the initial group of cells is chosen uniformly at random, so there is no bias towards any of the four polarisations. The PDE for the total column-averaged density is correspondingly initialised as constant $d=0.5$ in the interval $[160,240]$, with the density of the individual subpopulations also being constant at $d=0.125$ in this region. The parameters for the model are $P_{m}=1$, $P_{r}=0.2$, $\varphi=0.8$, $v=3$.

The agreement between the discrete and the continuous descriptions is generally very good for all four types of cell interactions, although we note that discrepancies are most noticeable for the type 3 interactions for which our assumption of independence of site occupancy is least valid. We also tested our continuous approximation for a different parameters values and found that the good agreement with the discrete model holds for a wide parameter range (results not shown). The good agreement is lost, however, as the value of persistence increases, i.e. $P_{r}\ll P_{m}$ and $\varphi\approx 1$ (see Fig. 9 (c) for an example). This disagreement is, in part, due to the significant spatial correlations induced by persistence in these models. Agents are highly likely to move to be adjacent to each other rather than aborting their movements (compare type 3 interactions to type 2 interactions, respectively). This tendency is ignored by our continuous approximations. We discuss this issue and potential improvements further in Section 5.

The hyperbolic nature of previous continuum models of persistence of motion has meant that initial conditions with steep gradients have been difficult to investigate. Due to the diffusive nature of our continuous model, we are now able to examine initial conditions which have steep density gradients and still maintain a good agreement with the discrete model. In particular, this allows us to consider the initial condition described above with only a central region uniformly populated.

To compare the behaviour of the different interaction mechanisms, we increase the initial total density to $d=0.9$ and we display numerical solutions of the four PDEs (LABEL:eq:PDE_0), (LABEL:eq:PDE_1), (LABEL:eq:PDE_2) and (LABEL:eq:PDE_3) at time $T=250$ (see Fig. 5). Non-interacting agents (type 0), lead to a faster agent spreading than any of the interacting types 1-3. Type 1 interactions cause slightly slower spread of agents: although focal agents can jump over their neighbours, a small number of type 1 movement events are aborted if there is a cell in the target site. Type 2 interactions lead to the slowest spreading due to the high proportion of aborted movements events in which the focal agent stays stationary. Agents interacting through the type 3 mechanism spread slightly faster than agents undergoing type 2 interactions: although a significant proportion of movement events are aborted when agents are immediately adjacent to the focal agent, some small movement events are permitted towards a near neighbour which would otherwise be aborted under the implementation of type 2 interactions.

In Fig. 6 we compare the ABM and the continuum approximation for the right- and left-moving subpopulations in order to evidence that the agreement between the models does not only hold at the population-level. The density profiles of the up- and down-moving subpopulations are indistinguishable and we have omitted them for simplicity. However the good of agreement between the discrete and continuum models also holds for these cases (results not shown). We find good agreement, even for large values of the reorienting parameter, $P_{r}$, which has previously been thought not to be the case (Treloar et al., 2011). Note that the appearance of loss of total mass is due to the fact that we are only visualising two of the four subpopulations in our two-dimensional model.

In the following subsections we outline some of the inherent features of these models which have largely been overlooked, but which must be considered if the model is to be used in real applications.

At the agent-level, we have modelled persistence of motion for interacting agents through excluding velocity-jump processes. Traditionally, such processes are associated with systems of advective PDEs which describes the model at a population scale. Our work continues the investigation of this type of model. First of all, we generalised the traditional velocity-jump model which allowed us to derive a system of diffusive equations from the ABMs (equations (LABEL:eq:PDE_0)-(LABEL:eq:PDE_1)-(LABEL:eq:PDE_2) and (LABEL:eq:PDE_3)). Moreover, our observations reveal some unusual phenomena (density spikes, anisotropy and aggregation) caused by the interplay of persistence and volume exclusion.

Despite the new diffusive PDEs correctly predicting the macroscopic behaviour of the ABM for a wide range of parameters, (including large values of reorienting rate) our findings highlight unrealistic behaviours for certain choices of parameters. Specifically, when the agent velocity is large and we implement one of the two more complex forms of agent interaction (types 2-3), the density profile of the ABM presents regular peaks (spikes) in density which are not captured by the corresponding continuum models. We also give evidence of an inherent anisotropy that occurs when we implement persistence in an on-lattice context in two or more dimensions Thompson et al. (2011). This phenomenon represents a problem when applying the model to experimental data, for which isotropy is usually a natural feature. One possibility for reducing the scale of this issue and still obtaining the macroscopic description would be to work on an hexagonal lattice. This would increase the number of preferential directions from four to six, making the anisotropy less evident, although not completely removing it. Alternatively, one could allow cells to move in diagonal directions as in Thompson et al. (2011) which may also serve to mitigate, but not completely remove anisotropy. Off-lattice models should have the advantage that they are not afflicted by anisotropy. However, the derivation of a corresponding macroscopic description becomes more complicated and sometimes intractable. Therefore, the problem of modelling persistence of motion in an on-lattice context at multiple scales, without incurring anisotropy in the lattice directions, remains an interesting challenge for future research.

Finally, the other main achievement of this work is that the new continuum approximation that we propose is capable of qualitatively reproducing the spontaneous aggregation driven by persistence and volume exclusion Slowman et al. (2016); Thompson et al. (2011); Sepúlveda and Soto (2016); Soto and Golestanian (2014); Peruani et al. (2006); Levis and Berthier (2014); Redner et al. (2013); Bialké et al. (2013). To our knowledge, this is the first time that such behaviour in the ABM has been replicated at the macroscopic level. The unintuitive consequence is that, in the case of strong persistence, the process of cell dispersion is initially slowed down by the aggregation phenomenon which constrains some of the cells in the internal region around the initial condition. Once the aggregates dissolve, the agents’ dispersion is effectively faster than the normal diffusion. When such aggregation occurs, the agreement between the models at the two scales is qualitative rather than quantitative. More work might be done in order to recover a better agreement. In particular, this might be achieved by including a higher order of spatial correlations in the continuous model as in Markham et al. (2013a, b), however this remains an open question.

The authors would like to thank the CMB/CNCB preprint club for constructive and helpful comments on a preprint of this paper.