Random walk models in biology

2.1 The simple isotropic random walk

Consider a walker moving on an infinite one-dimensional uniform lattice (i.e. a line split into discrete points). Suppose the walker starts at the origin (x=0) and then moves a short distance δ either left or right in a short time τ. The motion is assumed to be completely random, so the probabilities of moving both left and right are 1/2. After one time step, the walker can either be at a distance δ to the left or right of the origin, with probability 1/2 each. After the next time step, the walker will either be at a distance 2δ to the left or right of the origin (with probability 1/4 each) or will have returned to the origin (with probability 1/2). Note that, after an even (odd) number of steps, the walker can only be at an even (odd) distance away from the origin. Continuing in this way, the probability that a walker will be at a distance mδ to the right of the origin after n time steps (where m and n are even) is given by

(2.1)

(2.2)

Useful time-dependent statistics of this process are the mean location E(X_t) and the mean squared displacement (MSD) , defined as

(2.3)

2.2 A BRW with waiting times

It is also a standard procedure (Lin & Segel 1974; Okubo & Levin 2001) to derive the governing equation for a SRW using a difference equation, and subsequently solve to find solutions of the form (2.2). This method also allows more complex processes to be modelled as part of the random walk and we illustrate this in the following example by including both a preferred direction (or bias) and a possible waiting time between movement steps.

Consider a walker moving on a one-dimensional lattice, where, at each time step τ, the walker moves a distance δ to the left or right with probabilities l and r, respectively, or stays in the same location (‘waits’), with probability 1−l−r (the isotropic random walk in §2.1 has r=l=1/2 and there is no waiting time). Now, if the walker is at location x at time t+τ, then there are three possibilities for its location at time t: (i) it was at x−δ and then moved to the right, (ii) it was at x+δ and then moved to the left, and (iii) it was at x and did not move at all. Thus we have

(2.4)

(2.5)

(2.6)

(2.7)

It is worth noting that, owing to the waiting time between movements, the value of the diffusion constant D is smaller (k=l+r<1 in this case) than that in §2.1 and hence, as expected, diffusion is less rapid. Note also that, due to the way the limits are taken in (2.6), the term ϵ=r−l must tend to zero in order to obtain a finite drift rate u.

The solution of (2.7) subjected to the initial condition p(x, 0)=δ_d(x) (the Dirac delta function), which means that the walker is at x=0 at time t=0, is (Montroll & Shlesinger 1984; Grimmett & Stirzaker 2001)

(2.8)

The mean location and the MSD can be easily calculated by substituting (2.8) into (2.3). Alternatively, it is possible to calculate the moments directly from the governing PDE (2.7). This is a technique that is particularly useful in more complex models when the governing differential equation is known, but the solution for p(x, t) may be difficult or impossible to find (see §3.5). First, we multiply equation (2.7) either by x (to find mean location) or by x² (to find MSD) and integrate by parts. Then, using the definitions in (2.3), the fact that

(2.9)

(2.10)

2.3 The BRW in higher dimensions

A similar derivation to that in §2.2 can be completed using an N-dimensional lattice to give the standard drift–diffusion equation

(2.11)

(2.12)

(2.13)

• Download figure
• Open in new tab
• Download PowerPoint

(2.14)

It is interesting to note that, if two-dimensional movement is not restricted to a lattice (box 1), then diffusion in the direction of the bias is lower than that in the perpendicular direction (see ch. 2 of Codling 2003; Coscoy et al. 2007). This result is not always properly considered in spatial population dynamics models, where diffusion is usually assumed to be the same in all directions even when there is an average drift in a particular direction.

Box 1. Stepping off the lattice.

The simple random walk models discussed in §2 are mostly restricted to N-dimensional lattices so that there are only a finite number of choices of direction at each time step. A more realistic model allows for a continuous choice of direction. In two dimensions, this means that the walker is allowed to move in any direction θ on the unit circle. Linear statistical measures cannot be used because any angular value θ is only defined modulo 2π, so θ and θ+2π correspond to the same direction. Instead, circular distributions can be used to draw random angles for either the direction of movement (in an SRW or BRW) or turning angle (in a CRW or BCRW) at each step of a two-dimensional random walk.Useful moments of a circular distribution (with PDF f(θ)) are the mean cosine c and the mean sine s, defined as

The mean angle is given by θ₀=arctan (s/c) if c>0 and by θ₀=arctan (s/c)+π if c<0, and the mean vector length r as r=(s²+c²)^1/2, which gives a measure of how peaked the distribution is about the mean angle (Mardia & Jupp 1999).Commonly used circular distributions are (Batschelet 1981; Mardia & Jupp 1999) as follows:

1. The von Mises distribution

   
   where I₀ denotes the modified Bessel function of the first kind and order 0, defined by . The mean vector length for this distribution is r=I₁(κ)/I₀(κ).

2. The wrapped normal distribution

   
   The mean vector length for this distribution is , where sigma is the standard deviation of the linear normal distribution before wrapping.

3. The wrapped Cauchy distribution

   

The wrapped normal distribution is the standard normal distribution ‘wrapped’ around a unit circle, and is hence easy to interpret and simulate. The von Mises distribution is easier to deal with analytically, however, and is similar to the wrapped normal distribution with the same mean cosine, so is often used instead of the wrapped normal. The wrapped Cauchy distribution is more peaked, with fatter tails than the von Mises and wrapped normal (Mardia & Jupp 1999).In three-dimensional space, a second angular direction ϕ must also be specified, and spherical distributions defined (Fisher et al. 1987). Spherical distributions and data are difficult to deal with (which may be the reason why many movement studies are restricted to two dimensions), although it is possible to use a vector-based approach (e.g. Uttieri et al. 2007).

