Stochastic and deterministic modelling of cell migration
Abstract
Mathematical models are vital interpretive and predictive tools used to assist in the understanding of cell migration. There are typically two approaches to modelling cell migration: either microscale, discrete or macroscale, continuum. The discrete approach, using agentbased models (ABMs), is typically stochastic and accounts for properties at the cellscale. Conversely, the continuum approach, in which cell density is often modelled as a system of deterministic partial differential equations (PDEs), provides a global description of the migration at the population level. Deterministic models have the advantage that they are generally more amenable to mathematical analysis. They can lead to significant insights for situations in which the system comprises a large number of cells, at which point simulating a stochastic ABM becomes computationally expensive. However, finding an appropriate continuum model to describe the collective behaviour of a system of individual cells can be a difficult task. Deterministic models are often specified on a phenomenological basis, which reduces their predictive power. Stochastic ABMs have advantages over their deterministic continuum counterparts. In particular, ABMs can represent individuallevel behaviours (such as cell proliferation and cellcell interaction) appropriately and are amenable to direct parameterisation using experimental data. It is essential, therefore, to establish direct connections between stochastic microscale behaviours and deterministic macroscale dynamics.
In this Chapter we describe how, in some situations, these two distinct modelling approaches can be unified into a discretecontinuum equivalence framework. We carry out detailed examinations of a range of fundamental models of cell movement in one dimension. We then extend the discussion to more general models, which focus on incorporating other important factors that affect the migration of cells including cell proliferation and cellcell interactions. We provide an overview of some of the more recent advances in this field and we point out some of the relevant questions that remain unanswered.
Keywords: Cell migration, discretecontinuum equivalence, agentbased models, partial differential equation, collective behaviour.
1 Introduction
The process of cell migration plays an essential role in several developmental and pathological mechanisms. For example, cell migration orchestrates morphogenesis throughout the development of the embryo [Gilbert, 2003, Keller, 2005], and plays a crucial role in woundhealing [Maini et al., 2004, Deng et al., 2006] and immune responses [Madri and Graesser, 2000]. Migration also contributes to many pathological processes, including vascular disease [Raines, 2000] and cancer [Hanahan and Weinberg, 2000].
Many of these phenomena are highly complex and the collective behaviour of a set of interacting cells. At the individualcelllevel, a variety of mechanisms can be involved, including cellcell adhesion [Niessen, 2007, Trepat et al., 2009], attraction [Yamanaka and Kondo, 2014] and repulsion [CarmonaFontaine et al., 2008]. In other cases, cells can respond to a chemical gradient which regulates and guides their motility (chemotaxis) [Ward et al., 2003, Keynes and Cook, 1992]. In addition, when the migration occurs over a sufficiently long time, cell proliferation and death can also play important roles in the process [Mort et al., 2016]. An understanding the impact that these (and other) individualcell mechanisms have on global collective migration is important, since it can illuminate the origin of major diseases and suggest effective therapeutic approaches.
Over the past few decades, great progress has been made in understanding many aspects of cell migration [Ridley et al., 2003, Friedl and Gilmour, 2009]. However, we still lack a proper understanding of many of its underlying mechanisms. In particular, there are aspects which remain impenetrable to experimental biologists [Staton et al., 2004]. One of the major issues when studying cell migration is obtaining and interpreting experimental data. For example, a comprehensive controlled experiment in vivo or a realistic in vitro set up can be difficult or impossible to obtain. In other cases, it may be extremely complicated to investigate the microscopic origin of a macroscopic phenomenon, given the number of different actions that a single cell can perform [Westermann et al., 2003, Dworkin and Kaiser, 1985, Trepat et al., 2009, Tambe et al., 2011]. In this context, mathematical modelling has become a necessary interpretive and predictive tool to assist in the understanding of such complex phenomena [Noble, 2002, Simpson et al., 2006, Maini et al., 2004].
Mathematical models have been employed as tools for validation of the experimentally generated hypotheses. Moreover, due to the advances in computation of the last few decades, mathematical models are becoming more detailed and accurately parametrised which makes them capable of generating experimentally testable hypotheses and exploring new questions which are still not approachable from an experimental prospective [Tomlin and Axelrod, 2007]. For example, the use of stochastic mathematical models facilitates the investigation of biological systems in which randomness plays an important role and for which viewing only a single instance of the evolution of a system can be inconclusive [Lee and Wolgemuth, 2011]. In a mathematical modelling framework, the system can be repeatedly simulated using the same deterministic or randomised initial conditions. Equally, repeats with slightly altered initial conditions or parameters are useful for quantifying the sensitivity of systems in which a small fluctuation in the initial state of the system can have a significant effect on the outcome of the experiment, or in understanding the sensitivity of the system to certain parameter choices, respectively. In addition, the outcome of a mathematical simulation can be examined easily since all the variables of the system are explicitly accessible [Flaherty et al., 2007].
There are typically two approaches to modelling cell migration, depending on the scale of interest. At the level of an individual cell, stochastic, discrete agentbased models (ABMs) are popular. Each cell is modelled as a single individual (agent) with its own rules of movement, proliferation and death. Alternatively, to model collective cell migration at the population level, a deterministic, continuum partial differential equation (PDE) for the population density is normally used^{2}^{2}2Although deterministic ABMs [Kurhekar and Deshpande, 2015] and stochastic continuum models [Schienbein and Gruler, 1993, Dickinson and Tranquillo, 1993] have also be considered, the most common pairing in the literature involves stochastic ABMs and populationbased deterministic models. Therefore, these are the two modelling subtypes that we will consider in this review.. These two modelling paradigms have complementary advantages and disadvantages which we summarise in Table 1.
Generally speaking, ABMs are attractive because their macroscale behaviour is completely selfinduced, rather than being superimposed on a phenomenological basis [BenJacob et al., 2000, Shapiro, 1988]. This is particularly interesting when a complex structure emerges at the populationlevel [Othmer and Stevens, 1997, Thompson et al., 2011, 2012]. In fact, although the agents represent the driving force of such structure, typically they are unaware, as individuals, of the macroscopic configuration of the system, since their behaviour is dictated only by their local environment.
Another advantage of ABMs is that their formulation is generally more intuitive than continuum models. Each cell’s actions can be easily incorporated in the model by implementing appropriate rules for the corresponding agent, which mimic specific known biological behaviours. The parameters regulating of each of these rules can then be inferred directly from real observations which makes ABMs more easily relatable to experimental data.
Mort et al. [2016], for example, used a simple ABM to study the migration of mouse melanoblasts, the embryonic precursors of melanocytes, responsible for pigmentation (see Figure 1). The authors studied the effects of a family of mutations to the receptor tyrosine kinase, Kit. These mutations affect the success of the colonisation of the growing epidermis by the melanoblasts, leading to unpigmented regions of hair and skin. Mort et al. employed an ABM to study the interplay of movement and proliferation of melanoblasts and they carried out a statistical analysis on single cell trajectories in order to parametrise the model against experimental data. By simulating their ABM and comparing the results with experimental data, Mort et al. were able to show that bellyspot formation in Kit mutants is likely to be induced by a reduction in proliferation rate, rather than motility, a result which was contrary to the received wisdom in the experimental literature.
Deterministic models also have their advantages. For example, when the population size becomes large, a PDE model is generally desirable, since there exist a wide range of welldeveloped numerical tools for their rapid solution. This is in direct contrast to agentbased models which are typically coded up ad hoc and require multiple computationally intensive repeats in order to gather reliable ensemble statistics. Moreover, continuum models have the advantage that they are generally more amenable to mathematical analysis than ABMs, and this analysis can often lead to significant and general insights. For example, PDEs can be used to carry out stability analysis to determine the conditions which lead to pattern formation [Anguige and Schmeiser, 2009]. In other scenarios, one can use travelling wave analysis to obtain expressions for the speed of the cell invasion in terms of the model parameters [Murray, 2007].
Depending on the biological questions of relevance and available experimental data, either deterministiccontinuum or stochasticdiscrete models (or some combination of both) may be appropriate. Ideally, the individuallevel representation can be used to parametrise the model, and the continuum leveldescription should link the populationlevel results back to the parameters of interest.
The problem of connecting the parameters of the individuallevel model to those of a representative populationlevel description represents, therefore, a crucial step of the multiscale modelling process.
Advantages  Disadvantages  

Discrete/Stochastic 



Continuum/Deterministic 