2.4 The Fokker–Planck equation

It is also possible to extend the SRW in two (or more) dimensions to include movement probabilities that are spatially dependent. This results in the Fokker–Planck equation for p(x, y, t) (Okubo & Levin 2001). Suppose that an individual moves on a two-dimensional lattice. At each time step τ, an individual can move a distance δ either up, down, left or right with probabilities dependent on location, given by u(x, y), d(x, y), l(x, y) and r(x, y), respectively (with u+d+l+r≤1), or remain at the same location with probability 1−u(x, y)−l(x, y)−d(x, y)−r(x, y). We now use a difference equation as in §2.2, expand as a Taylor series and define the following parameters:

(2.15)

2.5 General diffusive properties and model limitations

Closer inspection of the diffusion equation solutions (2.2) and (2.8) shows that the following important property holds: p(x, t)>0 for all locations x and all positive times t. Hence, for any positive time (no matter how small), there is a positive probability of being at any finite location.

In practice, the probability of being at an exceptionally large distance away from the origin after an infinitesimal time step is extremely small. However, the fact that this property of ‘infinite propagation speed’ exists suggests that the limiting process involved in deriving the governing equation of the SRW model has some limitations. As shown in §2.2, to derive the governing equation of a simple BRW, we take limits assuming that δ²/τ=2D is constant as δ, τ→0 (where δ is the distance moved at each time step τ). Clearly, if this property holds then, in the same limit, we get δ/τ→∞, i.e. the effective instantaneous speed is infinite (Okubo & Levin 2001). The solution of the diffusion equation should hence be considered only as an asymptotic approximation, valid for large times (t≫τ), of equations that more accurately describe the correlations in movement and finite speeds that are present when considering movement at shorter time scales.

A more detailed discussion of the infinite propagation speed in diffusive processes can be found in Okubo & Levin (2001), while in §2.7 we present a simple model based on a velocity jump process (rather than a position jump process), which avoids this problem.

2.6 Random walks with a barrier

To model movement in a confined domain, one can introduce a repelling or reflecting barrier into the random walk: a walker reaching the barrier will automatically turn around and move away in the opposite direction. Similarly, to model movement where walkers leave the system upon reaching a given point, one can introduce an absorbing barrier. Models such as these are not only appropriate to model movement in space, but they are also suited to modelling development and growth where critical life stages are reached (Pitchford & Brindley 2001; Pitchford et al. 2005; Mullowney & James 2007) and to ‘integrate and fire’ models of nerve responses (Iyengar 2000). These are examples of first passage time (or first hitting time) problems, in which the distribution of the time taken to reach an absorbing barrier is of primary interest (Condamin et al. 2007).

The simple example of an absorbing barrier given below is adapted from an example by Grimmett & Stirzaker (2001). Suppose we have a one-dimensional random walk process that satisfies the drift–diffusion equation (2.7) for x>0. Suppose the walker starts at location x=x₀>0, and there is an absorbing barrier at x=0 such that, if the walker reaches the point x=0, it is removed from the system. The appropriate boundary and initial conditions in this case are

(2.16)

(2.17)

(2.18)

(2.19)

(2.20)

A similar analysis is possible for a reflecting barrier (Grimmett & Stirzaker 2001) and the SRW model can be extended or analysed in many other ways. Montroll & Shlesinger (1984) give a good review of the general theory of random walks and discuss a wide variety of ideas and problems.

2.7 CRWs and the telegraph equation

A CRW takes into account short-term correlations in the direction of movement. In most cases, this means that the walker is more likely to move in the same or a similar direction to its previous movement direction. This tendency to continue in the same direction is known as persistence (Patlak 1953; see for example figure 1d). By explicitly including persistence and a fixed speed of movement in the random walk process, the problem of infinite propagation speed discussed in §2.5 is avoided (see also the discussion in Turchin 1998). The location at each step of the random walk is no longer a Markov process (as it depends on the sequence of previous locations). Hence the usual framework for describing a CRW is a velocity jump process, in which the variable following a Markov process is the walker's velocity rather than the location (Othmer et al. 1988).

Consider a population of individuals moving either left or right along an infinite line at a constant speed v. Denote the density of right- and left-moving individuals at location x and time t by α(x, t) and β(x, t), respectively. The total population density is p(x, t)=α(x, t)+β(x, t). At each time step τ, each individual either changes direction and moves a distance δ in this new direction (with probability r=λτ), or moves a distance δ in the previous direction (with probability q=1−λτ). Hence, turning events occur as a Poisson process with rate λ. If we take a forward time step, then the number density of individuals at location x moving right and left, respectively, is given by

(2.21)

(2.22)

(2.23)

(2.24)

(2.25)

(2.26)

(2.27)

Note that, although (2.27) describes a correlated movement process, the random walk is globally unbiased in the sense that there is no overall preferred direction, simply a tendency for individuals to persist in their present direction of motion (a localized bias). Although the original process uses a fixed time step τ, the mean time between turning events, , is different since the telegraph process does not have a turning event at every time step.

Equation (2.27) can be solved given the initial conditions p(x, 0) and (∂p/∂t)(x, 0), but the full solution is quite complex (see Morse & Feshbach 1953, for details). Since there is a fixed speed v, the solution does not imply infinite propagation speeds (as found with the solutions (2.2) and (2.8) of the diffusion equation).

It is possible to use a similar method to that used in §2.2 to derive equations for the moments of the solution: E(X_t)=0 and

(2.28)

Box 2. Anomalous diffusion.

As discussed in §§2 and 3, the MSD , defined in (2.13), for a typical diffusive random walk is linearly related to time (or number of steps) as in (2.14). However, there are situations where dispersal may not be diffusive, and the MSD is not linear in time but instead has some other power-law relationship. Such situations are known as anomalous diffusion (e.g. Weeks et al. 1996). The key parameter is the value of μ that is observed in the relation . There are five possible situations to consider.

1. μ=0. This corresponds to a stationary process with no movement over the period of observation.

2. 0<μ<1. This situation is known as sub-diffusion since MSD increases at a slower rate than in the case of standard diffusion. Such situations typically occur when waiting times between steps are included in the models of movement (e.g. Weeks et al. 1996; although a model with waiting times may not always be sub-diffusive, see §2.2), or if the spatial domain is constrained in some way (Coscoy et al. 2007), e.g. with the presence of a barrier (as in §2.6).

3. μ=1. This is the standard relation between MSD and time for diffusive movement.

4. 1<μ<2. This situation is known as super-diffusion since MSD increases at a faster rate than in the case of standard diffusion (although not so fast as with ballistic movement). Such situations typically occur when the step lengths in the walk are drawn from a distribution with infinite variance. Such a process is known as a Lévy walk and has been extensively studied by physicists and more recently ecologists (Viswanathan et al. 1996) although its general applicability to animal movement is still open to debate (Benhamou 2007; Edwards et al. 2007).

5. μ=2. In this situation, the movement is described as ballistic or wavelike, and MSD increases quadratically with t. This corresponds to the absolute displacement (cf. MDD) increasing linearly with time, which is a standard property of a wave process. In such cases, the characteristic backtracking and random movement associated with diffusive processes is not present, and each individual effectively moves in a straight line (in a random direction) away from the origin for the whole time period. (Note that, as seen is §2.2, MSD scales with t² for large t in a BRW; here, we are concerned only with an unbiased movement process.)

These five cases are the only possibilities that can be observed for an individual moving with a finite speed; μ>2 does not correspond to anything meaningful in this context. It should also be noted that an observed path may appear to belong to several of the above categories, depending on the time period of the observation and the spatial scale used. For example, correlated movement may appear ballistic for small t, but diffusive at large t (with highly correlated movement taking a longer time to appear diffusive than movement with low correlation). Hence, care should always be taken before classifying any particular path as belonging to one of the cases above (Coscoy et al. 2007).

2.8 Diffusion limit of the telegraph equation

In §2.1, we showed that the diffusion coefficient of a SRW is D=δ²/(2τ), where δ is the distance moved at each jump and τ is the time step between jumps. Since the telegraph turning process is a Poisson process of intensity λ, the mean time between turning events is and the average distance moved between turning events is . Hence the effective diffusion coefficient for the telegraph process is indeed given by

(2.29)

Box 3. Hyperbolic or parabolic?

A linear second-order PDE can be written in the form

where u(x, y) is a function of both x and y (which could represent, for example, space and time). If ac−b²<0 then the PDE is classed as hyperbolic. By contrast, if ac−b²=0 then the PDE is parabolic and if ac−b²>0 then the PDE is elliptic.Hyperbolic and parabolic equations have been extensively studied and there is a large body of theory available, which deals with solving such equations (and also examines their relationship with each other), mainly relating to problems arising in physics. However, in our case it is interesting to note that, using the above definitions, the drift–diffusion equation (2.8) derived from a position jump processes is parabolic, while the telegraph equations, (2.27) and (2.30), and other velocity jump processes (see §3.5) are hyperbolic.The relationship between the limiting processes involved in the movement model and the governing hyperbolic and parabolic equations has been discussed in more detail by, for example, Othmer et al. (1988), Hillen & Othmer (2000), Hillen (2002) and Othmer & Hillen (2002).

2.9 The biased telegraph equation

The derivation of the biased one-dimensional telegraph equation is similar to the unbiased case, except that we use different turning probabilities depending on the direction of movement. Denoting the probability of turning by r₁=λ₁τ for right-moving individuals and r₂=λ₂τ for left-moving individuals, it can be shown that the governing equation, called the biased telegraph equation, is

(2.30)

Note that, in this example, the bias is introduced through the different rate of turning in each direction, which is a form of klinokinesis; this contrasts with the uncorrelated BRW in §2.2, where the bias came about through the probability of moving in each direction, which is a form of taxis (see §3.8).

2.10 The telegraph equation in higher dimensions

In §2.5, the solution to the two-dimensional diffusion process with or without drift was shown to be valid only as a long-time approximation. One can introduce correlation by completing a similar derivation as in §2.7, but working with a two-dimensional lattice rather than a line. There are now four possible directions of movement: right, left, up and down. Initially, we assume that the probability of turning is independent of the direction of movement (so there is no bias). As before, we assume a constant speed v. We split the population into individuals moving in each of the four directions α₁, …, α₄. At each time step τ, an individual can turn π/2 rad anticlockwise or clockwise (with probabilities λ₁τ and λ₂τ, respectively), turn π rad (with probability λ₃τ), or continue in the previous direction (with probability 1−(λ₁+λ₁+λ₃)τ).

Completing a similar analysis to §2.7 leads to a set of differential equations for each of the sub-populations moving in the four directions. Further manipulation similar to that in §2.7 then leads to

(2.31)

3.1 Mean squared displacement of CRWs