In this Chapter, we provide an overview of a range of techniques which can be used to connect ABMs of cell migration to macroscopic PDEs for average cell density. We review a series of stochastic and deterministic models which are capable of reproducing some of the key features of the behaviour of cells and describe how these two modelling regimes can be linked to form a multiscale mathematical framework.
In Section 2.1 we carry out a detailed derivation of deterministic models from ABMs which focus on cell movement. We present derivations in two different scenarios, depending on whether the spatial domain in which the ABM is defined is partitioned into a finite lattice (the onlattice case) or not (the offlattice case). In each of these situations, we first study the case in which cells move independently of others (noninteracting cells) and then we introduce the ability of cells to sense the occupancy of neighbouring regions of space and to avoid overlapping with other cells (interacting cells). All the derivations in Section 2.1 are carried out in one dimension and assuming a meanfield moment closure approximation. We discuss generalisations to higher dimensions and other closure approximations in Sections 2.2 and 2.3, respectively.
The remaining part of the chapter is devoted to reviewing a series of biologically relevant features that can be incorporated in the models of cell migration. In each case we highlight the implications of introducing a given behaviour both at the individual and populationlevel. More precisely, in Section 3.1 we present models which incorporate the ability of cells to reproduce by division. In Section 3.2.1 we present models in which cells can interact indirectly through an external signal, e.g. chemotaxis or slime following. Some examples of direct forms of cellcell interactions, such as adhesionrepulsion, pushing and pulling, are discussed in Section 3.2. In Section 3.3 we review a series of recent papers which model cells migrating on a growing domain. Finally, we discuss how to derive a macroscopic limit for ABMs which are not based on simple random walks and which are capable of representing persistence of motion, in Section 3.4. We conclude the Chapter with a short discussion and final considerations in Section 4.
2 Cell motility
In this Section we present a detailed description of the basic models of cell motility, which represent the fundamental basis for the majority of the spatially extended representations of the remaining part of the Chapter. We use the case of these models to illustrate standard approaches to deriving diffusive, deterministic, continuum representation from ABMs, based on occupancy master equations. We carry out the explicit derivations for onedimensional versions of the models in Section 2.1, and we discuss the generalisation to higher dimensions in Section 2.2. We consider two distinct cases, depending on whether the motility mechanism of the ABM is implemented on or offlattice. In each case, we focus on two variants of the model, both with and without crowding effects. We explain the standard techniques for deriving the corresponding deterministic descriptions at the populationlevel in each case. Note that, throughout the Chapter, we adapt the notation of models taken from the literature for consistency.
2.1 Connecting stochastic and deterministic models of cell movement
2.1.1 Onlattice models
Broadly speaking, onlattice ABMs can be classified as cellular automata in which a set of agents occupy some or all sites of a lattice. In general, these agents have a number of state variables associated with them and a set of rules prescribing the evolution of their state and position. There exists a great variety of forms of ABMs. The appropriateness of each representation depends on the phenomenon that is being modelled. For example, cellular Potts models have been used by Turner and Sherratt [2002] and Turner et al. [2004] in the context of cancer modelling and by Graner and Glazier [1992] in the context of cell sorting via differential adhesion. Othmer and Stevens [1997] modelled bacterial aggregation by using a positionjump process and Painter et al. [1999] adopted a similar approach to study the interplay of chemotaxis and volume exclusion. One of the advantages of using these onlattice models is that their formulation and analysis tend to be straightforward in comparison to their offlattice counterparts. Moreover, the presence of the grid considerably reduces the computational cost of simulations when a large number of agents is involved. However, since in the majority of scenarios the assumption that cells move on a discrete grid is not appropriate, the implementation of cell behaviours is usually phenomenological.
Consider a onedimensional domain, with periodic boundary conditions. An ABM on the domain comprises a set of agents positioned in and a range of stochastic rules which govern their evolution in time. The ABMs that we consider throughout this Section are all based on continuoustime positionjump processes [Othmer et al., 1988, Othmer and Stevens, 1997]. This means that the position of the agents undergoes series of sequential Markovian jumps i.e., at any given time, the evolution of the position of each agent depends only on its current position. An alternative approach is to use velocityjump processes, in which the Markov property applies to the velocity of the agents, instead of their position. We discuss this approach in Section 3.4.
Discretetime ABMs are popular in the literature [Simpson et al., 2007, 2009, Treloar et al., 2011, 2013], although, for the examples presented here, we treat time as a continuous variable^{3}^{3}3In general, when the time discretisation step is sufficiently small, the behaviour of the discretetime ABMs is similar to its continuoustime counterpart. [Othmer et al., 1988, Othmer and Stevens, 1997]. As time evolves, agents can attempt movement events which occur as a Poisson process with rate . In other words, each agent attempts to move after an independent exponentially distributed waiting time with parameter . When such attempts take place, we say that the agent has been selected to attempt a movement.
Consider the onlattice scenario, in which the domain is partitioned in intervals, each of length , whose centres are denoted by , respectively (see Figure 2 for a schematic illustration). We denote the position of the th agent at time as , hence for every and .
Noninteracting cells undergoing unbiased movement.
In this Section, we consider the basic case of noninteracting agents. In particular, the movement rates of each agent are independent of the other agents’ positions and unbiased. If, at time , an agent, , is selected to move, it attempts to move with equal probability, , to either one of its adjacent sites, . Equivalently, we can say that each agent moves to the right and to the left with rates , respectively. Notice that multiple agents can occupy the same lattice site. We aim to obtain a deterministic representation of the mean evolution of the agent density at a given time . Therefore, we assume the model is simulated until time for a large number, , of independent realisations, with each realisation identically prepared. In Figure 3 (a) we consider 30 onedimensional lattices with and . In each onedimensional lattice we populated the 20 central sites, , at random with agents on average (with multiple agents per site possible). In Figures 3 (c),(e) and (g) we show three snapshots of these 30 identically prepared simulations of the ABM, as time evolves. We denote by the number of agents which lie in the interval at time of the th repeat simulation, with . Namely
(1) 
If , we say that the site is empty and if we say it is occupied. We define the mean occupancy of site at time , averaged over the number of realisations, as
(2) 
By considering all the possible agents movements we can write down a conservation law for the average occupancy at time , where is a sufficiently small that the probability that two or more movements take place in the interval is . This reads
(3) 
The righthand side of equation (3) comprises three parts: the first term, which accounts for the occupancy of site at time , and two terms which determine the expected loss of average occupancy due to agents moving out of site and the gain due to agents moving into site , respectively.