The MSD, defined in N dimensions by (2.13), gives a measure of the spatial spread of the population with time and, owing to its relationship with the diffusion coefficient D via, for example (2.14), is of great importance to those studying dispersal in biological systems (Okubo & Levin 2001). Interestingly, however, many of the results discussed in this section were first derived through studies in molecular chemistry. The growth and space-filling properties of polymer chains and larger molecules have been modelled as a CRW by, for example, Tchen (1952) and Flory (1969), who both derived equations for MSD. Tchen (1952) also demonstrated the important result that the location coordinates after a large number of steps of a CRW are normally distributed.

Taylor (1921) derived the following expression for the MSD in a one-dimensional correlated walk (see also Tchen 1952; Flory 1969; Hanneken & Franceschetti 1998; Okubo & Levin 2001). Suppose the walker takes a series of steps y_j (j=1, …, n) of constant length (|y_j|=δ). A correlation is explicitly introduced between the directions of successive steps (although this correlation propagates in a gradually diminishing way through the Markov process). As the number of steps n becomes large, the MSD for this discrete process tends asymptotically to

(3.1)

Kareiva & Shigesada (1983) were the first to set up and analyse a generalized model of a two-dimensional CRW for animal movement that included a variable step length and a general angular distribution for the direction moved at each step (see also Skellam 1973; Nossal & Weiss 1974; Lovely & Dahlquist 1975; Hall 1977; Dunn 1983; Marsh & Jones 1988). The CRW consists of a series of discrete steps of length L_j and direction Θ_j. The length L_j of the jth move and the turning angle Φ_j=Θ_j+1−Θ_j are assumed to be random variables with no autocorrelation or cross-correlation (and no correlation between step length and step direction). The CRW thus consists of a series of independent draws from the step length PDF, p(l), and the turning angle PDF, g(ϕ) (box 1), for each step (i.e. the process is a first-order Markov chain, see Grimmett & Stirzaker (2001)). We define the mean cosine c and mean sine s of the turning angle as (see also box 1)

(3.2)

Using this model, Kareiva & Shigesada (1983) derived the following equation for the MSD after n steps

(3.3)

The general result (3.3) reduces to a much simpler form in particular cases. For example, if there is no persistence (i.e. the walk is uncorrelated), then g(ϕ) has a uniform density and both c and s are zero, and (3.3) reduces to (cf. (2.14) with no drift, u=0 and with N=2). As mentioned, a more realistic case is when organisms exhibit equal probabilities of turning clockwise or anticlockwise, so g(ϕ) is symmetric about ϕ₀=0 (although it is worth noting that this is not always the case—for example, bacteria can exhibit an inherent rotational bias due to the handedness of their flagellar motor). In the case of equal turning probabilities, we get s=0 and (3.3) may be written in terms of the coefficient of variation b of the step length L (b₂=E(L²)/(E(L))²−1)

(3.4)

In the case of random walks with bias, it can be more complicated to derive expressions for MSD. It is possible to write down the MSD of an uncorrelated BRW using a similar approach to the above. Consider a walk consisting of steps of length L_j and direction Θ_j (note the important difference between the direction of movement Θ_j and the turning angle Φ_j=Θ_j+1−Θ_j), whose mean sine is zero and whose mean cosine is q. After n steps, the MSD is given by (Marsh & Jones 1988; Benhamou 2006)

(3.5)

An equation for MSD in a biased velocity jump process can be generated for both the one-dimensional case and higher dimensions, using a generalized transport equation (see §3.5). However, these are special cases of a BCRW where the turning events occur as a Poisson process. In the case of a BCRW with a fixed time step between turning events (where the spatial step lengths can be either fixed or variable), it remains an open problem to calculate a direct equation for MSD.

3.2 Mean dispersal distance of unbiased CRWs

The MSD, defined in (2.13), is of interest to ecologists and biologists owing to its relation to the diffusion coefficient D via (2.14). The MSD is also a statistic that is reasonably mathematically tractable, as illustrated in §3.1. However, because the MSD deals with the squared dispersal distance, it has been suggested that a more intuitive statistic is the mean dispersal distance (MDD; Bovet & Benhamou 1988; McCulloch & Cain 1989; Wu et al. 2000; Byers 2001). The MDD of a dispersing population is defined in N dimensions as

(3.6)

Consider a two-dimensional, unbiased, CRW with mean step length E(L), and a symmetrical about zero (i.e. zero mean sine) probability distribution g(ϕ) for the turning angle at each step. After a sufficiently large number of steps n, the location coordinates X_n and Y_n are independently normally distributed with equal variance. Bovet & Benhamou (1988) used this to derive the following approximate relationship between MDD and MSD:

(3.7)

As discussed in §2.7, it is possible to derive differential equations for the MSD of a BCRW where the times between turning events are distributed as a Poisson process but, owing to the presence of an absolute value in (3.6), this approach cannot be used to derive an equation for the MDD of a BCRW. This remains an open problem but, for particular cases, it may be possible to use the result of Bovet & Benhamou (1988), together with an equation for the MSD derived from a moment closure method (see Codling 2003).

3.3 Tortuosity of CRWs

The tortuosity of a path describes the amount of turning in a given space or time. Clearly, tortuosity is related to the MSD and MDD: highly tortuous paths will spread out in space slowly (small MSD), while straight paths will spread out in space quickly (high MSD). Hence, it can be useful to measure and study the tortuosity of observed paths in order to understand the processes involved, estimate the area searched by an organism and predict spatial dispersal. Several measures of tortuosity are available but most have some limitations.

The straightness index (sometimes called the net-to-gross displacement ratio) is a relative measure that compares the overall net displacement G of a path with the total path length T (Batschelet 1981). For example, if a random walk starts at location (x₀, y₀) and, after n steps with lengths l_j (j=1, …, n), ends at (x_n, y_n), then the straightness index is given by

(3.8)

Unfortunately the straightness index is not a reliable measure of tortuosity of a CRW, because the mean of G corresponds to the MDD and hence increases with the square root of n. Consequently, the ratio G/T tends to zero as the number of steps increases (Benhamou 2004). Thus, when observing a CRW and measuring the straightness index, a different result will be found depending on the number of steps considered (the total time observed or the total path length used). It is, therefore, very difficult to compare the tortuosity of different CRWs using this method, unless they all consist of a similar number of steps.

The tortuosity of a CRW corresponds to the amount of turning associated with a given path length and, in some way, measures its long-term diffusion potential. On this basis, Benhamou (2004) defined the path sinuosity in terms of the mean step length and MSD after a large number of steps as . Using equation (3.4) simplified by removing the second term on the right-hand side, which becomes negligible in the long term, gives

(3.9)

Equation (3.9) can be applied directly to actual paths when animals naturally move in a discrete way (e.g. a bee flying from one flower to another), provided the basic requirements of a CRW are respected (independence of step lengths and turns). For animals moving in a continuous way, the path is usually discretized when recorded (i.e. the raw data consist of a set of locations rather than of a continuous track), and can be rediscretized for analysis purposes. (Re)discretization has been shown to alter the turning angle and step length distributions, and hence involves the use of corrected formulae (Bovet & Benhamou 1988; Benhamou 2004; Codling & Hill 2005a). Note that rediscretization can be done either spatially (Bovet & Benhamou 1988), which is easier to fit with the analysis above, or temporally (Codling & Hill 2005a), which is perhaps more natural for experiments in the field.

It has also been proposed to measure the tortuosity of animals' random search paths by a fractal dimension. This would be useful for actual paths that can be reliably represented using fractioned Brownian motion, which is a fractal movement model where the persistence is related to the fractal dimension through the use of a parameter called the Hurst coefficient. To our knowledge, however, it has not yet been shown that this kind of model can provide a better representation of animals' random search paths than CRW, whose true fractal dimension is equal to 2 for two-dimensional movements (i.e. they will eventually fill the entire plane) for any positive sinuosity (for S=0, the CRW reduces to a straight line, with fractal dimension equal to 1). Applying the classical ‘divider’ method, initially developed to measure the fractal dimension of fractal lines, the CRW provides pseudo-fractal dimension values (Turchin 1998; Benhamou 2004), which correspond to indirect measures of the mean cosine of turns c (Nams 1996, 2005). These pseudo-fractal values cannot reliably estimate the path tortuosity of a CRW, which depends not only on the mean cosine of turns c, but also on the mean step length E(L) and coefficient of variation b.

3.4 Bias in observed paths

The next step up in model complexity from a CRW is a BCRW. We will deal with modelling approaches for BCRWs in §§3.5 and 3.6, while in this section we will consider the ways in which bias may be detected in an observed path (see also Coscoy et al. 2007). Assuming the simplest possible environment and behaviour, there are two main ways in which individuals may respond to a signal and hence introduce bias into their movement. Firstly, there may be a fixed sensory gradient such that the preferred direction is always the same for all individuals at all locations in space (e.g. micro-organisms moving under the influence of gravity). Secondly, there may be a target or point source fixed in space, such that the preferred direction is towards a fixed point (e.g. animals searching for a food source). Note that the first scenario can be considered as a special case of the second where the target is fixed at infinity.

In order to parametrize continuum models of bioconvection in micro-organisms (Kessler 1986; Hill & Pedley 2005), Hill & Häder (1997) analysed the paths of swimming algae undergoing gyrotaxis (upward swimming due to a gravitational torque) or phototaxis (swimming towards a light source). In both situations, it was assumed that the preferred absolute direction of movement was independent of location (i.e. a ‘target at infinity’), but that the turning angle was dependent on the most recent direction of movement. To help inform experimental observations, Hill & Häder (1997) set up a random walk on the unit circle as follows: at each time step τ, the walker makes a small turn δ clockwise with probability a(θ) or anticlockwise with probability b(θ), or continues in the same direction with probability 1−a(θ)−b(θ). Following a similar method to that used in §2.2 and taking the limit δ, τ→0 such that δ²/τ is constant yields the Fokker–Planck equation for the PDF p(θ, t) of travelling in direction θ at time t

(3.10)

By rediscretizing the observed paths and binning the data, Hill & Häder (1997) estimated using a similar method to Bovet & Benhamou (1988; as is essentially the same measure as sinuosity, see §3.3), and also the functional form of μ₀. For gyrotaxis, a sinusoidal reorientation model was found to fit well: μ₀(θ, t)=−sin (θ−θ₀(t))/B, where θ₀(t) is the preferred direction and B is the average time taken to reorient to the preferred direction. If both and θ₀ are assumed to be constant then the steady state of the Fokker–Planck equation (3.10) is a von Mises distribution (box 1). This sinusoidal response was predicted by Kessler (1986) and can be explained by the gravitational torque that acts on the individual alga: the torque is zero if the alga is moving in the preferred direction or in the opposite direction, and is the greatest when the alga is moving perpendicular to the gravitational force. For phototaxis, a better fit was given by a simple linear response: μ₀(θ, t)=−(θ−θ₀(t))/B. This case is more generally applicable since there are only a few physical situations likely to result in a sinusoidal response. However, if both and θ₀ are allowed to be variable, then the steady state of the Fokker–Planck equation (3.10) is more complicated (Hill & Häder 1997).