If we rearrange equation (3), divide through by and take the limit as , we obtain the Kolmogorov equation (or continuoustime master equation) of the process which reads
(4) 
Now we Taylor expand the terms about the point to the second order, i.e. we use the following approximations:
(5a)  
(5b) 
By substituting equations (5) into equation (4) and taking the limit as while holding constant, we recover the diffusion equation for the continuous approximation of the average occupancy function, :
(6) 
where is the diffusion coefficient defined as
(7) 
We refer to Figure 3 (i) for a comparison between the average occupancy of the ABM and the diffusion equation (6). The results confirm the good agreement between the two models.
The derivation of equation (6) from a simple random walk, as outlined above, is a well known result [Murray, 2007, Deutsch and Dormann, 2007, Codling et al., 2008]. Alternatively, we could have carried out the derivation for the probability density function, , of finding a single agent at position at time . This would result in a macroscopic PDE of the same form as equation (6). The equation for can be obtained by dividing both sides of equation (6) by the total number of agents, .
The diffusion equation (6) is a classical PDE, also known as Fick’s diffusion equation or the heat equation depending on the application. Extensive discussions on this type of equation can be found in Crank [1979] or Welty et al. [2009]. The existence of an analytic solution of equation (6) depends on the imposed initial and boundary conditions. For example, if we assume an infinite domain and that all agents are initialised in the same position:
(8) 
equation (6) admits the fundamental solution given by
(9) 
In the basic model outlined above, agents can only move to their nearest neighbour sites. A generalisation of this model, in which agents have the ability to perform nonlocal jumps is studied by Taylor et al. [2015a]. Specifically, in the model of Taylor et al. [2015a] agents are allowed to jump up to sites away from their current position. For , each length jump occurs with rate in the right and leftdirections, respectively. When the motility is unbiased, i.e. for every , the authors show that the behaviour of their agents in the continuum limit evolves according to the diffusion equation (6), but with diffusion coefficient given by
(10) 
This implies that, when the following condition is satisfied,
(11) 
the local ABM and the nonlocal ABM are described by the same macroscopic PDE. There are an infinite number of possible combinations of s which satisfy condition (11). However, Taylor et al. [2015a] suggest that choosing
(12) 
is particularly appropriate, since it preserves the wellknown property of diffusive processes, that the meansquared displacement scales linearly with time [Codling et al., 2008, Othmer et al., 1988]. A comparison of the simulations of the two AMBs confirms the accuracy with which the nonlocal models typically correspond to their local equivalent. The results also reveal a significant reduction in the average simulation time of the nonlocal ABM compared to the local ABM. This timesaving potential, given that condition (11) is satisfied, is highlighted by Taylor et al. as an important application of the nonlocal models in order to reduce the computational cost of stochastic simulations.
However, it should be noted that, when dealing with steep gradients in agent numbers, the nonlocal model loses accuracy in comparison to the local model. This inaccuracy stem from the choice of the transitional rates (12) which match the terms of the local model only up to the second order of the expansion. In order to address this issue, the authors propose a spatially extended hybrid method in which regions containing steep gradients in agent numbers are dealt with using the local representation and regions in which gradients are more shallow using the nonlocal representation. This allows the acceleration of simulations afforded by the nonlocal representation, whilst maintaining the accuracy associated with the local method.
In the second part of the paper, Taylor et al. derive a general class of boundary conditions for local and nonlocal ABMs, corresponding to classical boundary conditions in the deterministic, continuum, diffusive limit.
For a general nonlocal ABM, the authors study a class of firstorder reactive boundary condition known as Robin conditions which, for the leftboundary, can be written as
(13) 
Here corresponds to a purely reflective boundary and corresponds to a purely absorbing boundary. In order to obtain the corresponding stochastic boundary condition at the individuallevel, Taylor et al. allow agents which attempt length jumps that result in them hitting the boundary, to be absorbed (i.e. removed from the domain) with absorption probability . If no absorption occurs, the agents reach the boundary and are then reflected in the opposite direction for the remaining number of steps of the jump. By taking a diffusive limit from the corresponding occupancy master equation, the authors determine the expression of the absorption rates in terms of the parameter, , of the Robin boundary condition. The absorption rates are given by (13):
(14) 
where is the lattice step, is the maximum jump length and is defined as in equation (10). In particular, when a local ABM is considered, the absorption rate, , for agents at the site adjacent to the boundary is given by
Finally, Taylor et al. [2015a] extend their study to the case in which Neumann boundary conditions are imposed at the population level. For the left boundary this condition is given as
(15) 
where represents the influx into the system at this boundary. The corresponding ABM implementation requires new agents to enter the domain. For a general nonlocal model with maximum jump length , these agents are positioned in the nearest sites to the boundary. By postulating a discrete master equation, which respects conservation of the total influx, , and taking a diffusive limit as , the authors derive an expression for the rate of introduction of agents into the th nearest site:
(16) 
which in the local case (), becomes
Excluding cells
We now incorporate volume exclusion in the model following the approach of Simpson et al. [2009]. We modify the basic ABM of the previous Section, by introducing a specific form of agentagent interaction which prevents two or more agents from occupying the same lattice site. We initialise the domain by populating the central interval, with excluding agents, i.e. with the property that any given site can be occupied by at most one agent. In the example in Figure 3 (b) in each row we populated sites with agents, on average, with a maximum of one agent per site. When an agent, , is selected for a movement event, one of its two neighbouring sites, , is selected with equal probability, . In order for the movement to be successful, however, the new selected site has to be empty. When the selected site is occupied, the movement is aborted and the selected agent remains in its position. This movement rule prevents agents from moving into occupied sites. Therefore, the excluding property of the initial condition is preserved as the time evolves. In particular, this implies that , defined in equation (1), takes only two values, and . See Figures 3 (d),(f) and (h) for snapshots of simulations of the ABM.
For a given a realisation of the model, , we make the usual momentclosure approximation that the occupancies of different sites, and with , are independent, i.e.:
(17) 
for every . Equation (17) is know as the meanfield assumption, and it provides a good approximation in many scenarios, however it becomes invalid when spatial correlations play an important role in evolution, for example when proliferation or attractive forces between agents are involved (see Sections 3.1 and 3.2, respectively). We discuss the validity of this approximation and possible alternative approaches to the mean field assumption in Section 2.3.
By assuming the independence captured by equation (17), we can write down the transition rates from a given site as
Note that the terms keep into account the possibility of abortion of the movement due to volume exclusion. Hence the occupancy master equation at position reads:
If we divide by , rearrange and let , we obtain
(18) 
for . By summing equations (18) over each repeat and dividing by the total number of realisations, , we recover equation (4) for the average occupancy, . Remarkably, we obtain exactly the same macroscopic representation as in the case of noninteracting agents, which is given by the canonical diffusion equation (6). In Figure 3 (j) we show a comparison between the average occupancy of the ABM with excluding agents and the corresponding diffusion equation at increasing times.
Notice that this inclusion of volume exclusion on a lattice does not affect the populationlevel dynamics. In other words, since the diffusion equation (6) arises as the limit equation of the standard random walks of noninteracting agents, the effect of the local interaction between the agent in the simple exclusion process disappears as we let . Since agents are indistinguishable, the situation in which two agents occupy adjacent positions and block each other’s movement is equivalent to the scenario in which the two neighbouring agents swap their positions in the noninteracting random walk. If multispecies agents are considered, such equivalence no longer holds and the exclusion property leads to different continuous equations [Simpson et al., 2009].
Taylor et al. [2015b, 2016] study alternative approaches to implementing volume exclusion in compartmentbased models at different spatial scales. They consider a coarsegrained representation of volume exclusion (as opposed to finegrained representation described above, in which at most one agent can occupy a given site) in which finegrained sites are amalgamated together into a single coarsegrained site with capacity and length . This coarsegrained representation is referred to as a partiallyexcluding ABM. See the schematic in Figure 4 for an illustration. Agents can perform jumps between neighbouring compartments with a rates proportional to and which scales linearly with the proportion of available space in the target compartment.
Taylor et al. [2015b] consider a uniform regular lattice. By comparing the fullyexcluding model, , and the partiallyexcluding model, , the authors show that the mean and variance of the number of agents in each compartment is the same in both representations. In other words, it is possible to provide a consistent description of the effects of volume exclusion across different spatial scales with these coarse and fine representations. While the usage of a coarsegrained, partiallyexcluding approach leads to a significant time saving in their example simulations, it also requires movement events to occur across a wide range of spatial and temporal scales [Walpole et al., 2013].
In a more recent work, Taylor et al. [2016] extend the study of these partiallyexcluding models to nonuniform lattices in which the sites’ carrying capacities can vary across the domain. In particular, they develop a set of hybrid methods which allow the interfacing of regions of partiallyexcluding sites with regions of fullyexcluding sites. The advantage of these hybrid methods is that they allow for the study of complicated scenarios, in which the accuracy of the fully excluding model is required in some regions of space, but the partiallyexcluding model can be exploited in other regions, allowing a considerable reduction of computational cost in comparison to the ubiquitous fullyexcluding model [Taylor et al., 2016].
Simpson et al. [2009] study other important extensions of the fundamental volumeexcluding ABM outlined towards the start of this Section. In particular, the authors incorporate the possibility of having multiple subpopulations of agents with different motility parameters and a deterministic directional bias.
The ABM of Simpson et al. [2009] is defined in discrete time and on a twodimensional lattice with lattice step . An exclusion property is implemented, requiring that a lattice site can be occupied by at most one agent at the time. Agents move according to a biased random walk with the bias modulated by a parameter .
In the first part of the paper, a single population of identical cells is considered. Simpson et al. [2009] derived an advectiondiffusion continuum approximation of the model by writing down the occupancy master equation and letting the bias parameter scale with the spatial step . The simulations of the ABM are averaged over the columns of the lattice and over many realisations. The resulting density profiles are compared with continuum descriptions showing good agreement with the corresponding PDE.
In the second part of the paper, the model is extended to describe the migration of subpopulations or species of cells. The cells of each species move according to a biased random walk, with the motility rate and bias intensity , with , depending on the species. Using a similar derivation as in the onespecies case, the authors derive a system of advectiondiffusion equations describing the evolution of the occupancy of the different species, . When the bias is turned off, the system reads
(19) 
where
In this case, the exclusion property leads to nonlinear diffusivity for the individual species. However, in the case in which all species have the same motility parameters as each other, by summing all the equations, unsurprisingly, simple diffusion is recovered for the total population.
To compare the two levels of description, the authors consider a population of agents formed by two species initialised in adjacent regions with different initial densities. The results highlight a spontaneous aggregation in one of the species’ density profiles in the continuous model. From this observation, Simpson et al. [2009] conclude that the single species densities do not obey any maximum principle, i.e. the monotonicity of the density profile is not preserved in time, although this is true for the total population.
2.1.2 Offlattice models
In offlattice ABMs the positions of the agents are represented in continuous space. This improvement adds realism to the model, since the movement of real cells is not constrained to a discrete grid. Moreover, the continuous framework introduces a larger variety of possible actions which cells can perform. For example, when modelling the migration of cells in two dimensions, the offlattice framework allows both the distance and the direction of the movement to be a continuous variable, rather than being restricted to a discrete set of values, as in the onlattice case. This extension is consistent with many biological observations of cell migration in which cells are not restricted to a lattice [Stokes and Lauffenburger, 1991, Plank and Sleeman, 2004]. However it has the disadvantage that it makes the mathematical analysis more complicated and sometimes intractable.
Noninteracting cells undergoing unbiased movement.
We consider an ABM on a onedimensional domain, , as in Figure 5, with periodic boundary conditions. Each agent, , is defined by its position at time , , and it occupies the interval , where is the analogue of the cell’s radius. The model is defined in continuous time and agents perform unbiased jumps of fixed distance in one of the two directions. The rate and the distance of these jumps are denoted by and respectively. See the schematics in Figure 5 for an illustration of the model. A set of agents are initially located uniformly at random in the interval (see Figure 6 (a)).
For now, we assume agents move independently of other agents’ positions, as in Othmer et al. [1988]. Consequently, a given point, , can be occupied by more than one agent simultaneously. In this case, obtaining a macroscopic description of the model can be achieved by employing a similar method to the corresponding onlattice case (Section 2.1.1). We consider identically prepared simulations of the model and aim to write down the master equation for the average occupancy of position at time , , which is defined as
(20) 
where
(21) 
The master equation for then reads
(22) 
where is chosen sufficiently small that at most one movement event can take place in . By Taylor expanding the terms to the second order about the point and taking the limit , while keeping a non zero constant, we recover the diffusion equation (6), with the diffusion coefficient defined as
(23) 
In other words, when agents are not interacting, the offlattice framework of the ABM does not change the resulting populationlevel representation which was obtained for the onlattice case. In Figure 6 (i) we show a comparison of the average agent density and the corresponding diffusion equation as time evolves. The results confirm a good agreement between the stochastic and the deterministic models.

Volumeexcluding cells
Incorporating volume exclusion in an onlattice model is a natural extension of the simple multiple occupancy model. However, in an offlattice framework there are several ways that the effects of volume exclusion can be incorporated. In this section we follow the approach of Dyson et al. [2012], however we highlight that other approaches can lead to slightly different results both at the individual and populationlevels.
Dyson et al. [2012] consider an exclusion mechanism in which an attempted move is aborted if it would lead to the overlap of two agents, i.e. the corresponding centres are closer than . A schematic in Figure 7 shows an example in which an agent (blue) attempts to move in the rightwards direction and an examples of an agents (grey) which would obstruct the movement. In order to write down the occupancy master equation for the average number of agents, we need to compute the transition rates, . In particular, we need to determine the probabilities that an agent located at position at time , moves to the right and leftdirections, at time . By using the continuous form of the meanfield assumption (17), we can write
(24a)  
(24b) 
where and are the number agents in the exclusion zones and , respectively, which impede the movement. See the interval highlighted in pink in Figure 7 for an illustration right exclusion zone.
Hence the master equation for the occupancy, , reads
(25) 
where is defined as in equation (21). In order to compute , we reduce to the case in which there is at most one other agent obstructing the movement by assuming . If we consider a number of agents, , sufficiently large and we use the continuous version of the meanfield approximation (17), we can write
(26a)  
(26b) 
Notice that, due to the assumption on , the two integrals on the righthand side of equations (26) assume values in and so the two transition probabilities defined in equations (24) are meaningful. We can Taylor expand and to second order to obtain
(27a)  
(27b) 
By substituting equations (27) into equation (25) and taking the limit as , we obtain
(28) 
for . We can now sum equations (28) over all values of and divide by the total number of realisations, , to obtain an equivalent equation for the average occupancy ) using the assumption that all agents are identically uniformly distributed initially we find
(29) 
The dependence on the agents’ size can be explained by noting that larger agents will collide more often. Dyson et al. [2012] identify that when is large compared to , the diffusion coefficient decreases to the point at which it may be negative, which might suggest the occurrence of cell aggregation. It is also important to notice that, if the agents are initialised on a lattice with step , and we choose , the ABM is equivalent to the onlattice ABM with excluding agents. However, for such choice of , the derivation of equation (29) breaks down, which explains why simply setting in equation (29) does not recover the simple diffusion equation as might be expected. By taking the limit as , while keeping a nonzero constant, we arrive at a nonlinear diffusion equation:
(30) 
with
and is defined as in equation (23).
Notice that the exclusion property results in a nonlinearity in the diffusion coefficient of the populationlevel equation. In particular, the term in equation (30) leads to faster diffusion where the average occupancy is large and slower diffusion where the average occupancy decreases to zero, in which case the diffusion coefficient reaches its minimum value, . In Figure 6 (j) we compare the average agent density with the numerical solution of equation (29) at four subsequent times. The results confirm the good agreement between the ABM and the corresponding continuum equation.
2.2 Higher dimensions
Although the detailed derivations of the previous Section are carried out for a onedimensional interval domain, discretecontinuum equivalence frameworks can be extended to higher dimensions.
In fact, for models which do not account for crowding effects, either on or offlattice, the resulting macroscopic description can be obtained in a similar manner to the onedimensional case [Othmer et al., 1988, Deutsch and Dormann, 2007, Codling et al., 2008]. For the regularsquarelattice and the offlattice models, the resulting deterministic description is the natural generalisation of equation (6) which is given by
(31) 
where represents the Laplacian operator and the diffusion coefficient is given by
where is the dimension. Equation (31) is an isotropic PDE, i.e. it is not biased in any spatial direction. That such isotropic equation can be derived from an ABM which is defined on a regular, anisotropic lattice is perhaps surprising. In other words, the intrinsic individuallevel anisotropy of the ABM in the lattice directions vanishes in the diffusive macroscopic description [Othmer et al., 1988, Deutsch and Dormann, 2007, Codling et al., 2008]. This same property has been demonstrated to not hold for models of cell movement which are based on velocityjump processes [Gavagnin and Yates, 2018]. We refer the reader to Section 3.4 for a more detailed discussion of such models.
Incorporating volume exclusion in higher dimensions does not lead to substantial changes for the onlattice models in which case equation (31) still holds [Simpson et al., 2009]. However, higher dimensions significantly increase the complexity of offlattice ABMs with crowding effects. For example, Dyson and Baker [2014] study an extension of their previous onedimensional model in two and three dimensions. In their ABMs, cells are represented as circular or spherical agents of radius . To include volume exclusion, agents movements which would lead to an overlap of agents are aborted. In Figure 8 we reproduce a schematic of the volume exclusion property for the ABM of Dyson and Baker [2014] in two dimensions. Notice that, to compute the probability of finding obstructing agents, it is necessary to integrate the occupancy function over the grey shaded region, , of Figure 8. This makes the computation complicated and intractable from a mathematical point of view. To overcome this problem, Dyson and Baker [2014] suggest a simplification of the calculation by extending the integral to the two blue regions, . Clearly, if the distance of each jump, , is sufficiently small, this is a reasonable approximation which significantly simplifies the analytical calculation.
By adopting this simplification, Dyson and Baker [2014] carry out the derivation of diffusive PDEs for the average agent occupancy, leading to equations of the form of (29). In all cases, the nonlinear diffusion coefficients are increasing functions of the agents’ radii, , and decreasing functions of the jump distance, .
Finally, the authors consider a system comprising multiple species of agents in which the size and the movement rates of the agents depend on the species to which they belong. The authors consider an example with two species, one comprising agents that are large and slow moving, and the other smaller but quicker. A system of PDEs for the densities of the two species is derived and the results show that the species with smaller agents is less affected by volume exclusion. Interestingly, the effect is reduced in the area where the two species coexist, since the majority of the area is occupied by smaller agents.
2.3 Higherorder closure approximations
In general, when agentagent interactions are included in an ABM, as in the example of volumeexclusion above, the positions of the agents are not independent. When this happens, the evolution of the average agent density depends on the distribution of agent pairs, which also depends on the distribution of agents triples and so on. Therefore, if we aim to derive a deterministic representation of the model, we have to deal with an infinite system of unclosed equations for each of the spatial moments of the agents’ distribution. In order to overcome this problem, a moment closure approximation is necessary in order to make the system amenable to mathematical analysis.
The meanfield closure, given by equation (17), represents the easiest form of moment closure which is typically used [Simpson et al., 2009, Cheeseman et al., 2014, Dyson et al., 2012, Dyson and Baker, 2014]. Although such rudimentary approximation provides good results in a wide range of scenarios, when the spatial distribution of the agents is strongly correlated, the meanfield closure can lead to under or overestimations of the total agent density, resulting in a poor agreement between the stochastic ABM and its deterministic representation.
Baker and Simpson [2010] investigate the role of spatial correlation in exclusion processes and its effect on the agreement between the averaged agentbased dynamics and the meanfield approximation.
They consider an excluding ABM on two and threedimensional lattices. Agents can move to neighbouring lattice sites, with rate , and proliferate by placing daughter agents into neighbouring lattice sites, with rate . They derive a general set of master equations for the point distribution functions, , (of which the 1point distribution function is simply the density, the 2point distribution function is the pairwise occupancy etc). This infinite hierarchy of master equations is unclosed: the differential equations for the point distribution functions depend on the point distribution functions. Therefore, a closure approximation is required in order to solve for the lowerorder distribution functions. The authors compare the averaged density for the discrete model with the firstlevel moment closure (the meanfield approximation), in which neighbouring sites are assumed independent (see equation (17)). They also compare to a secondlevel moment closure (the Kirkwood superposition approximation), which takes into account pairwise spatial correlations. On the square lattice, the distance between two lattice sites increases irregularly (, , , ) and the neighbours for each site, therefore, have to be calculated separately for each distance. This fact means that the number of ODEs for the correlation functions becomes intractable very quickly. The system of ODEs, therefore, needs to be truncated at a maximum distance , beyond which the sites are considered independent. Different values of are compared and the results suggest that, in two dimensions, the system can be truncated at without losing accuracy. For the threedimensional case, the cutoff can be reduced to . The inclusion of correlations in the model, even if closed at level two, provides a significant improvement in the approximation to the agentbased model.
Baker and Simpson also investigate the effects of motility, birth and death events on spatial correlations. As the proliferation parameter increases with respect to the (fixed) motility parameter, spatial correlations play a more important role and the meanfield prediction appears to overestimate the growth of the population. This is due to the fact that cell motility can not break up the correlations caused by the appearance of new agents close to their parents sufficiently quickly. This leads to cluster formation which, in turn, reduces the number of successful proliferation events, slowing population growth. When agent death is included in the model, a counterintuitive effect appears. We may naively expect that deaths would decrease spatial correlation and so increase the agreement between the meanfield and the agentbased models. Instead the opposite happens. A high death rate has a similar effect to having a more sparsely populated initial seed. It provides opportunities for correlations to build up where previously sites occupancies were uncorrelated.
Markham et al. [2013a] continue the work of Baker and Simpson [2010] by deriving a deterministic, continuum analogue of a spatially extended, onlattice, agentbased model. As in the previous paper, they look at the influence of including spatial correlations in the continuum model and determine how this improves the agreement, in comparison with the meanfield model, which assumes absence of spatial correlations.
The ABM considered is the same as in Baker and Simpson [2010], i.e. a volumeexclusion process on a lattice on which individuals can move, proliferate and die with rates , and , respectively.
Baker and Simpson [2010], truncate the pairwise correlation functions at a maximum distance, , resulting in a system of ODEs. The aim of Markham et al. [2013a], however, is to derive a tractable deterministic PDE for the evolution of the pairwise spatial correlation function. Specifically, as the lattice step is small relative to the size of the domain, one can Taylor expand the correlation functions up to second order in and then move to radial coordinates. Finally, assuming that remains constant as , one obtains a reactiondiffusion PDE.
A trivial analysis shows that, as the motility parameter, , increases, the diffusion coefficient of the PDE for the correlation function increases. This agrees with the intuitive prediction that higher rates of movement break up clusters of agents more effectively. On the other hand, the reaction term of the PDE is negative and decreases as the proliferation rate, , increases. This means that increasing proliferation leads to a decrease in correlations. This fact appears to contradict the results of Baker and Simpson [2010]. Nevertheless, Markham et al. [2013a] provide an informal explanation for this phenomenon.
The results of spatially extended simulations show good agreement between the solution of the PDE and the discrete model. In particular, the PDE approximation behaves similarly to the ABM which, in the case of nonzero death rate, is significantly lower than the meanfield prediction. Moreover, an improvement in the agreement between the PDE and the ABM is obtained if the movement rate increases or if higherdimensional domains are considered.
Markham et al. [2013a] find that there exists a region of parameter space for which the deterministic models (either with or without correlations) can not replicate the ABM dynamics. This phenomenon corresponds to high values of the death rate which leads the population of the discrete model to eventually go extinct.
Markham et al. [2013b] extend their previous spatial correlation model [Markham et al., 2013a] to an ABM with heterogeneous agents. They provide other examples in which spatial correlations play a crucial role in predicting the populationlevel behaviour correctly. Agents are divided into species and can move, proliferate and die with rates , and respectively, where denotes the specie of the selected agent. With the same idea as their previous work, Markham et al. [2013b] focus on the evolution of the pairwise correlation, closing the system at the level above. Firstly they derive a set of ODEs for the system in which the pairwise correlation is divided into autocorrelation, i.e. between agents of the same species, and crosscorrelation, if the agents belong to different species. Subsequently they obtain a set of PDEs by Taylor expanding and taking the limit as the lattice step goes to . A key assumption throughout the paper, is that the pairwise correlation depends only on the distance between the agents, i.e. it is translationally invariant and isotropic. For this reason, the agents are always assumed to be spread uniformly across the domain initially.
The results for two species show a good agreement between the PDE model and the ABM in comparison to the agreement between the meanfield approximation and the ABM. The authors assume agent death is negligible and consider two species of cells with the same proliferation rates and different movement rates. While the logistic dynamics of the meanfield model predicts that the densities of the two species will converge to the same steady state, in the ABM, the species with the greater movement rate reaches a higher density at the equilibrium. They also show that with specific values of proliferation rates, the results can be the reversed, i.e. in the meanfield approximation, the density of one species at the equilibrium is greater than the other specie’s density, while the ABM predicts that, at equilibrium, the two species reach the same density. Remarkably, the PDE model incorporating agentagent correlation shows a good agreement with the ABM in all these scenarios.
3 Model extensions
The models summarised in the previous section, represent the fundamental basis for the vast majority of spatially extended models of cell migration. However, during the process of migration, cells perform a variety of other actions and interactions whose role can dramatically impact upon the macroscopic behaviour of the system [Niessen, 2007, CarmonaFontaine et al., 2008, Ward et al., 2003, Keynes and Cook, 1992, Trepat et al., 2009, Tambe et al., 2011]. Here we provide a brief overview of some of the most important extensions and modification of the standard ABMs and demonstrate their effects on the corresponding macroscopic models.
3.1 Cell Proliferation
In all the models considered so far, the ability of cells to divide and produce daughter cells was not included. In reality, in certain circumstances (e.g. tumour invasion and wound healing) the role of cell proliferation is crucial for the dynamics of the system [Sherratt and Chaplain, 2001, Maini et al., 2004].
Simpson et al. [2007] propose an ABM which takes into account cells’ motility and proliferation. The ABM is defined on a twodimensional lattice with a simple exclusion property, meaning that each lattice site can be occupied by at most one agent at the time. The model advances in discrete time. At each time step, every agent can move and proliferate with probabilities and , respectively. Agents move according to a simple, unbiased random walk to one of the their four nearest sites (von Neumann neighbours). When a proliferation event occurs, the proliferating agent moves to one of its eight surrounding sites (Moore neighbours) and the new offspring is displaced to the site diametrically opposite. If either of the two sites is already occupied, the proliferation event is aborted. Both motility and proliferation events are limited by a carrying capacity . In particular, if an agent attempts to move or proliferate into a site with a number of surrounding neighbours greater than , the event is aborted.
Firstly, the authors simulate an invasion wave and investigate some of the features of the ABM. They relate and compare the behaviour of their ABM with the traditional deterministic Fisher equation [Fisher, 1937]:
(32) 
where is the diffusion coefficient, is the growth rate and is the carrying capacity. Note that the average speed of invasion of the deterministic Fisher wave is known to be .
The results of the ABM simulations show that, consistent with the continuum equation, the wave speed of the invasion increases with motility parameter, , and proliferation probability, . However, the speed in the ABM is more sensitive to variation in proliferation than in motility, in contrast to the corresponding deterministic description which predicts equal sensitivity. As in the continuum setting, the speed of invasion is found to be independent of the carrying capacity.
To establish a connection between the ABM and the traditional continuum model given by equation (32), Simpson et al. compare the behaviour of the two levels of description for scenarios in which either proliferation or motility are active, but not both. When agents can only proliferate, the total density evolves in a logistic manner and the authors identify a linear relationship between the parameters of the agentbased and populationlevel models. When agents can move but not proliferate, the results show that the diffusivity of the model is a decreasing function of the background density. Nevertheless, a linear dependence is found between the motility probability of the ABM, , and the diffusion coefficient at zero agent density.
Finally, the ABM is used to investigate the individual cells’ trajectories within the invasion wave. The authors use synthetic data to obtain statistical information on the average direction of movement. The results compare the behaviour of the most advanced cell in the invasion with a second cell close behind the wave front. Despite the cells moving according to an unbiased random walk, the crowding effects inhibit the motility into highly populated regions and this induces a bias towards the direction of the invasion in both the cells considered. By using the individual agent trajectories, this bias can be quantified. Consistent with experimental data [Druckenbrod and Epstein, 2007], the observations on the ABM show a more evident bias for cells near the wave front in comparison to cells behind the wave front.
Cheeseman et al. [2014] continued the study of cell proliferation by modelling the spatial and temporal dynamics of different cellular lineages within an invasion wave. The authors defined an ABM model on a twodimensional lattice with a volume exclusion property. The model is initialised with a population of cells located on one end of the domain (see Figure 9 (a)). Agents are labelled according either to their generation number or the lineage to which they belong. The results for the ABM simulation show a clear spatial organisation in the distribution of the agent generation number. However, there is a large individual variability in the spatial distribution of a single agent’s linage tracing (see Figure 9). In order to reproduce the dynamics of the agent generation number, a set of nonlinear diffusion equations are derived from the ABM and are studied in one dimension. In particular, by using the same approach as Simpson et al. [2009], Cheeseman et al. derive a system of conservation of mass equations describing the density profile for each generation . In order to investigate the lineage tracings, they develop a GenerationDependent GaltonWatson (GDGW) process and they look at the Lorenz curves, which describe the proportional contributions of each of the initial seed cells’ lineages to the final population [Lorenz, 1905].
The results from agentbased and PDE models show that the invasion wave is composed of spatially regular and predictable generations. As the invasion occurs, the older generations reach a steady state behind the travelling wave. The structure is more apparent when averaging over more realisations of the ABM model. Indeed, both the mean and the variance of the generation number increase linearly with the distance from the location of the initial group of cells. Conversely, the individual agent lineages exhibit a clear asymmetry in their contribution to the final population. A small group of cells (superstars) generate the majority of the total population. The GDGW model shows a good agreement with the ABM results. This fact is interesting, since the GDGW process ignores some of the spatial and temporal correlations between individual lineages. Nevertheless this phenomenon can be partially explained by noticing that the competition between agents affects the generation densities, , which are used in the GaltonWatson process.
As highlighted in Section 2, choosing between on and offlattice frameworks can lead to differing macroscopic descriptions for excluding ABMs incorporating motility only. The same can be said for models of cell proliferation. Plank and Simpson [2012, 2013] provide a pioneering comparison between onlattice and offlattice models of proliferating and migrating cells.
The authors develop a new offlattice discretetime ABM for migration and proliferation of cells in a twodimensional domain. Agents are represented as incompressible circles of diameter and the total number of agents at time is denoted . At each time step, agents are selected independently and are given the chance to move with probability . A selected agent attempts to move a fixed distance, , in a given direction, , which is chosen uniformly at random in . To include crowding effects, the movement is aborted if any of the other agents lie within a distance of the line segment which connects the current location to its potential new location. Once all the motility events have been attempted, another agents are selected independently to attempt proliferation events with probability . A proliferating agent attempts to divide into two daughter agents whose positions are chosen diametrically opposite each other at distance from the original agent. The chosen proliferation event takes place only if it does not cause overlap between agents.
The authors derive a meanfield approximate ordinary differential equation for the spatially averaged agent density, . In doing this, they assume that the population of agents is homogeneous in the domain, which is known to be a poor assumption for large values of proliferation [Baker and Simpson, 2010, Markham et al., 2013b].
Plank and Simpson [2012] show a comparison between their offlattice model and a standard onlattice model [Simpson et al., 2010] for a spatially uniform initial condition. For small values of agent density, the absence of significant agentagent interactions leads to similar behaviours for the on and offlattice models. However, the two types of approach show a substantial difference for high densities. The onlattice model leads to a faster growth of the population density. The density eventually reaches the maximum value of unity when all the lattice sites are occupied. Conversely, it is impossible for agents in the offlattice model to reach the theoretical maximum agent density. This because the agents are not perfectly aligned, and at high densities become jammed in more realistic, yet irregular configurations. The appearance of a natural carrying capacity in the offlattice model, as opposed to the artificial value induced by a onlattice approach, suggests the offlattice model is a more suitable representation for biological applications.
Plank and Simpson [2013] further investigate the behaviour of their onlattice and offlattice models in scenarios in which the spatial variability of the cell density profile plays an important role. In particular, they initialise the domain by uniformly populating only the lefthand side region with a relatively low density. They then proceed to study the speed and the shape of the resulting invasion waves in the two models.
Their results highlight that the two models behave differently behind the invasion front. However, the models’ behaviours are similar at the leading edge. These findings are consistent with the previous observation, that the effects of crowding are more evident in the offlattice model than in the onlattice model [Plank and Simpson, 2012]. Specifically, in the invaded region, where the agent density is higher, the population size of the offlattice model approaches carrying capacity more slowly than the analogous onlattice model. However, at the front of the invasion, the low value of agent density makes the two models almost indistinguishable. This implies that, in the longterm, the speeds of the invasion fronts in the onlattice and offlattice models are the same.
Finally, Plank and Simpson [2013] carry out a leastsquares parameter estimation in order to fit the two continuum representations to density profile data simulated using the offlattice ABM. The solutions for both the onlattice and offlattice PDEs are in good agreement with the simulated data and it is hard to distinguish between the two fitted curves. The authors present this as an example that highlights the difficulty of choosing the correct mathematical framework when modelling real experimental data.
Typically, cell proliferation is represented in ABMs as a Poisson process with a certain rate , which means that each agent’s interdivision times are exponentially distributed [Mort et al., 2016, Simpson et al., 2007, Treloar et al., 2013, Turner et al., 2009]. One of the advantages of this approach is that, due to the memoryless property of the exponential distribution, the process can be efficiently simulated using the popular Gillespie algorithm [Gillespie, 1977].
However, assuming exponentially distributed cell interdivision times is not biologically realistic [Golubev, 2016]. In particular, the monotonicity of the exponential distribution implies that the most likely time for a cell to divide is immediately after its own creation. Recently, Yates et al. [2017] have proposed a novel model for incorporating nonexponentially distributed cellcycle times, based on a multistage scheme. The authors divide the cell cycle into stages. The waiting time for an agent to pass from the th phase to the th phase is exponentially distributed with rate . When an agent exits from the th stage, it divides into two daughter agents which are initialised in the first phase. This can be summarised by the following chain of reactions:
In general, the resulting probability distribution of the total cell cycle is a hypoexponential distribution. Note that this general implementation involves independent parameters (), one for each stage transition and, if is large, then this may lead to issues of parameter identifiability. In order to reduce the number of free parameters while maintaining the advantage of using a multistage representation, Yates et al. [2017] focus on two cases: the case in which all the transition rates are identical, for , and the case in which all the transition rates are identical apart from one, for . The corresponding distribution for the total cell cycle in these two cases are the Erlang distribution and exponentially modified Erlang distribution, respectively. The authors show these distributions to be both biologically plausible and computational feasible.
In order to investigate how the multistage representation of the cellcycle affects the total population growth, Yates et al. [2017] implemented two models (one spatial and one nonspatial). Firstly, they modified the model of Turner et al. [2009] in which the spatial position of the cells is considered unimportant. The authors investigate the stage distribution of agents at the steady state for large times. By writing down a system of ODEs describing the average proportion of agents at each stage , , for , Yates et al. [2017] show that the number of agents in each stage is not proportional to the average length of that stage. In particular, for identical transition rates, as the total number of stages becomes large, the proportion of agents at the first stage approaches twice that of agents at the last stage.
The authors also modify the spatially extended ABM of Baker et al. [2010] to include their multistage proliferation scheme. The agents are initialised uniformly on a lattice with periodic boundary conditions and a standard volume exclusion property. Agents attempt movements with rate . To facilitate the comparison with the traditional model with exponentially distributed waiting time of rate , the authors consider the case of identical transition rates for . Hence the average waiting time required for progress through all the stages is independent of the number of stages. In particular, the two models, with and without multistage scheme, have the same average attempted proliferation waiting time.
The simulations of the ABM highlight that, for lowdensity colonies, the multistage scheme leads to a slower population growth compared to the model with exponentially distributed proliferation. However, under the realistic assumption that agents that fail to proliferate due to crowding remain at the last stage (as opposed to returning back to the first stage of the cycle), a clear proliferating rim of (grey) cells can be seen with the bulk of cells being kept at stage (see Figure 10). Agents located in this rim, reattempt division after aborted events more quickly than in the singlestage cell cycle model since they only have to wait a time on average. Therefore, the effective average interdivision time for cells with a multistage cell cycle at the proliferating rim of a cluster decreases in comparison with cells to a singlestage cell cycle.
——————————————–
3.2 Cell interactions
Cells can undergo a great variety of interactions in response to their environment and more specifically, their neighbouring cells [Cai et al., 2006, Dworkin and Kaiser, 1985, Trepat et al., 2009, Tambe et al., 2011]. Typically, these can be divided in two broad categories: indirect and direct [Othmer and Stevens, 1997].
When cells have the ability to detect the presence of an external signal and change their behaviour according to its concentration, we call these indirect interactions. This can affect the movement speed or the turning rate (kinesis) [Cai et al., 2006], it may induce a directional bias (taxis) [Painter and Hillen, 2002] or it may involve a combination of these [Erban and Othmer, 2004]. In general, the external signal can depend on the cells themselves. For example, it can be directly produced by a moving cell as it moves, forming a trail behind it, as it has been shown for Myxobacteria [Dworkin and Kaiser, 1985]. Notice that this type of interaction does not involve direct cellcell contact, communication is mediated only through the external signal. Direct interactions are based on the ability of cells to detect and interact with surrounding cells. These include contact interactions such as volume exclusion [Abercrombie, 1979] and adhesionrepulsion [Trepat et al., 2009, Tambe et al., 2011].
The introduction of these interactions in models of cell migration is crucial, since they often constitute the driving forces that generate spatial structure at the population scale. In this Section we focus our attention on a set of important cell interactions: the indirect response to a chemical signal, direct adhesionrepulsion forces between cells and the ability of cells to push and pull each other.
3.2.1 Chemotaxis
Othmer and Stevens [1997] derive a general class of PDEs from a set of ABMs incorporating a simple signal following mechanism. The aim of their paper is to determine whether stable bacterial aggregation can arise as result of a simple chemotactic response.
The main ABM considered in the paper consists of a system of noninteracting agents which move on a onedimensional lattice with unit step size, according to a continuoustime, nearestneighbour random walk. As time evolves, agents produce an attracting signal, , which sits on an embedded lattice of half the step size.
The authors distinguish three types of model depending on the information used to compute the transition rates. These are the local model, if the transitions depend only on the signal concentration at the current site, ; the barrier model, in which cells only sense the value of the signal at the shortrange, nearestneighbour sites, ; and the gradientbased model, in which the transition rates depend on the longrange, nearestneighbour differences, . In the barrier and gradientbased cases an additional distinction is made between normalised and unnormalised transition rates. For unnormalised rates, larger values of lead to a faster total movement rate, whereas, for normalised rates, the average waiting time at each site does not depend on the signal concentration. In all cases, by using a limiting argument, the authors derived the corresponding continuous diffusive approximations for the average agent density.
Othmer and Stevens carry out a stability analysis of the solutions of the continuous barrier model with normalised transition rates. The analysis is repeated with three different models for the evolution of the signal: linear growth, exponential growth and saturating growth. The results show that linear growth can only lead to uniform agent density profiles, even when starting from a single peak in the initial distribution of agents (collapse). On the other hand, when the signal grows exponentially, the system tends to converge to a single peaked distribution (blowup). Finally, when a saturation in the production of is assumed, stable aggregation can occur as a result of the interplay between the production of the signal and the shortrange chemotactic response. The formation of such aggregates strongly depends on the decay rate of the signal and on the initial condition. In particular, a large decay rate or a small initial value of can lead to aggregation or blowup, whereas, in the absence of signal decay, the agent density will eventually collapse to uniformity.
3.2.2 Adhesionrepulsion
Anguige and Schmeiser [2009] investigate the emergence of aggregation through the mechanisms of cellcell adhesion and volume exclusion. They predict, under their model, that the spontaneous aggregation of cells is not possible if the initial cell density throughout the domain is too low, regardless of the intensity of the adhesion force.
To describe cell migration, Anguige and Schmeiser [2009] consider a discretespace, continuoustime random walk model on the unit interval. Multiple agents can occupy the same compartment up to a carrying capacity, , and they move at random with given rates to one of the two nearest neighbour compartments. In order to include volume exclusion in the model, the transition rates in both directions are decreased linearly with the density of the target site. If a site is fully occupied, no transitions into that site can occur. To incorporate adhesive forces, the rate of moving in a particular direction is decreased linearly with the density of the adjacent site in the opposite direction, in proportion to an adhesion parameter, .
By Taylor expanding and neglecting terms of third (or higher) order it is possible to derive the diffusive limit of the discrete model as the number of lattice sites goes to infinity. This macroscopic model is a nonlinear diffusion equation as equation (30) with
for which the homogenous density profile is the only steady state. Anguige and Schmeiser [2009] find that there is a critical value of the adhesion parameter . For the lowadhesion regime (), pattern formation is not possible in either the discrete or continuum descriptions.
When the adhesion force is large (), the results of the discrete model show complex behaviours such as pattern formation and spatial oscillations in density. Patterns in the discrete model are found to be transient and metastable. All clusters which arise eventually coalesce to single stable aggregate or to two aggregates separated by a single trough. For the macroscopic PDE representation, the authors identify an interval of unstable values of density for which the diffusion coefficient of the nonlinear PDE takes negative values making the continuum model illposed.
To overcome the issue of the illposedness in the PDE and to obtain a reasonable continuum description, Anguige and Schmeiser [2009] propose to include more terms from the Taylor expansion of their discrete model. These higher order corrections result in a fourthorder diffusion equation reminiscent of the CahnHilliard equation [Sun and Ward, 2000]. The presence of a viscosity terms in the revised equation allows the authors to prove it has a wellposed initial value problem.
Thompson et al. [2012] continue the work of Anguige and Schmeiser [2009] on the modelling of cellcell adhesion at multiple scales. The authors consider a modification of the ABM of Anguige and Schmeiser [2009] in order to study the interactions with a second species of cell. In this new model, agents are of two types, and , and the there are four different coefficients ( and ) governing the intra and interspecies adhesion forces, respectively. Thompson et al. [2012] show that the model is capable of reproducing three configurations (complete sorting, engulfment and mixing), which have been observed previously both in experiments and continuous models [Armstrong et al., 2006]. In Figure 11 we report an example of these three configurations depending on the choice of the adhesion parameters. Specifically, if interspecies adhesion is small and intraspecies adhesion is large the system reaches a structured configuration (cell sorting) in which agents of the same species tend to cluster together (see left column of Figure 11). If interspecies adhesion is small and one of the intraspecies adhesions is larger than the other, then engulfment of the species with the larger adhesion occurs (see middle column of Figure 11). Finally, when interspecies adhesion is large and the intraspecies adhesion is small, then mixing occurs, i.e. agents of the two species are uniformly distributed across the domain (see right column of Figure 11).
Recently, Binny et al. [2015] have studied a onedimensional offlattice ABM (later extended to twodimensions [Binny et al., 2016]) in which both the rate and the direction of the movement of each agent are influenced by the configuration of the neighbouring agents. In their model an agent, , moves with a rate, , which is computed by
where is an intrinsic motility rate, independent of the other agents’ positions, and is a kernel function weighting the strength of interaction between agents positioned a distance from each other. Binny et al. [2015] consider the case in which the kernel is Gaussian:
(33) 
where and determine the intensity and the range of the interaction. Notice that if , agents tend to move more often when they are surrounded by other close agents. Conversely, for , the presence of close neighbours inhibits the motility. When a movement even takes place, the moving agent performs a jump of random length and direction. The length of jump is drawn at random from a Laplace distribution (i.e. short steps are more likely to be taken than long jumps), independent of the other agents’ positions. For the direction of the movement, the authors incorporate a directional bias, , such that the presence of neighbouring agents can affect the final direction. They defined the directional bias for an agent, , as
where is a Gaussian kernel function of the same form of equation (33) with intensity and rage given by and , respectively. If , agents are biased to move away from highly concentrated regions, whereas if , agents tend to move towards one another.
To study the resulting spatial structure from a deterministic prospective, Binny et al. [2015] develop a populationlevel model in terms of the first two spatial moments of the ABM. By using the Kirkwood superposition approximation (see Section 2.3 for a discussion about alternative approaches to the meanfield approximation), the authors derive a closed system of two ordinary differential equations (ODEs) describing the evolution of the average density of single agents, and of pairs of agents separated by a distance , .
In order to test the accuracy of their deterministic representation against the numerically simulated individualbased model, the authors compare the pair correlation functions (PCFs) [Illian et al., 2008] of the ABM with the PCF predicted by the spatial moment model. The pair correlation function can be expressed in terms of the first two spatial moments as
(34) 
Notice that corresponds to the case of absence of spatial structure, whereas denotes evidence of spatial correlation and implies anticorrelation, at length .
Their findings show that the moment model provides a good approximation of the spatial structure predicted by the ABM for a wide range of parameter choices. The moment model underestimates or overestimates the spatial structure only for large values of intensity and . Binny et al. [2015] suggests that this is due to higher order interactions which are not incorporated in the deterministic representation.
Binny et al. [2016] continue the study of agentagent interactions by extending the previous model to two dimensions and including agent proliferation and death. In the new model, agents are allowed to proliferate with a rate which is dependent on their neighbourhood and defined similarly to of equation (33). The rate of agent death is assumed independent of the location of other agents. Following Binny et al. [2015], the authors derive a deterministic model for the first two spatial moments of the ABM. After comparing the performance of four different types of moment closure, the authors select an asymmetric power2 closure as the most appropriate [Murrell et al., 2004, Raghib et al., 2011] and they adopt it throughout the paper. The results show that the good agreement between the stochastic and deterministic models, found in [Binny et al., 2015] is preserved when proliferation and death are incorporated in the model. Unsurprisingly, the presence of nonuniform spatial structure leads to both population growth and the average cell density differing significantly from the predictions of a meanfield model that ignores agentagent interactions.
Finally, Binny et al. [2015, 2016] investigate the role of neighbourdependent movement and proliferation in the formation of spatial structure. In general, movement tends to break up spatial structure because new agents can move out of clusters generated by shortrange proliferation [Baker and Simpson, 2010]. When movement is biased away from neighbouring agent, , this effect appears more clearly, since agents undergo directed movement out of clusters. The same dispersal of spatial structure is found for the neighbourdependentinhibition of the motility rate, but it occurs in a less efficient way. The authors also highlight that the crowdinginduced inhibition of proliferation is capable of counteracting the formation of clusters due to shortrange proliferation.
3.2.3 Pushing and pulling
Yates et al. [2015] study how cellcell pushing affects cells’ dispersal and proliferation by considering on and offlattice ABMs. The authors begin by considering a simple ABM incorporating exclusion on a twodimensional lattice, as in Simpson et al. [2007]. They modify the basic model in order to incorporate the ability of cells to push their neighbours [Ewald et al., 2008, Schmidt and Friedl, 2010]. Pushing is implemented according to four different mechanisms with variable degrees of complexity and realism.
Consider an agent at position which attempts a movement into an occupied site to its right. In the most basic pushing mechanism, the authors let the moving agent push the agent at position into the site , with a given probability, , with the pushing agent taking the pushed agents’ original position, . The pushing is effective only if the target site of the pushed agent, , is empty. The entire movement is aborted otherwise. For the next two mechanisms, the assumption on the location of the target site of the pushed agent is relaxed. For one of these mechanisms the pushed agent is allowed to move into each of its empty neighbouring sites with equal probability. If all neighbouring sites of the pushed agent are occupied, then the whole event is aborted. In the other mechanism, the target site is chosen uniformly at random from amongst the pushed agent’s three nearest neighbours, but if the chosen site is found to be occupied then the entire event is aborted. In the final case considered, the moving agent can attempt to push up to other agents in a straight line in a chosen direction.
In all the cases considered, the authors derive diffusive PDE approximations for the total agent density. For the basic pushing mechanism, for example, the resulting nonlinear PDE takes the form of equation (30) with density dependent coefficient given by
(35) 
In most scenarios, the comparison between the ABM and the corresponding PDE is good. However, the more complicated pushing mechanisms, such as the linear pushing of multiple agents, are found to introduce strong spatial correlation in the occupancies of adjacent sites. Since the derivation of the continuum models is based on a meanfield approximation, the presence of such spatial correlation leads to a divergence between the ABM and its deterministic counterpart.
Yates et al. [2015] also investigate the effect of introducing pushing in an offlattice setup by modifying the onedimensional model of Dyson et al. [2012] (see Section 2.1). In this model, a pushing event is attempted with probability when a movement would lead to an overlap of agents, which are represented as interval of length . The moving agent can displace the pushed adjacent agent far enough for it to complete its movement, but only if this will not produce an overlap with a third agent. In this case the entire movement would be aborted. By using a similar approach to that of Dyson et al. [2012], a nonlinear diffusive PDE of the form (30) is derived with
(36) 
and the numerical solution of the PDE shows good agreement with the averaged density of the ABM.
In all the cases considered, both on and offlattice, Yates et al. [2015] find that the introduction of pushing behaviour leads to densitydependent diffusion coefficients in the macroscopic descriptions. For the onlattice models, pushing leads to a faster diffusion, due to the larger number of possible movements for each agent in comparison to the simple exclusion model. However, in the offlattice scenario, increasing the tendency to push leads to a decreasing diffusivity. To explain this behaviour, Yates et al. [2015] underline that when a pushing event occurs in the offlattice model, the moving agent and the pushed agent end up being adjacent, even if they were not in contact prior the movement. This aggregation effect leads to a slower dispersal of agents corresponding to a lower diffusion in comparison to the nonpushing case.
There is much evidence in the literature to suggest that cells can pull each other while they are undergoing collective migration [Bianco et al., 2007, Ghysen and DamblyChaudière, 2007]. Chappelle and Yates [2018] explore the effect of allowing cellcell pulling in both on and offlattice models of cell migration.
Firstly, the authors specify a simple onlattice excluding ABM in which agents are capable of pulling each other. For example, if an agent is moving to the right and another agent is occupying the neighbouring site to its left, the moving agent will pull the neighbouring agent rightwards, with probability .
By writing down the occupancy master equation of the model and taking appropriate limits, the authors obtain a nonlinear diffusion equation for the average agent density, which takes the form of equation (30) with
Pulling is found produce a quadratic densitydependence in the diffusion coefficient of the macroscopic PDE. A comparison with the corresponding model for agents which are able to push each other (see equation (35)) highlights that the density dependence of diffusivity due to pushing is larger than in the pulling case.
In order to study the effect of pulling in offlattice models, Chappelle and Yates [2018] use a similar approach to that of Dyson et al. [2012] (see Section 2.1). The onedimensional model in Dyson et al. [2012] is extended by introducing a pulling distance proportional to the moving distance, , such that, if an agent is chosen to move in a given direction and there is another agent whose centre is within distance in the opposite direction, both agents move in the given direction. The corresponding PDE is of the form of equation (30) with
which predicts that the effect of pulling is to decrease the effective diffusivity in the offlattice model, leading to a slower dispersion at the macroscopic level. The apparent contradiction with the onlattice counterpart is explained in the paper by using a similar argument to Dyson et al. [2012].
3.3 Growing domains
Considering cell motility on a purely stationary domain may often be unrealistic for biological processes: cell growth, division and movement itself will cause the size of the tissue to change dynamically in many biologically plausible situations [Rogulja and Irvine, 2005, Wolpert et al., 2015]. Surprisingly, it is only in the last decade that attention has been given to investigating how this phenomenon affects models at multiple scales. Baker et al. [2010] represents a pioneering example of the incorporation of domain growth into models of cell migration.
For the motility scheme, Baker et al. [2010] follow the approach of Othmer and Stevens [1997]. At the individualcelllevel they consider a continuoustime, discretespace ABM. Agents are initialised on a onedimensional lattice with zeroflux boundary conditions. They have the ability to sense the concentration of a signal profile, , at the their current lattice site and at their two nearest neighbour sites. Agents can jump to immediately neighbouring lattice sites. These jumps are regulated by two transition rates that are linear combinations of the signalling molecule concentration at the agent’s current and immediately sites. The authors consider four types of transition rate comprising different linear combinations of these concentrations (local, nonlocal, average and difference). If no signal profile exists then agents can either diffuse randomly (jumping with constant rates independent of their position) or implement densitydependent transition rates. For each case, using a master equation, the authors derive populationlevel descriptions, which comprise an advectiondiffusion PDE for the average agent density, , as a function of position, , and time, .
The inclusion of domain growth in the model is carried out in two phases. Firstly, Baker et al. [2010] consider the case of exponential growth. On the individuallevel this implies that each lattice site divides at a constant rate per unit time. When the site divides, a new daughter site is added adjacent to the parent site. The agents in the parent site are divided between these two sites according to a symmetric distribution and the sites to the right of the daughter site (and their contents) are all shifted one site’s width to the right. A PDE description of cell density is derived both from a phenomenological, continuum perspective, using a conservation of matter argument [Crampin et al., 1999], and from the master equation of the ABM. In order to derive the corresponding PDE from the master equation, the authors use a moment closure approximation on the number of agents in each lattice site. This approximation is based on the assumption that movement occurs on a faster timescale than domain growth.
Secondly, the authors consider a more general type of domain growth which is densitydependent. As before, upon being chosen to undergo a growth event, a lattice site is divided into two daughter sites, but now with a densitydependent rate, . A continuous approximation is derived using a further moment closure assumption  that the mean of a nonlinear function of the agent density can be expressed as the same nonlinear function of the mean agent density. The authors prove that domain growth is linear for the case of linear densitydependence on domain growth, while it must be evaluated for numerically in the general case.
The paper shows comparisons between the simulations of the individualbased models and the numerical solution of the PDEs. In addition, the predicted macroscopic domain length is compared with the average value from simulations. Results are displayed for the cases of the nongrowing domain, constant growth, linearly densitydependent growth and quadratically densitydependent growth. In each case, the results show a good agreement between the discrete model density and the continuous approximation. However, the predicted value of the domain length in the continuum approximation is underestimated in comparison to the true value from the stochastic model in the case of the nonlinear growth model.
The last part of the paper incorporates a signalling profile which agents can sense and respond to. Specifically, an exponentially decreasing concentration of signalling molecule, , is placed in the domain and agents are assumed to interact with it, both via the local and nonlocal schemes. Results show good agreement between the simulation and the deterministic prediction for both linear and logistic domain growth. As a last application, Baker et al. [2010] consider linearly densitydependent growth with both local and nonlocal sensing again demonstrating good agreement between the two modelling regimes.
Yates et al. [2012] generalise the existing ABMs of Baker et al. [2010] to a nonuniform lattice allowing the process of cell division on the underlying domain to be modelled in a more realistic manner. The domain is allowed to extend in a manner which can be made arbitrarily close to continuous, as opposed to the discrete increments by which the domain is extended in Baker et al. [2010].
The authors consider two distinct ways to discretise a onedimensional domain in a nonuniform way: the Voronoi partition method and the intervalcentred method. In the first case, the agent positions, for , are chosen first and the edges of the intervals are the bisectors of these points. In the second case, the edges of the intervals are specified first and the agents lie at the centre of these intervals. As in Baker et al. [2010], no exclusion property is implemented and multiple agents can occupy the same lattice site.
In the first part of the paper the domain is fixed and agents move following a positionjump process on the irregular lattice. The transition rates of an agent in interval are computed from the mean first passage times that a Brownian particle, initialised at , takes to hit one of the neighbouring particle positions ( or ). Using the master equation for agent densities, Yates et al. [2012] derive a macroscale description of the model in the diffusive limit. For the Voronoi partition, the density evolves according to the diffusion equation and the comparison between the simulations and the PDE shows a good agreement. In the case of the intervalcentred partition, the agreement between the individuallevel behaviour and the populationlevel diffusion equation is poor since the transition rates are not inherently linked to the size of the intervals.
The remainder of the paper focuses on the introduction of domain growth to the models of agent migration using the Voronoi partition. Yates et al. [2012] start by defining a deterministic scheme for domain growth. Every time an agent moves, every interval grows an amount proportional to its length and proportional to the time step between movement events. When an interval reaches a threshold length, it splits into two daughter intervals. The boundaries of such daughter intervals are chosen in order to preserve the Voronoi property of the domain partition. Each agent is redistributed to a new interval with probability proportional to the overlap between the new intervals and the old intervals. The authors also consider an alternative growth mechanism in which, with a give rate, an interval is chosen at random to grow with probability proportionally to its length. When a growth event occurs, both the selected interval and an adjacent interval grow in order to preserve the Voronoi property. For both growth schemes, a populationlevel description can be obtained through the master equation for the densities. The resulting PDEs are of the same form as in [Baker et al., 2010].
Finally, a series of comparisons between simulations and PDEs confirm the good agreement between the microscopic and macroscopic models for the Voronoi domain partition, but a poorer agreement for the intervalcentred domain partition, as expected.
In addition to their investigations into cellcell adhesion and volume exclusion on a static domain (see Section 3.2.2) Thompson et al. [2012] study the interplay of the same properties on a growing domain using a modification of the ABMs of Baker et al. [2010] and Yates et al. [2012]. The authors incorporate a partially excluding property into the model of domain growth of Baker et al. [2010] by allowing the carrying capacity of each compartment to be proportional to its size (see Figure 4 for a schematic illustration). When a growth event occurs, a compartment is chosen at random and its carrying capacity is increased by unity. Commensurately, its size it is also increased. When a compartment reaches a predefined threshold, it is split into two compartments each with carrying capacity set to half of the previous value. To compensate for the unequal compartment sizes generated by growth, the transition rates are amended as in [Yates et al., 2012]. A populationlevel description of the model is derived, but it is only valid for short timescales, specifically until the first splitting event occurs.
The results demonstrate that domain growth decreases the chance of cell clustering occurring by cellcell interactions. For high values of adhesion and with small growth rates, clusters can still appear over short timesales. However, as the growth rate increases, the probability of clusters appearing decreases and all initial cell clusters are eventually destroyed.
Hywood et al. [2013] suggest a modified version of the model of Baker et al. [2010] to represent tissue growth. The authors initialise a onedimensional lattice with a set of contiguous nonoverlapping agents (tracers agents). These agents are inactive (i.e they do not perform any jumping movement between sites) and their role is only to mark the position of the site in which they are located as the underlying lattice grows. The authors focus on the case of exponential growth with constant rate, , which is implemented similarly to Baker et al. [2010]. If a marked site splits into two daughter sites, the corresponding tracer agent moves to the right daughter site (see Figure 12 for an illustration). The usage of tracer agents follows from the previous work of Binder and Landman [2009] in which, contrastingly, domain growth is implemented in a deterministic manner.
By coupling the position of the tracers to a system of noninteracting random walkers and writing down the corresponding occupancy master equation for the tracers, Hywood et al. [2013] derive a formula for the infinitesimal mean and variance of the underlying stochastic process, and , respectively. These expressions can then be used as coefficients for an advectiondiffusion PDE (FokkerPlanck equation (FPE)) describing the spatiotemporal evolution of the average occupancy of the tracer agents, :
(37) 
Yates [2014] extends the work of Hywood et al. [2013] to more general scenarios in which the rate at which sites divide is timedependent, . This extension allows the incorporation of a variety of biologically realistic mechanisms of domain growth. By using similar steps, Yates [2014] obtains a series of PDEs for the average occupancy of the tracer agents of the same form as equation (37). In particular, the author provides the expressions of and for the case of linear growth, , generalised logistic growth, , and Gompertzian growth, .
Finally, Yates [2014] studies the implications of implementing site death in the model which can be done by removing sites with a given rate, . When this occurs, all the site to the right of the removed site are moved to the left in order to fill the vacant space created by the site’s removal. If a tracer agent occupied the removed site, it remains in its position. Notice that this new mechanism can lead to multiple tracer agents occupying the same site. See Figure 12 for an illustration. The comparison between the ABM simulations and the corresponding PDE preserves a good agreement when site death is introduced, even for cases of net domain shrinkage, i.e .
The results of Yates [2014] highlight the danger of neglecting death events in the model, even when the net growth rate is positive, . In other words, although the mean growth rate can be estimated correctly by using purely growing dynamics (with growth rate given by ), the second and higher moments of the process will be incorrect. For example, the higher value of variance leads to faster diffusivity of the agents when site death is explicitly incorporated in the model.
3.4 Persistence of motion
A fundamental assumption of modelling cell movement with a positionjump process (such as those implemented in Section 2), is that the cells’ positions are subject to a series of Markovian jumps in space and, therefore, the direction of motion between sequential jumps is totally uncorrelated. In reality, experimental observations suggest that many types of cell show a tendency to preserve their direction of motion for some time before reorienting [Berg and Brown, 1972, Hall, 1977, Gail and Boone, 1970, Wright et al., 2008], even when the movement remains unbiased in the long term. This behaviour is known as persistence of motion or of direction [Patlak, 1953]. The standard framework to model persistence of motion in an ABM based is to employ a velocityjump process [Othmer et al., 1988, Othmer and Hillen, 2000, 2002, Codling et al., 2008, Campos et al., 2010] in which the variable which is performing markovian jumps is the velocity rather than the position of the agents [Othmer et al., 1988].
A simple example of an velocityjump process ABM in one dimension was analysed by Goldstein [1951] and later by Kac [1974] and Othmer et al. [1988]. The model can be formulated as follows. Consider an ABM in which agents have an assigned direction of motion, either right or left. As time evolves they can either move, with rate , or change their direction, with rate . The two events are assumed to occur independently. When a movement event takes place, the moving agent performs a jump of distance, , in its assigned direction. We can interpret as the two possible values of the velocity of an agent, depending whether it is moving in the right or left direction, respectively. Notice that this formulation is valid for both on and offlattice models. However, in an onlattice framework, the distance, , has to be an integer multiple of the lattice step .
Let and denote the average occupancy of agents in position at time with associated direction to the right and left, respectively. For a system of noninteracting agents, the corresponding occupancy master equations can be written as
(38a)  
(38b) 
By Taylor expanding the two terms and about to first order and taking the limit with kept constant, one obtains the system of advective PDEs given by
(39a)  
(39b) 
where . By adding equations (39) and differentiating with respect to we can write
(40) 
Similarly, by subtracting equation (39a) from equation (39b) and differentiating with respect to , we obtain
(41) 
Finally, substituting (40) into equation (41) and recalling equations (39) we can write