This ad hoc method used by Hill & Häder (1997) was tested using simulations by Codling & Hill (2005a) and found to be valid (allowing for smoothing errors not accounted for by Hill & Häder (1997) and assuming that sinuosity is low enough to avoid ‘wrapping’ problems leading to an underestimate of (Bovet & Benhamou 1988; Benhamou 2004)). Note that this method allows for both the functional form of the bias response to be estimated and the relevant parameters (mean reorientation time and angular diffusivity) to be quantified.

The situation is more difficult when the preferred direction changes with spatial location (and is thus individual dependent), a common scenario when animals are moving towards a target in space. The main problem in such a scenario is distinguishing between localized bias due to forward persistence and true biased movement towards a target. This is particularly true if there is more than one target or if the target moves in space. Benhamou (2006) suggested a new procedure based on the backward evolution of the beeline distance from the end of the path (the goal) to each animal's preceding locations. This procedure is efficient, as it requires only a small sample of short paths for detecting a possible orientation component, but is not perfect as there remains a relatively high (approx. 30%) probability of misidentifying a true CRW as a BCRW. This type I error can be reduced by considering together a number of paths assumed to be of the same kind.

The above scenarios assume that bias is introduced through reorientation towards the preferred direction. However, there are many other ways to model mechanisms that produce a directional drift, including the case where the mean turning rate is null and the directional bias comes about through variations in (see Benhamou 2006, and §3.8). For example, the classical run-and-tumble behaviour observed in chemotactic bacteria is usually modelled through a low rate of turning when moving in the preferred direction (runs) and a high rate of turning (tumbles) otherwise (Berg 1983). We discuss the relation between these various mechanisms in §3.8.

3.5 The transport equation and general hyperbolic models of movement

In §2.7, it was shown how the telegraph process could be used in one dimension to model a CRW and how this resulted in a hyperbolic governing equation (box 3). A similar process does not produce a closed equation in two or more dimensions, but it is still possible to work with a generalized hyperbolic governing equation to model the CRW and BCRW (velocity jump processes) in higher dimensions. Othmer et al. (1988) introduced the idea of a governing mass balance equation, the linear transport equation, which can be used to describe the movement and reorientation of cells (and animals; see also Alt 1980).

Let p(x, v, t) be the density function for individuals moving in N-dimensional space, where x∈^N is the location of an individual and v∈^N is its velocity. The total number density of individuals at location x, regardless of velocity, is given by integrating p over all possible velocities

(3.11)

(3.12)

In general, we are interested in the first few velocity moments (as from these we can calculate the statistics of interest such as E(X_t) and E(R_t²)), including the number density n(x, t) introduced in (3.11), and the average velocity u(x, t), which is defined by

Integrating (3.12) over v gives an evolution equation for n in terms of u

(3.13)

(3.14)

It is worth noting that equation (3.12) describes a general process, of which many of the basic random walk models discussed previously are special cases. For example, in one spatial dimension, and under the assumption that individuals move with constant speed v, there are only two possible velocities: +v and −v. With the additional assumption that the turning frequency λ is constant, it can be shown that (3.13) and (3.14) reduce to equations (2.25) and (2.26) derived in §2.7, which lead to the one-dimensional telegraph equation.

The transport equation also provides a natural extension of the basic telegraph process to two dimensions, without the need to restrict the population to a lattice as in §2.10. Retaining the assumptions of constant turning frequency λ and constant speed v, an individual's velocity may be described simply by the angle θ between its direction of motion and the positive x₁-axis. The appropriate density function is now p(x, θ, t). The transport equation (3.12) reduces to

(3.15)

Othmer et al. (1988) used the following example to illustrate how a random walk in an external field can be modelled using (3.15). Suppose individuals are moving with a taxis-inducing gradient in the direction θ₀=0 (i.e. the positive x₁-direction), under the assumption that the gradient influences only the reorientation kernel T(θ, θ′). Suppose also that the reorientation kernel T(θ, θ′) is the sum of a symmetric probability distribution h(ϕ), where ϕ=θ−θ′, and a bias term k(θ) that results from the taxis-inducing gradient. Since the gradient is directed along the x₁-axis, the bias term takes its maximum at θ=0 and is symmetric about θ=0. Thus

(3.16)

Assuming that all individuals start at the origin (0, 0), with initial directions uniformly distributed around the unit circle, one can derive a system of differential equations, known as moment equations, for the statistics of interest. For the above choice of T(θ, θ′), this system is straightforward to solve (Othmer et al. 1988; Codling 2003). The mean velocity and mean location are given by

The choice of reorientation kernel in (3.16) is crucial in the above analysis as it results in a closed system of differential equations for the statistics of interest. Codling & Hill (2005b) used an arguably more realistic reorientation kernel based on a single symmetric distribution with mean turning angle dependent on the direction of movement (see §3.4). However, this results in a ‘cascade’ of higher moment equations and further assumptions to close the system need to be made.

A general discussion of the issues relating to the moment closure of systems of moment equations resulting from transport equations is given in Hillen (2002). The transport equation can be developed further to allow more complex scenarios to be modelled (spatially dependent parameters, etc.) although this is likely to make finding a closed-form solution more difficult. The hyperbolic movement model has been used in the place of diffusion models to create reaction–transport systems (Hillen 1996, 2002), which can be used in place of the classical reaction–diffusion models (Turing 1952; Murray 1993) used in pattern formation and developmental biology.

3.6 Reinforced random walks

The types of processes whereby walkers modify the chemical environment of themselves and of other individuals in the population have led modellers of cell locomotion to employ the theory of RRWs, which allows the walker to modify (reinforce) the transition probabilities associated with the grid points, or interval, it traverses (Davis 1990). The most common way of representing a one-dimensional RRW is the so-called master equation

(3.17)

There are many possible models for the transition rates T^± in terms of the concentration w(x, t) of a control substance, such as the ‘local model’, ‘barrier model’ and ‘normalized barrier model’ proposed by Othmer & Stevens (1997). For a specific choice of transition rates, the continuum limit of the master equation (3.17) can often be found by a similar method to that used to obtain equation (2.5), and may usually be written in the form

(3.18)

In the local model of Othmer & Stevens (1997), the transition rates depend only on the local concentration of control substance: for some constant ρ and positive function F, called the transition probability function. Taking the limit δ→0 and ρ→∞, such that D=ρδ² is a constant, leads to a continuum limit equation (3.18) with diffusivity d(w)=DF(w) and chemotactic sensitivity Χ(w)=−D(dF/dw). Hence, as well as modulating the rate of random motility, the control substance provides a directional bias because the cell is more likely to move from an area of high motility to low motility than vice versa.

In the barrier model, the transition rates depend on the concentration of control substance in the interval to be traversed: . Hence, the cell can sense concentrations at x±δ/2 (so δ is the effective sensory range of the cell) and uses these to decide its direction of movement. The continuum limit equation is (3.18) with diffusivity d(w)=Dτ(w) and chemotactic sensitivity Χ(w)=0. Hence, the control substance has an effect of random motility, as in the local model, but there is no directional bias. This is expected because the transition probability from x to x+δ, determined by the control substance concentration at x+δ/2, is the same as the transition probability from x+δ to x, so there is no directional bias.

In the normalized barrier model, the transition rates are as in the barrier model, but are normalized so that T⁺+T⁻ is a constant,

(3.19)

The actions of the control substance under each of these three models are summarized in table 1. Note that the normalized barrier model represents a control substance that is an attractant or repellent, whereas the barrier model represents an inducer or inhibitor. The local model allows both the motility and the directional bias to be varied. However, the constraint on the relationship between d(w) and Χ(w) in the local model means that the substance must either be an inducer and a repellent, or an inhibitor and an attractor, which is not always realistic. Motivated by this limitation, the normalized barrier model may be generalized to allow for a variable mean waiting time (VMWT). This may be achieved in two ways. Firstly, the transition rates take the form (3.19), but ρ is permitted to be a function of the control substance concentration w (Plank 2003). This results in a diffusivity of d(w)=Dλ(w) and a chemotactic sensitivity of , thus allowing both non-directional (random) and directional effects of the control substance to be specified independently. A second possible method is to combine the normalized and unnormalized barrier models by taking the transition rates as

(3.20)

Partial differential equations of the form (3.18) have been studied analytically for certain choices of the functions Χ(w) and d(w) and of the dynamics for w. For example, Keller & Segel (1971) obtained travelling wave solutions and Rascle & Ziti (1995) obtained similarity solutions in the case where chemotactic sensitivity is inversely related to w, and the control substance is consumed linearly by the cells. Othmer & Stevens (1997) considered three main cases of cell behaviour, which they termed ‘aggregation’ (meaning solutions in which the cell density aggregates to a finite positive value in certain regions), ‘blow-up’ (cell density tends to infinity in finite time) and ‘collapse’ (cell density tends to zero everywhere). Analogous results for the discrete RRW were originally shown by Davis (1990), and were linked to the continuum-level results by Plank (2003). Analytical conditions for the cases of aggregation, blow-up and collapse were developed by Levine & Sleeman (1997) and further theoretical work on mixed parabolic–hyperbolic equations was done by Takase & Sleeman (2002).

3.7 Applications of reinforced random walks

The theory of RRWs has been presented in terms of one-dimensional walks for simplicity, but is straightforward to extend to higher dimensions. Furthermore, the actions of multiple control substances (w₁, …, w_n) can be modelled by taking the transition probability function τ to be of the form F(w)=F₁(w₁) … F_n(w_n) (Levine et al. 2001). Under the normalized barrier model, this leads to a sum of n chemotactic terms in the continuum limit equation (3.18).

RRWs have been used to model the movement of myxobacteria (Othmer & Stevens 1997) and the migration of endothelial cells during tumour-induced angiogenesis (growth of new blood vessels; Levine et al. 2001; Sleeman & Wallis 2002; Plank & Sleeman 2003). These studies adopted the normalized barrier model, chiefly because this models a control substance that offers a directional stimulus (which is of prime importance for successful angiogenesis), but has no direct effects on cell motility. Nevertheless, the growth factors regulating angiogenesis do also influence cell motility and an extension to a VMWT model, which allows for both directional and non-directional effects, was proposed by Plank et al. (2004). A non-lattice RRW, similar to that outlined in §3.4, was subsequently used (Plank 2003; Plank & Sleeman 2004) to model angiogenesis in two and three dimensions. In addition to choosing an appropriate functional form for the reorientation rate μ₀(θ), one must define the preferred direction θ₀. In the case of a chemoattractant, the natural choice is to define θ₀ to be in the direction of increasing chemical concentration (i.e. in the same direction as ∇w), which may be a function of the cell's location x and time t.

In models such as these, the control substance(s) may be assumed to be in a steady state, or the random walk model for cell movement may be coupled to one or more time-dependent PDEs for the concentration of the control substance(s). If the cells modulate their chemical environment by producing or consuming control substances, these PDEs will contain source/sink terms depending on cell density. In the case of angiogenesis, the basic cellular motion driving capillary growth is modelled by a random walk, but further rules must be specified to allow for capillary branching and looping (anastomosis). In addition to the transition probabilities, each cell has a probability, which may depend on local control substance concentration, of branching to create two daughter capillaries. If a capillary tip collides with an existing vessel, an anastomosis is formed and the colliding cell ceases to take part in the random walk.

An alternative to the RRW approach described above is to take a continuum Fokker–Planck equation of the form (3.18) as the starting point. This may be discretized using standard techniques to obtain a finite-difference equation of the form (2.4) (with the obvious extension to higher dimensions), which can be used to identify the relevant transition probabilities of moving left and right, and staying still. This method was developed by Anderson et al. (1997) to model nematode movement and has since been used to study tumour angiogenesis (Anderson & Chaplain 1998) and tumour cell invasion (Anderson et al. 2000), and to distinguish between the inhibitory and repellent effects of a signalling molecule called ‘Slit’ in experimental work (Cai et al. 2006). The random walk model for angiogenesis has also been coupled to a flow model to study the effects of blood flow and nutrient delivery in a nascent capillary network (McDougall et al. 2006).

Both the ‘bottom-up’ approach of the RRW and the ‘top-down’ approach of discretizing a PDE have the same governing equation (3.18) in the continuum limit. However, the two methods will, in general, result in random walks with different transition probabilities, illustrating the fact that, usually, there is not a unique random walk model corresponding to a given continuum equation. The RRW method has the advantage that the transition probabilities are derived mechanistically from the underlying biology, using a transition probability model of the type summarized in table 1, rather than via a mathematical discretization and normalization procedure (see Plank & Sleeman 2004, for details).

3.8 Biological orientation mechanisms

There are two main types of mechanisms for movement in response to a stimulus. Kinesis refers to the situation where the organism samples only the stimulus intensity at a single point and modulates its speed of movement (orthokinesis or O-kinesis) or its path sinuosity (klinokinesis or K-kinesis) accordingly. The terms O-kinesis and K-kinesis were first defined by Gunn et al. (1937) in terms of linear and angular speeds, but the redefinition of K-kinesis in purely spatial terms (Benhamou & Bovet 1989), by using sinuosity rather than rate of change of direction, has led to a much clearer view of the properties of kineses. By contrast, taxis is where the organism is able to detect a preferential direction of movement, and bias its turns accordingly, without necessarily altering its overall speed of movement or rate of turning.

The random walks described above are all couched in terms of transition rates, which are implicitly assumed to be under the full control of the migrating organism. This raises questions concerning the information that the organism requires in order to be able to exert this control, and the probable realism of assuming that such information is available. In the case of the barrier-type models (see §3.6), it is clear that the walker needs at least two sensors in order to compute the transition rates ‘on the spot’ (i.e. without moving); this type of mechanism is referred to as tropotaxis. The same is also true for the non-lattice random walk (§3.5) if the walker needs to resolve the target direction θ₀, for example by measuring ∇w. The same taxis may also be produced if the organism has only one sensor (or many sensors that are too close together to allow it reliably to detect differences), but moves its body in various directions to sample the local variations in stimulus; this is termed klinotaxis. In its simplest form, taxis does not incorporate a correlation between successive step directions, and hence corresponds to a BRW: individuals directly settle their local direction with respect to the target direction. More realistic taxis models can be obtained in the form of BCRW by introducing such a directional correlation (Benhamou & Bovet 1992; Benhamou 1994). The result of the taxis is then determined by the relative weights of the local directional bias (persistence, which controls the motility) and the global directional bias (goal attractiveness, which controls the advection), and by the level of random noise in the system.

By contrast, the local model in §3.6 comes under the category of kinesis, as opposed to taxis, because the walker measures only the stimulus intensity at a single point in order to compute the transition rates. Although the probabilities of moving in different directions are always equal at any given time, this mechanism results in an effective directional bias towards areas of low motility, illustrating the fact that a kinetic mechanism can produce a directional bias.

There are two working modes of kinesis (Benhamou & Bovet 1989): absolute (A) and differential (D). The movement is controlled, in A mode, with respect to the local stimulus intensity experienced at any given location or, in D mode, with respect to the change in stimulus intensity experienced during a step (e.g. through sensory adaptation). Both AO-kinesis and AK-kinesis are space-use mechanisms for exploiting patchy environments, enabling the animal to reduce locally the diffusion coefficient (D=v/S²) of its movement (by decreasing mean speed v or increasing the sinuosity S) in the most suitable areas (Benhamou & Bovet 1989). Unlike AO-kinesis, however, AK-kinesis is able, when applied to a gradient field, to generate a slight drift in the gradient direction based on simple differences in angular diffusivity , with a null mean turning rate (μ₀=0). This property was used by Jamon & Bovet (1987) to account for the homing behaviour of mice. However, the advection component is very weak (and accordingly most animals get lost), so that AK-kinesis cannot be considered as an efficient orientation mechanism. By contrast, DK-kinesis can reach 90% of the efficiency of a pure taxis (BRW; although taxis with persistence can outperform it in a noisy environment because persistence can be used to smooth the noise (Benhamou & Bovet 1992)). Indeed the so-called chemotaxis of bacteria (Alt 1980; Berg 1983) is mainly a DK-kinesis where the sinuosity is modulated through step length rather than turning angle variance (called a run-and-tumble mechanism). Finally, DO-kinesis seems to have no biological applications.