Modelling nematode movement using time-fractional dynamics

ARTICLE IN PRESS

Journal of Theoretical Biology 248 (2007) 212–224 www.elsevier.com/locate/yjtbi

Modelling nematode movement using time-fractional dynamics Simona Hapcaa,, John W. Crawforda, Keith MacMillanb, Mike J. Wilsonb, Iain M. Younga a

SIMBIOS, University of Abertay Dundee, Dundee DD1 1HG, UK School of Biological Sciences, University of Aberdeen, Aberdeen AB24 3UU, UK

b

Received 19 December 2006; received in revised form 3 May 2007; accepted 3 May 2007 Available online 10 May 2007

Abstract We use a correlated random walk model in two dimensions to simulate the movement of the slug parasitic nematode Phasmarhabditis hermaphrodita in homogeneous environments. The model incorporates the observed statistical distributions of turning angle and speed derived from time-lapse studies of individual nematode trails. We identify strong temporal correlations between the turning angles and speed that preclude the case of a simple random walk in which successive steps are independent. These correlated random walks are appropriately modelled using an anomalous diffusion model, more precisely using a fractional sub-diffusion model for which the associated stochastic process is characterised by strong memory effects in the probability density function. r 2007 Elsevier Ltd. All rights reserved. Keywords: Correlated random walk; Stochastic process with memory; Conditional probability; Mean square displacement

1. Introduction The study of movement patterns has been productive in explaining a range of ecological questions related to foraging strategies (Bell, 1991), dispersal behaviour (Stamps et al., 1987), and community interactions (Banks et al., 1987). However, because detailed studies of organism movement over long distances are technically difficult, the need to develop quantitative methods for the description and prediction of such movements is important (Levin, 1987). Phasmarhabditis hermaphrodita is a slug parasitic nematode capable of killing several species of slugs and in particular Deroceras reticulatum, the most widespread pest slug in UK (Wilson et al., 1993). Therefore, understanding their dispersal abilities has practical importance for biological control of pest species (Shapiro and Glazer, 1996) and contributes to knowledge of the dynamics of soil animal communities. There are very few models of nematode movement in the literature. Croll and Sukhdeo (1981) treat the nematode as a random walker and Anderson et al. (1997b) used a linear diffusion equation, perturbed by an additional term to simulate the foraging Corresponding author. Tel.: +44 1382 308695; fax: +44 1382 308117.

E-mail address: [email protected] (S. Hapca). 0022-5193/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2007.05.002

behaviour. Even in the wider context of modelling animal movement, diffusion models, assuming a random movement pattern of the organism, have been a common starting point for efforts to quantify the animal dispersion in different environments. The basic assumption of this model being that the animal moves in a completely random manner with no correlation between successive steps, this leads to a uniform distribution of turning angles and independence between successive step lengths. For large scale, sometimes the pure random walk model can be appropriate (Kareiva and Shigesada, 1983), but in many other cases it failed due to the simplistic assumption about the movement pattern of individuals, which does not take into account the fact that many entities, including animals, exhibit correlated movement (Bovet and Benhamou, 1988). A random walk is said to be correlated when the turning angle between successive steps in a path is not uniformly distributed in the interval (p, p), such that the subsequent direction of movement depends on the previous directions. The simplest example of correlated random walk is given by a very short memory, i.e., the step i depends only on the previous step i1 and is independent of the other previous steps; therefore, two successive turning angles are independent and hence uncorrelated. When the turning angle distribution has a bias in either direction, the correlated

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

random walk presented above can be used to model the movement pattern characterised by loops (Bell, 1991; Blanche et al., 1996; Anderson et al., 1997a; Conradt et al., 2000; Bengtsson et al., 2004). In many applications this leads to an oversimplification, since the resulting loops are orientated in one direction only. A solution to overcome this problem is to add more memory to the correlated random walk, which was done by introducing correlation between the turning angles (Blanche et al., 1996; Wu et al., 2000; Preisler et al., 2004). Wiktorsson et al. (2004) went further and added correlations between turning angles and the variable speed. By using Cartesian coordinates instead of polar coordinates, they presented a linear Gaussian autoregressive model, of p order, to describe the movement of a soil-living insect, Prophorura armata, where the parameter p expresses the memory of the process. The present work uses digital recordings of the movement of slug parasitic nematode, P. hermaphrodita as a point of departure for the development of a mathematical model that quantifies the nematode movement in homogeneous environments. More precisely, the model aims to provide the probability density function (PDF) and the associated diffusion equation in order to predict the nematode dispersal over time and space. Image analysis of nematode trails shows that, there is a significant individual variation within the population, from nematodes that move very fast and make large loops to nematodes that move slowly and often fluctuate between backward and forward movement. By computing the correlation between different steps, we also found that nematode movement is characterised by a long-range correlation, in terms of step lengths and turning angles. In other words, the subsequent step depends on the previous steps. Therefore, the starting point of the model is a 2d correlated random walk on a rectangular grid to simulate nematode trail, and the associated long time memory stochastic process, whose PDF is computed using a large number of trail simulations. The profile of this predicted probability function is compared with some theoretical diffusion models, whose parameters are derived from data computing the mean square displacement (MSD) of the real movement. 2. Materials and experimental methods

213

2.2. Experimental treatments and image capture The experimental arena consisted in a 9-cm-diameter Petri dish containing 1.2% technical agar (Oxoid Ltd., Hampshire, UK). One active nematode was picked from the water where it had been rehydrated, and transferred, by means of an eyelash mounted on a syringe pick, to the centre of the agar plate. An Axio MRc Zeiss camera attached to a Leica microscope MZ16 and connected to the computer monitored the movement of nematode with an image field size of 3 cm  4 cm. The position of the nematode was recorded at 8-sec intervals over a period of 15 min, which was equivalent to 100 frames. Preliminary investigations revealed that the nematode activity on the Petri dish remained constant for about an hour. The experiments were done under a constant temperature of 19 1C and 12 replicates were prepared. 2.3. Image analysis Images were processed using Axio Vision 3.1 software, then analysed using image analysis software incorporating a tracking algorithm (Image-Pro PlusTM v.5, Media Cybernetics Inc., Maryland, USA), in order to obtain the x and y coordinates at each time point. MatLab v.7.3 was used to plot the final digitised nematode trail and to compute the distance travelled between frames, the value of the turning angle, and the MSD from the point of departure as a function of time. Minitab v.14.13 was used to compute the correlation between successive steps. The nematode speed was determined by calculating the average distance between successive positions, and divided by the 8-sec time interval. An 8-sec interval was chosen on the basis that in this time the nematodes travel an average distance of 1 mm: approximately equal to a body length. Shorter timescales would have introduced artifactual error due to changes in the body shape during movement, and longer timescale would have reduced the amount of data available for analysis. Next, the turning angle distribution was generated by computing the angle between successive steps of the trail. The angular range chosen was from p radians to p radians, and all angles were computed with respect to the previous angle. Negative values correspond to the right turns and positive values to left turns.

2.1. Organisms used 3. Model description Dauer larvae of P. hermaphrodita were obtained from Becker Underwood (Littlehampton, UK) in a formulated product (Nemaslugs) consisting of partially dehydrated dauer larvae bound by clay particles. One hour prior to the experiments, the nematodes were added to water and cleaned by repeated sedimentation and washing. Typical body dimensions of the neatodes used were 1 mm in length and 0.03 mm in width.

3.1. A stochastic process associated to nematode movement We denote by O the set of nematodes whose dispersal on the plate is investigated, and by (Xt, tX0) a twodimensional stochastic process on O defined by Xt(o) ¼ (Xt1(o),Xt2(o)), the Cartesian coordinates in 2d of the position of nematode oAO at time tX0. Nematode activity

ARTICLE IN PRESS 214

S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

on the Petri dish was observed to remain constant during the recording time, as the nematode speed on the first half (s1 ¼ 0.108 mm/s) of the recording time was not significantly different to the nematode speed on the second half (s2 ¼ 0.11 mm/s) of the recording time (t-test, P ¼ 0.87). Therefore, we assume that the process (Xt, tA(0, T) has stationary increments, where T is a positive real value denoting the time the nematode activity on the Petri dish remains constant. We aim to compute the PDF pt(x),xAR2, tA(0, T), of the stochastic process XT, and to identify the associated diffusion equation. A very useful tool in this direction is given by the MSD of the process (Wiktorsson et al., 2004), which is a function of time defined by MSDðX t Þ ¼ EðjjX t  X 0 jj2 Þ X 1 jjX t ðoÞ  X 0 ðoÞjj2 ; ¼ cardðOÞ o2O

(a) correlation between the step lengths rðRn ; Rm Þ ¼ rR ðn  mÞ, (b) correlation between the turning angles: rðyn ; ym Þ ¼ ry ðn  mÞ, (c) cross-correlations between the turning angle and the step length: rðyn ; Rm Þ ¼ ry;R ðn  mÞ or rðRn ; ym Þ ¼ rR;y ðn  mÞ,

ð1Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where for any x 2 R2 ; jjxjj ¼ ðx1 Þ2 þ ðx2 Þ2 ; denotes the norm of x in R2. If the pattern of the MSD has a power-law form: MSDðX t Þ  K a ta ,

where mon, there are four correlation coefficients associated to a 2d stochastic process:

(2)

then the process Xt generates an anomalous diffusion (Metzler and Klafter, 2000, 2004) of sub-diffusion type for 0oao1 and super-diffusion type for 1oao2. If a ¼ 1 then MSD follows a linear pattern in which case Xt is the Brownian motion. Examples of anomalous diffusion following a power-law pattern that are considered in the paper as potential candidates for our model are the timefractional dynamics (TFD) and the fractional Brownian motion (FBM). Details on these two stochastic processes can be found in Appendix A. 3.2. Correlation between steps—general theory Let t40 be a real value corresponding to 8-sec real time, the time period between successive frames then tn ¼ n  t; n ¼ 1; 2; . . . ; 100 represents the time points for which the nematode position X n ¼ X tn was recorded. Denote by Y n ¼ X n  X n1 ; n ¼ 2; . . . ; 100 the steps of the process Xt corresponding to the time interval (tn1, tn) of length t, and let us introduce the variable step length Rn ¼ jjY n jj; n ¼ 2; . . . ; 100, being the distance travelled during the interval (tn1, tn) and the variable turning angle yn, n ¼ 3,y,100 between successive steps. As the process Xt has stationary increments we deduce that for all n ¼ 2,y,100, the steps Yn have the same distribution. Consequently, all the step lengths Rn and all the turning angles yn have the same distribution. If the steps Yn are independent, then Xn is a Markovian process. Otherwise, Xn is called a stochastic process with memory, presenting correlations between steps. Therefore, a method to decide on the nature of the process Xn, is to compute the correlation between steps. Given two time points n and m

the four equalities being consequence of the assumption that Xn has stationary increments. The smallest value k ¼ nm for which all four correlations become 0 characterises the memory of the process (Rangarajan and Ding, 2003). 3.3. Statistics on the step length Rn and on the tuning angles yn Distributions of the step length Rn and of the turning angle yn were computed using data from the 12 nematode trails. In accordance with the assumption that the process has stationary increments, then Rn and yn will have the same distribution, say R and y respectively, for all steps n ¼ 2; . . . ; 100. These common distributions were computed by averaging over the individual step length and turning angle distributions respectively for each of the 12 nematodes. The distribution of the step length (Fig. 1A) shows that some nematodes move faster than others and that in 8-s interval time they moved an average distance of 1 mm that represents one nematode body length. Also, unlike a random walk, the turning angle distribution (Fig. 1B) is peaked around the value 0, corresponding to a more forward movement. Correlation and cross-correlations between step lengths and turning angles were computed with time lags up to 10 steps. If we denote the time lag by k, k ¼ 1,y,10, then the correlation was computed between the step n+k and the step n, over all the step length or turning angle values corresponding to n ¼ 2,y,90 in all the 12 trails. We obtained a strong correlation (Po0.001) between different step lengths (Fig. 2A) with the correlation coefficient, rR, decreasing approximately linearly with time lag k. Also, the turning angles were strongly correlated (Po0.001) for k ¼ 1, 2, but then seemed to stabilise around ryE0.1 for larger k. The step length and turning angle distribution were also significantly correlated (Po0.05) although the crosscorrelation coefficient was considerably smaller (Fig. 3A and B).

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

step length-turning angle correlation

0.45

0.04

0.4

0.02

frequency

0.35 0.3

0

0.25

-0.02

0.2

0

1

2

3

4

5

6

7

8

9

10

9

10

-0.04

0.15

lag between steps

0.1

turning angle-step lenght correlation

0.05

0.12 0.1 0.08 0.06 0.04 0.02 0

0 0

0.5

1.5

1

2

2.5

mm 0.25

0 0.2 frequency

215

1

2

3

4 5 6 7 lag between steps

8

Fig. 3. Coefficient of cross-correlation rR,y between the step length and the previous turning angles (A) and coefficient of cross-correlation rR,y between the turning angle and the previous steps lengths (B) depending on the lag between steps.

0.15 0.1 0.05 0 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 radians

1

1.5

2

2.5

3

Fig. 1. Step length distribution (A), turning angle distribution (B).

step length correlation 1 0.8 0.6 0.4 0.2 0 1

2

3

4 5 6 7 lag between steps

8

9

10

8

9

10

turning angle correlation 0.2 0.15 0.1 0.05 0 1

2

3

4 5 6 7 lag between steps

Fig. 2. Coefficient of correlation rR between steps length (A) and coefficient of correlation ry between turning angles (B) depending on the lag between steps.

3.4. Building a correlated random walk on the square lattice From the statistical analysis above, we deduce that (Xn)nX0 is a stochastic process with memory. In particular, we obtained a strong correlation between steps meaning that each nematode had a preferential step length, which

produces a large individual variation amongst the population from very passive individuals (Fig. 4A) to those that move relatively fast (Fig. 4C and D). An account of this individual variation was the main challenge of the model. Nematode movement was simulated on a square lattice using a correlated random walk with conditional probabilities derived from data. Trail simulations were done on a square lattice of 0.5 mm unit corresponding to half of the nematode body length. On this lattice, the nematode can stay in the same position or can move left, right, forward, and backward, with respect to the previous directions and with a step length varying from 0 to 5 lattice units, with the maximum step length according to Fig. 1A. The direction of the movement is numbered as in Fig. 5, and the probability that the nematode moves in one of these directions and with a certain step length depends on the previous steps. These conditional probabilities were computed by projecting the 12 nematode trails on the square lattice as in Fig. 6. Therefore, we have associated a position j on the lattice with each position i of the real trail. The projected positions, the corresponding steps lengths, and the turning angles were computed in MatLab6.1 and used to build the matrices of the conditional probabilities. In this way, we have associated a correlated random walk ðX n Þn¼1;...;100 with the stochastic process ðX n Þn¼1;...;100 , on the square lattice that will be used to simulate the nematode trail. Based on the statistical results we know that ðX n Þn¼1;...;100 is a stochastic process with memory, for which a state Xn depends on the previous states up to a certain state Xmin{1,nk}, where k measures the memory of the process. Therefore if k ¼ 1, ðX n Þn¼1;...;100 is a Markovian process and for k ¼ 99, the present state of the system depends on its entire history. Suppose that for a fixed kA{2, 99}, we have built the correlated random walk up to the instant n1, then we know the positions. X¯ 1 ; . . . ; X¯ n1

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

216

Fig. 4. Examples of nematode trails on agar plate (solid line) and the corresponding projection on the square lattice (dot line).

3=forward

0=not move 2=left

4=right

1=back Fig. 5. Nematode directions on the lattice.

and we can deduce the corresponding step lengths ¯ 2; . . . ; R ¯ n1 and turning angles y¯ 3 ; . . . ; y¯ n1 . Then the R probability of the walker to visit next, the lattice point xn on the grid will depend on the previous positions and can be computed as follows: PðX¯ n ¼ xn jX¯ n1 ¼ xn1 ; . . . ; X¯ 1 ¼ x1 Þ ¼ PðX¯ n ¼ xn jX¯ n1 ¼ xn1 ; . . . ; X¯ nkþ1 ¼ xnk Þ ¯ n ¼ rn ; y¯ n ¼ an jR ¯ n ¼ rn1 ; y¯ n1 ¼ an1 ; . . . , ¼ PðR ¯ nkþ1 ¼ rnkþ1 Þ R ¯ n1 ¼ rn1 ; y¯ n1 ¼ an1 ; . . . , ¯ n ¼ rn jy¯ n ¼ an ; R ¼ PðR ¯ nkþ1 ¼ rnkþ1 Þ R ¯ n1 ¼ rn1 ; y¯ n1 ¼ an1 ; . . . , Pðy¯ n ¼ an jR ¯ nkþ1 ¼ rnkþ1 Þ. R

Fig. 6. Projection of the real trail on the square lattice: real trail (dot line), projected trail (solid line).

in Eq. (3) are computed using data from the projections of the real trails on the square lattice: ð3Þ

The ri ; i ¼ 1; . . . ; 100 above represent the step lengths on the grid and they can take values from 0 to 5 square units, according to Fig. 1A, and ai ; i ¼ 1; . . . ; 100 represents the turning angles on the grid, which take values from 0 to 4 as described in Fig. 5. The last two conditional probabilities

¯ n1 ¼ rn1 ; y¯ n1 ¼ an1 ; . . . ; R ¯ nkþ1 ¼ rnkþ1 Þ Pðy¯ n ¼ an jR ¯ n1 ¼ rn1 ; . . . ; R ¯ nkþ1 ¼ rnkþ1 Þ Pðy¯ n ¼ an ; R ¼ ¯ ¯ n1 ¼ rn1 ; yn1 ¼ an1 ; . . . ; R ¯ nkþ1 ¼ rnkþ1 Þ PðR ¼

Nðan ; rn1 ; an1 ; . . . ; ankþ2 ; rnkþ1 Þ , Nðrn1 ; an1 ; . . . ; ankþ2 ; rnkþ1 Þ

ð4Þ

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

where for any r2 ; . . . ; rkþ1 2 f0; . . . ; 5g and a3 ; . . . ; akþ1 2 f0; . . . ; 4g, Nðakþ1 ; rk ; . . . ; a3 ; r2 Þ represents the number of times the nematodes have taken these turning angles and step lengths values as a succession, and can be computed as follows: Nðakþ1 ; rk ; . . . ; a3 ; r2 Þ ¼

100 X X

217

six length probability ranges j X ¯ kþ1 ¼ ljy¯ kþ1 ¼ an ; . . . ; R ¯ 2 ¼ rnkþ1 Þ PðR PRR ðjÞ ¼ l¼0

j ¼ 0; . . . ; 5. PRy ð1Þ ¼ 0.

wðy¯ j ðoÞ¼akþ1 Þ

o2O j¼kþ1

wðR¯ j1 ðoÞ¼rk Þ . . . wðy¯ jkþ2 ðoÞ¼a3 Þ wðR¯ jkþ1 ðoÞ¼r2 Þ

ð5Þ

and

With the direction an ¼ i and the step length rn ¼ j chosen as above, the walker can move on the square lattice to the next state Xn. 4. Model application and results

Nðrn1 ; an1 ; . . . ; a3 ; r2 Þ ¼

4 X

4.1. Numerical simulations of nematode movement

Nðan ; rn1 ; . . . ; a3 ; r2 Þ.

an ¼0

Here w represents the indicator function and is equal to 1 if the quantity between parentheses is true and 0 otherwise. A similar formula can be written for the conditional probabilities of the step length ¯ n ¼ rn jy¯ n ¼ an ; . . . ; R ¯ nkþ1 ¼ rnkþ1 Þ PðR ¯ n ¼ rn ; y¯ n ¼ an ; . . . ; R ¯ nkþ1 ¼ rnkþ1 Þ PðR ¼ ¯ ¯ n1 ¼ rn1 ; . . . ; R ¯ nkþ1 ¼ rnkþ1 Þ Pðyn ¼ an ; R ¼

Nðrn ; an ; . . . ; ankþ2 ; rnkþ1 Þ , Nðan ; rn1 ; . . . ; ankþ2 ; rnkþ1 Þ

ð6Þ

with the corresponding Nðrn ; an ; . . . ; ankþ2 ; rnkþ1 Þ. Recalling the fact that the process (Xt)tX0 and implicitly the random walk ðX¯ n Þn¼1;...;100 , have stationary increments, we note that the conditional probabilities above do not depend on the instant n but only on the previous k1 steps. As a consequence we have: ¯ n1 ¼ rn1 ; y¯ n1 ¼ an1 ; . . . ; R ¯ nkþ1 ¼ rnkþ1 Þ Pðy¯ n ¼ an jR ¯ k ¼ rn1 ; . . . ; R ¯ 2 ¼ rnkþ1 Þ, ¼ Pðy¯ kþ1 ¼ an jR and a similar equality for the conditional probability of the step length in formula (6). Having calculated the conditional probabilities, we can determine the step n of the correlated random walk. First the direction is generated by producing five angular probability ranges: PRy ðiÞ ¼

i X

¯ k ¼ rn1 ; y¯ k ¼ an1 ; . . . ; R ¯ 2 ¼ rnkþ1 Þ, Pðy¯ kþ1 ¼ ljR

l¼0

i ¼ 0; . . . ; 4.

It is obvious than PRy(4) ¼ 1 and let us set PRy(1) ¼ 0. We produce then a random number between 0 and 1, and determine in which of the intervals ½PRy ði  1Þ; PRyþ1 ðiÞ; i ¼ 0; . . . ; 4 this number falls. The number i selected gives the value of the next direction an. The larger that one of the five intervals is, the more likely is that it will be selected (Anderson et al., 1997b). The same method is used to determine the next step length n, by computing the

Following the model description above, numerical simulations of nematode trails on the square lattice were produced (Fig. 7). However, due to computational constraints and also to the limited number of real trails used to compute the conditional probabilities of the correlated random walk, it was technically impossible to take into account the entire memory of the movement. Based on the statistical results and on many preliminary numerical simulations, the longest memory we could account with the existing data was given by k ¼ 5 in formula (3). Therefore, we deduced from formula (4) that the next turning angle depends on the previous three turning angles and on the previous four step lengths. Also in formula (6) the next step length depends on the four turning angles and the four previous step lengths. Applying the method presented in the previous paragraph, we produced numerical simulations of trails (Fig. 7) in which, due to the conditional probabilities derived from data we were able to preserve the large individual variation, the simulations describing movement of passive individuals (Fig. 7A) but also very active individuals (Fig. 7C and D). 4.2. Quantifying the nematode dispersal in two dimensions This was achieved by computing PðX n ¼ xn jX 1 ¼ x1 Þ; i.e., the probability that a nematode, starting from a point x1 visits a certain point xn of the square lattice at instant nA1,y,100. Theoretically, this can be done using formula (3) together with a recurrent procedure. In practical terms, the limitations of computer storage make this prohibitive due to the large numbers of possible combinations of step length over the time lag. Therefore, it is easier to compute this probability by simulating a large number of realisations of the nematode trail, and counting the number of times that each point of the lattice is visited at instant n. In our case, we have simulated 10-million nematode trails on a 200  200 square lattice, with the starting point x1 ¼ (100, 100) and for 100-time steps corresponding to 15-min real time. The number of trail realisations used, seemed to be sufficient, since the simulated PDF (Fig. 8) satisfies the radial property with respect to the starting point x1, i.e., p¯ t ðxÞ ¼ functiont ðjjx  x1 jjÞ, for any lattice

ARTICLE IN PRESS 218

S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

Fig. 7. Numerical simulations of nematode trail on the square lattice.

Fig. 8. Numerical simulations of nematodes dispersal corresponding to 5, 10, and 15 min real time. White colour means high density of nematodes and dark colour low density.

point x of the square lattice. Further, the simulated PDF profile was computed using the formula: profilet ðrÞ ¼

X 1 p¯ ðxÞ, cardfx; rpjjx  x1 jjor þ 1g rpjjxx jjorþ1 t 1

where r ¼ 1,y,99 lattice units. We then needed to identify the type of diffusion to which this profile (Fig. 9) belongs. As mentioned in Section 2.1, the MSD defined by formula (1) is a very useful tool. Fig. 10 compares the MSD of the real movement computed over the 12 real trails with the MSD obtained from the simulated movement. It is clear that the correspondence between the two functions is poor. This lack of correspondence is due to the fact that the simulated trails use patterns found in the real trails averaged over the whole trail. Therefore, to compute the appropriate MSD for the real movement, it is necessary to reconsider the real trails as follows. Each trail is split into

six smaller overlapping trails of length equal to 50 steps. The first of the smaller trails covers the interval from step 1 to 50. The second covers the interval between steps 11 and 61, and the third from step 21 to 71, and so on until the final interval, which covers the steps from 51 to 100. Thus we obtained 12  6 ¼ 72 different samples of the real nematode trails of only 50 time steps that were used to compute a new MSD associated to the real movement (Fig. 11). Comparing observed MSD with that simulated, we found that the two MSDs are in good agreement and both correspond to a power law with exponent a ¼ 0.84. A detailed analysis of the whole MSD derived from the simulations (Fig. 12, corresponding to 100 time steps) suggested that it can be described by Eq. (A.8), corresponding to a TFD with finite memory. Fitting (A.8) to the data, we find d ¼ 33 and the MSD follows a power-law pattern with exponent a ¼ 0.84 for tp33 time steps and then becomes linear. The full-fitted equation is of its form

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

3

350

x 10-3

MSD (real movement) 5min 10min 15min

MSD (simulation)

300

Power-law (MSD)

square lattice units

2.5 2 PDF(x)

219

1.5 1

250 200 150 100 y = 11.61x0.84 R2 = 0.963 (real movement) R2 = 0.998 (simulation)

50 0.5

0

20

40 60 x-distance from origin

80

51

46

41

36

31

26

21

16

100

Fig. 9. Profile of the probability density function p¯ t ðxÞ obtained from the simulated trails corresponding to 5, 10, and 15 min real time.

time steps Fig. 11. The new mean square displacement (MSD) associated to the real movement, the MSD of simulated movement over 50 time steps and the power-law regression associated.

700

700 MSD (real movement)

100

y = 6.03x 2 R = 0.994

91

82

0 73

Fig. 10. Mean square displacement of the real movement and of the simulated movement.

y = 11.61x0.84 R2 = 0.995

64

time steps

200

55

99

92

85

78

71

64

57

50

43

36

29

22

15

8

0

300

46

100

37

200

400

28

300

500

1

400

19

500

1

MSD (simulation) Power-law (MSD) Linear (MSD)

600

MSD (simulation)

10

600

square lattice units

square lattice units

11

6

1

0 0

time steps Fig. 12. Mean square displacement of the simulated movement over 99 time steps, the power-law regression, and the linear regression associated.

11:61 minft; 33g0:84 þ 5:57 maxft  33; 0gÞ; R2 ¼ 0:999. For comparison, Fig. 12 shows the MSD profiles for a powerlaw regression (11.61t0.84, R2 ¼ 0.995) and a linear regression (6.03t, R2 ¼ 0.994). Applying formula (A.8), we can derive the diffusion coefficient KTFD,0.84 ¼ 11.61  G(1.84)/4 ¼ 2.74 corresponding to the TFD and from formula (A.11) the diffusion coefficient KFBM,0.84 ¼ 11.61/4 ¼ 2.90. Formula (A.11) can also be applied for a ¼ 1 to compute the diffusion coefficient KND ¼ 6.03/4 ¼ 1.51 corresponding to the normal diffusion (Figs. 11 and 12). To test these different theoretical models (TFD), the TFD with finite memory (TFD-FM), the FBM and the normal diffusion (ND)) we computed the probability over time, of the nematode to be at distance r from the starting point, denoted by Pt(r). Due to the radial property of the PDF corresponding to all these types of diffusion, Pt(r) can be computed using the following formula: Pt ðrÞ ¼ 2pr  PDFt ðxÞ;

r ¼ jjxjj,

where the PDF for the TFD, sub-diffusive case and TFD with finite memory, is the solution of the fractional equations (A.1) and (A.6), respectively, presented in Appendix A that can be solved numerically by implicit finite-difference approximations (Appendix B). For the PDF of the normal diffusion and the FBM, formula (A.10) in Appendix A was used. All these theoretical probabilities were compared in Fig. 13A and B with the probability derived from the simulated correlated random walk model. We deduce that the simulated movement is governed by the TFD with finite memory, as the corresponding probability fits the simulated profile best. 5. Summary and conclusions In this paper, we have developed a theoretical model to better understand the nematode dispersal in homogeneous environments. It is an individual-based model that uses

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

220

t=5min

t=15min

t=10min 0.06

0.06 simulation TFD TFD-FM

probability

0.06

TFD TFD-FM

0.04

0.04

0.04 0.02

0.02

0.02 0

0 0

20

40

60

80

100

0 0

20

40 60 80 distance from origin

t=5min

100

20

FBM ND

40

60

80

100

t=15min 0.06

simulation

0.06

0

t=10min 0.06

probability

simulation TFD TFD-FM

simulation

0.04

simulation

simulation

FBM ND

FBM ND

0.04

0.04 0.02

0.02

0.02

0

0 0

20

40

60

80

100

0 0

20

40 60 80 distance from origin

100

0

20

40

60

80

100

Fig. 13. (A, B) Comparison of different probabilities Pt(r), obtained from simulated trails (red), time-fractional dynamics (TFD, blue), time-fractional dynamics with finite memory (TFD-FM, black), fractional Brownian motion (FBM, green) and normal diffusion (ND, dark-green) corresponding to 5, 10, and 15 min real time.

digitised data on nematode trails to predict the dispersal in time and space of the whole population. The main challenge of the model was a large individual variation observed amongst the population from very passive individuals to individuals that move very fast, and also the fact that the movement was characterised by strong correlations between steps, each nematode having a preferential step length and correlated turning angles. Statistical analysis of the data inspired us to simulate nematode movement on the square lattice by 2d correlated random walks whose conditional probabilities are derived from the data. The most important finding is that the stochastic process associated with these correlated random walks is one with memory governed by the time-fractional sub-diffusion equation with finite memory (A.6). This finite memory characterisation of the simulated trails is a consequence of the construction of the correlated random walk on the square lattice, where due to technical limitations it was not possible to consider the whole memory of nematode trail, and where we assumed that a given position of the system would depend only on five previous positions. In reality, the observed trails showed long-range correlations between steps, especially between the steps length. Therefore, we believe the population dispersion is characterised by the time-fractional Eq. (A.1), where the whole memory of the system is taken into account. The parameters involved in these equations have been derived from the MSD of simulated tails, which is

related to the MSD associated to the real trails (Fig. 11). In many situations, the MSD has been shown to be an important tool for identifying the type of diffusion associated with a certain stochastic process (Metzler and Klafter, 2000, 2004; Wiktorsson et al., 2004). In our case, the MSD of the simulated trails follows a power-law pattern, and further comparison of the associated probability shows that the stochastic process is governed by a time-fractional equation of sub-diffusion type. This subdiffusive motion can find explanation in the individual movement patterns for which two main behaviours have been identified: ‘‘looping’’ and ‘‘reversals’’ (often fluctuating the backward and forward movement, Anderson et al., 1997a, b). For the nematodes that adopt the ‘‘reversals’’, these can be considered as no-motion periods or ‘‘imaginary traps’’ since in this case the individuals do not move very far from the point they were released. The same may happen to the nematodes that adopt a preferential turning angle, producing loops in the movement and for which there is a large probability to visit the same point for many times. Following this idea, the reason that the process is characterised by ao1, and consequently is a sub-diffusion process might be related to the fact that our model relies on data of trails from 12 different nematodes as being representative from the whole population. Due to the large range of individual behaviour within the population, the model result obtained in this paper might be a consequence of the fact that among the 12 nematodes the

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

proportion of those that like to stay close to the application point is greater than the proportion of nematodes that move away from the origin. We believe that with a different population of nematodes (in which the proportion of nematode that persists in preserving the same direction and consequently move away from the origin is greater than the proportion of nematodes that does not move very far from the starting point), it may be possible to obtain superdiffusive behaviour (1oao2) (Appendix A). The important thing about these two different types of diffusion is that both can explain the behaviour of organisms in which the movement presents a long-range correlation between steps. This long-range correlation is the reason of the anomalous nature of nematode movement and is in agreement with the fact that the associated stochastic process is characterised by strong memory effects on the level of the PDF. In the last decades, fractional dispersion has been shown to be a promising alternative for modelling anomalous diffusion in many biological systems from bacterial motion (Levandowsky et al., 1997) and soil larvae (Zhang et al., 2007) to other animals such as the spider monkey (Ramos-Fernandez et al., 2004), jackals (Atkinson et al., 2002), or the famed flight of an albatross (Viswanathan et al., 1996). Our study proves that nematode movement also generates anomalous diffusion that can be modelled by TFD, and which can also explain some previous results obtained in soil related to nematode dispersal that cannot be described by normal diffusion. The main advantage of the method is that it relies on real data, and therefore can be applied to other organisms for which movement can be recorded. Acknowledgements We would like to thank Dr. Xiaoxian Zhang for his useful suggestions and discussions related to time-fractional dynamics and their numerical resolution. This work was funded by the Biotechnology and Biological Sciences Research Council (Grant No. BBSB01065). Appendix A. Some examples on anomalous diffusion In this paragraph, we briefly introduce a few types of anomalous diffusion that have been developed in the last decade to model the dynamic of different biological systems. In the anomalous regime, an important indicator is the deviation of the MSD from the linear dependence MSD(Xt)EKt on time. In the next, we like to pay particular attention to those for which the MSD follows the power-law pattern in formula (2). A.1. Time-fractional dynamics The most used approach to describe fractional diffusion is given by the continuous time random walk (CTRW)

221

model (Metzler and Klafter, 2000, 2004; Metzler and Nonnenmacher, 2002). The CTRW model is based on the idea that the length of a given jump, as well as the waiting time elapsed between two successive jumps are drawn from a joint PDF c(x, t) called the jump PDF. If the jump length and the waiting time are independent, one finds the decoupled form c(x, t) ¼ l(x)w(t); if both are coupled, i.e., c(x, t) ¼ l(x)p(t/x) ¼ w(t)p(x/t) then a jump of a certain length involves a time cost. In the decoupled version, different types of CTRW processes can be R1 categorised by the mean waiting time T ¼ 0 twðtÞ dt and R þ1 the mean square of the jump length S2 ¼ 1 x2 lðxÞ dx. If both mean waiting times and mean square of the jump length, T and S2 are finite, then the long time limit corresponds to Brownian motion and infinite value of one of these moments results in anomalous behaviour. If ToN and S2 ¼ N (e.g., PDF of jump length has a heavy powerlaw tail) the CTRW model yields to the super-diffusive Levy flight described by a diffusion equation with fractional space derivatives (Metzler and Klafter, 2000, 2004). Due to the finiteness of T, this process is of Markovian nature while the fact that the jump length distribution does not posses a moment of second order causes the divergence of the MSD. If T ¼ N (e.g., the PDF of the waiting time has a power-law tail) and S2oN, the CTRW model leads to sub-diffusion and a detailed analysis of this process is presented below. In the last case when both T and S2 are diverging, corresponding to long jumps and long rests, depending on this competition we can get sub- and super-diffusive scalings and the diffusion equation associated will have both fractional time (memory) and fractional space derivatives. Here again the MSD diverges (Metzler and Klafter, 2000, 2004). Another process that can be associated to the CTRW model is the Levy walk that corresponds to the coupled version of the jump PDF c(x, t) ¼ l(x)p(t/x) and that has received a particular attention in the case when the conditional probability pðt=xÞ ¼ 12dðjxj  vtg Þ is introduced with a generalised velocity v, which penalises the long jumps (corresponding to S2 ¼ N) such that the overall process is super-diffusive and has a finite MSD described by the power-law pattern (2) with 1oao2. In what follows, particular attention will be given to the sub-diffusion process for which it has been shown that the associated PDF satisfies the following fractional diffusion equation (Metzler and Klafter, 2000): qt pt ðxÞ ¼ 0 Dt1a K a Dx pt ðxÞ with 0oao1,

(A.1)

also known as TFD equation. Although it cannot be derived from the CTRW model (Barkai, 2002), Eq. (A.1) can be formally extended to the super-diffusive case for 1oao2 (Metzler and Nonnenmacher, 2002) in which case the process is governed by a fractional wave equation: q2t pt ðxÞ ¼ 0 D2a K a Dx pt ðxÞ with t

1oao2,

(A.2)

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

222

where for any 0obo1, 0 Dbt is the Riemann–Liouville operator (Samko et al., 1993) defined by b 0 Dt f ðt; xÞ ¼

1 q Gð1  bÞ qt

Z

t 0

f ðs; xÞ ðt  sÞb

ds.

(A.3)

According to formula (A.3), the integro-differential nature of Riemann–Liouville fractional operator 0 D1a for t 0oao1 (or 0 D2a for 1oao2), with the internal kernel t Kðt; sÞ ¼ ðt  sÞa1 (Kðt; sÞ ¼ ðt  sÞa2 , respectively), ensures the non-Markovian nature of the associated stochastic process ðX t Þt40 , for which the MSD can be computed explicitly: 4K a MSDðX t Þ ¼ ta for 0oao1 (A.4) Gð1 þ aÞ

denote the length of this memory by d with dA(0, T). Then using formulae (A.1)–(A.3) we deduce the corresponding diffusion equations for 0oao1: Z Ka q t Dx ps ðxÞ ds qt pt ðxÞ ¼ GðaÞ qt 0 minft  s; dg1a 8 R K a q t Dx ps ðxÞ > > GðaÞ qt 0 ðtsÞ1a ds for 0ptpd; > > < R Dx ps ðxÞ Ka Ka q t D p ðxÞ þ GðaÞ ¼ qt td ðtsÞ1a ds GðaÞd1a x td > > > > : ðA:6Þ for dotoT; and for 1oao2: Z Ka q t Dx ps ðxÞ ds Gða  1Þ qt 0 minft  s; dg2a 8 R Ka q t Dx ps ðxÞ > > Gða1Þ qt 0 ðtsÞ2a ds for 0ptpd; > > < R Dx ps ðxÞ Ka Ka q t ¼ ds 2a Dx ptd ðxÞ þ Gða1Þ qt td Gða1Þd ðtsÞ2a > > > > : for dotoT; ðA:7Þ

q2t pt ðxÞ ¼

and MSDðX t Þ ¼

4K a a t GðaÞ

for

1oao2,

(A.5)

where Ka is the generalised diffusion coefficient that appears in formulae (A.1) and (A.2). Let us give the proof of formula (A.4). First, from the definition of the MSD, formula (1), can be rewritten using the PDF as Z MSDðX t Þ ¼ jjxjjpt ðxÞ dx. R

2

Therefore, by multiplying formula (A.1) by ||x||2 and integrating over R2 we obtain d MSDðX t Þ dt Z  Z Ka d t 1 2 ¼ jjxjj Dx pt ðxÞdx ds. GðaÞ dt 0 ðt  sÞ1a R2 The second integral in the above formula can be easily computed integrating by parts and using the fact that the PDF, pt has an exponential decay (Metzler and Klafter, 2000) when jjxjj ! 1. Therefore,  Z Z q jjxjj2 Dx pt ðxÞ dx ¼ ðx21 þ x22 Þ p ðx1 ; x2 Þ qx21 t R2 R2  Z q þ 2 pt ðx1 ; x2 Þ dx1 dx2 ¼ 4 pt ðxÞ dx ¼ 4. qx2 R2 Combining the last two formulae, we obtained that ðd=dtÞMSDðX t Þ ¼ 4K a ta1 =GðaÞ, from which formula (A.4) follows directly. This proof can be easily adapted for formula (A.5). In both cases described by Eqs. (A.1) and (A.2), the associated process (Xt)t40 becomes a stochastic process with memory for each a given state of the system depends on the whole history. Numerical solutions of Eqs. (A.1) and (A.2) are presented in Appendix B. In what follows, we would like to generalise this concept by introducing the stochastic process with limited memory, for which the memory of the system is finite in time. Let us

that are obtained by truncating the kernel K(t, s), so that at instant t the system depends on the previous states up to the instant td. For this process with finite memory, the associated MSD is given by MSDðX t Þ ¼

4K a 4K a minft; dga þ Gð1 þ aÞ GðaÞd1a  maxft  d; 0g

for

0oao1

ðA:8Þ

and MSDðX t Þ ¼

4K a 4K a minft; dga þ GðaÞ Gða  1Þd2a  maxft  d; 0g2

for

1oao2.

ðA:9Þ

The proof of (A.8) and (A.9) can be easily obtained following the same method as for the proof of (A.4) presented previously. A.2. Fractional Brownian motion Unlike the TFD, the 2d FBM (Xt)t40 is a stochastic process for which the corresponding diffusion equation cannot be expressed in terms of partial differential equations (it is possible only for the one-dimensional case, Baudoin and Coutin, 2007). It is defined as a stationary stochastic process, for which the increments Y t ¼ X tþt  X t ; t40 follow the Gaussian distribution Nð0; 2K a ta Þ with 0oao2, where Ka is the corresponding diffusion coefficient (Revuz and Yor, 1999; Enriquez, 2004; Metzler and Klafter, 2004). Therefore, the PDF in 2d is given by pt ðxÞ ¼

1 expðjjxjj2 =ð4K a ta ÞÞ. 4pK a ta

(A.10)

If 0oao1 then the behaviour of FBM is antipersistent, and persistent for 1oao2 (Mandelbrot, 1982). In both cases, the increments are strongly correlated and so (Xt)t40

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

is a stochastic process with memory for which the MSD is given by MSDðX t Þ ¼ 4K a ta

for

0oao2.

(A.11)

Appendix B. Numerical solutions of the time-fractional equation The sub-diffusion problem (A.1) with a boundary condition can be solved numerically using implicit Euler finite-difference approximations (Raviart and Thomas, 1983). In our case, the following two-dimensional fractional diffusion equation with homogeneous Dirichlet boundary condition is considered: qt u ¼ 0 D1a K a Dx u; ðt; xÞ 2 ð0; TÞ  A; t uðt; xÞ ¼ 0; ðt; xÞ 2 ð0; TÞ  qA; uð0Þ ¼ dx1 ;

(B.1)

x 2 A;

where T corresponds to 15 min real time, A ¼ (0,200)  (0,200) is the same square that we have used for the simulation and x1 ¼ (100,100) is the stating point. The Riemann–Liouville fractional operator 0 Dt1a is the same defined in formula (A.3). Let us take m and n two positive integers and let us discretise the space using the following square lattice: ðxi;j Þ0pi;jp2m ; xi;j ¼ ðiDx; jDxÞ ¼ ð100i=m; 100j=mÞ, and the time by ðtk Þ0pkpn ; tk ¼ k Dt ¼ kT=n. Then, the lefthand side can be approximated by qt uðtk ; xi;j Þ 

uðtk ; xi;j Þ  uðtk1 ; xi;j Þ . Dt

(B.2)

As for the right side we have: Z

tk

Dx uðs; xi;j Þ ds ðtk  sÞ1a Z tk  Z tk1 Dx uðs; xi;j Þ Dx uðs; xi;j Þ 1 ¼ ds  ds GðaÞDt 0 ðtk  sÞ1a ðtk1  sÞ1a 0 ! Z k k1 Z tl tl X X 1 Dx uðs; xi;j Þ Dx uðs; xi;j Þ ¼ ds  ds 1a GðaÞDt l¼1 tl1 ðtk  sÞ1a l¼1 tl1 ðtk1  sÞ ! k k1 X Dx uðtl ; xi;j Þ X Dx uðtl ; xi;j Þ 1   1a 1a GðaÞ l¼1 ðtk  tl Þ ðt l¼1 k1  tl Þ " k1 X Dx uðtk ; xi;j Þ 1 þ Dx uðtl ; xi;j Þ ¼ 1a GðaÞ ðtk  tk1 Þ l¼1  # 1 1   ðtk  tl Þ1a ðtk1  tl Þ1a " k1 X 1 Dx uðtk ; xi;j Þ þ Dx uðtl ; xi;j Þ  1a GðaÞðDtÞ l¼1  # 1 1 .   ðB:3Þ ðk  lÞ1a ðk  1  lÞ1a

1a Dx uðtk ; xi;j Þ ¼ 0 Dt

1 q GðaÞ qt

0

223

So, gathering formulas (B.1)–(B.3) we obtain uðtk ; xi;j Þ Dx uðtk ; xi;j Þ uðtk1 ; xi;j Þ Ka  Ka þ ¼ 1a Dt Dt GðaÞðDtÞ GðaÞðDtÞ1a   k1 X 1 1  Dx uðtl ; xi;j Þ  , ðB:4Þ ðk  lÞ1a ðk  1  lÞ1a l¼1 where the Laplacian can be approximated also by finite differences: Dx uðtk ; xi;j Þ 

uðtk ; xiþ1;j Þ þ uðtk ; xi1;j Þ þ uðtk ; xi;jþ1 Þ þ uðtk ; xi;j1 Þ  4uðtk ; xi;j Þ ðDxÞ2

for 1pi, jp2m1 and 1pkpn. The solution is equal to 0 on the boundary, for any tk. As for the initial data, we shall take it equal to 1/(Dx)2 in the starting point x1 ¼ xm,m and 0 elsewhere on the other points of the grid. By solving system (B.3) we get an approximation for the solution of problem (B.1). From the analytical point of view, the solution of problem (A.1) in infinite domain can be expressed explicitly in terms of Fox functions (Samko et al., 1993; Metzler and Klafter, 2000), which helps to derive the asymptotical behaviour: rffiffiffiffiffiffiffiffiffiffiffi ð1aÞ=ð2aÞ   1 1 2 jjxjj ð1aÞ=ð2aÞ ffiffiffiffiffiffiffiffiffi ffi p pt ðxÞ  4pK a ta 2  a a K a ta " #   2  a aða=ð2aÞÞ jjxjj ð2=ð2aÞÞ pffiffiffiffiffiffiffiffiffiaffi  exp  ðB:5Þ 2 2 K at pffiffiffiffiffiffiffiffiffiffi valid for jjxjj4 K a ta . The homogeneous Dirichlet boundary condition does not affect the solution of problem (B.1) much with respect to the Cauchy solution of problem (A.1), which was considered on whole space R2, since thanks to (B.5), the Cauchy solution became very close to 0 for x far enough from the starting point. In Fig. 14, the profile of the numerical solution of problem (B.1) is compared to the asymptotic profile described by formula (B.5). The fractional equation with finite memory (A.6) can be numerically solved using the same procedure. The only difference appears in formula (B.4) when summing over the time previous to the instant tk. In this case, we have to start summing from kk0 to k1 as soon as k4k0+1, where k0 corresponds to d, the memory of the process that we take into account in the integro-differential operator in formula (A.6). As for Eqs. (A.2) and (A.7) corresponding to 1oao2, the same numerical scheme can be employed with the difference that we have the second derivate in time and therefore we need two initial conditions, the first one is u(0) that can be choose as previously and the second u(t1) needs a particular attention, as in order to preserve the mass it has to be positive, radial with respect to the starting point x1, and of mass equal to 1. As stated in Metzler and Klafter (2004), the analytical solution of Eq. (A.2) also satisfies the

ARTICLE IN PRESS S. Hapca et al. / Journal of Theoretical Biology 248 (2007) 212–224

224

x 10-3

2

x 10-3

1.5

x 10-3

3 PDF(x)

1.5 2

1

1

1

0.5

0.5

0

0 0

50

100

0 0

50 100 distance from origin

0

50

100

Fig. 14. Comparison between the numerical solution of problem (B.1) (solid line) and the asymptotical approximation (B.5) (dot line) for a ¼ 0.84 and KTFD,0.84 ¼ 2.90, corresponding to 5, 10, and 15 min real time.

asymptotical behaviour described by formula (B.5), this time with 1oao2. References Anderson, A.R.A., Sleeman, B.A., Young, I.M., Griffiths, B.S., 1997a. Nematode movement along a chemical gradient in a structurally heterogeneous environment. 1. Experiment. Fundam. Appl. Nematol. 20, 157–163. Anderson, A.R.A., Sleeman, B.A., Young, I.M., Griffiths, B.S., 1997b. Nematode movement along a chemical gradient in a structurally heterogeneous environment. 1. Theory. Fundam. Appl. Nematol. 20, 165–172. Atkinson, R.P.D., Rhodes, C.J., Macdonald, D.W., Anderson, R.M., 2002. Scale-free dynamics in the movement patterns of jackals. OIKOS 98, 134–140. Barkai, E., 2002. CTRW pathways to the fractional diffusion equation. Chem. Phys., 13–27. Baudoin, F., Coutin, L., 2007. Operators associated with a stochastic differential equation driven by fractional Brownian motions. Stochastic Process. Appl. 117 (5), 550–574. Bell, W.J., 1991. Searching Behaviour: The Behavioural Ecology of Finding Resources. Chapman & Hall, London, 358pp. Bengtsson, G., Nilsson, E., Ryden, T., Wiktorsson, M., 2004. Irregular walks and loops combines in small-scale movement of a soil insect: implications for dispersal biology. J. Theor. Biol. 231 (2), 299–306. Banks, H.T., Kareiva, P.M., Murphy, K., 1987. Parameter-estimation techniques for interaction and redistribution models, a predator–prey example. Oecologica 74, 356–362. Blanche, S., Casas, J., Bigler, F., Janssen-van Bergeijk, K.E., 1996. An individual-based model of Trichogramma foraging behaviour: parameter estimation for single female. J. Appl. Ecol. 33, 425–434. Bovet, P., Benhamou, S., 1988. Spatial analysis of animal movement using a correlated random walk model. J. Theor. Biol. 131, 419–433. Conradt, L., Bodsworth, E.J., Roper, T.J., Thomas, C.D., 2000. Nonrandom dispersal in the butterfly Maniola jurtina: implications for metapopulation models. Proc. R. Soc. London, Ser. B—Biol. Sci. 267, 1505–1510. Croll, N.A., Sukhdeo, M.V.K., 1981. Hierarchies in nematode behaviour. In: Zuckermann, B.M., Rohne, R.A. (Eds.), Plant Parasitic Nematodes, vol. 3. Academic Press, New York, pp. 227–251. Enriquez, N., 2004. A simple construction of the fractional Brownian motion. Stochastic Process. Appl. 109, 203–223. Kareiva, P.M., Shigesada, N., 1983. Analyzing insect movement as a correlated random walk. Oecologia 56, 234–238. Levandowsky, M., White, B.S., Schuster, F.L., 1997. Random movements of soil amebas. Acta Protozool. 36, 237–248. Levin, S.A., 1987. Ecological and evolutionary aspects of dispersal. In: Hallam, T., Levin, S.A. (Eds.), Mathematical Topics in Population

Biology, Morphogenesis and Neurosciences. Lecture Notes in Biomathematics, vol. 71. Springer, Berlin, pp. 80–87. Mandelbrot, B.B., 1982. The Fractal Geometry of Nature. W.H. Freeman et Co., New York, 468pp. Metzler, R., Klafter, J., 2000. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep. 339, 1–77. Metzler, R., Klafter, J., 2004. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J. Phys. A: Math. Gen. 37, 161–208. Metzler, R., Nonnenmacher, T.F., 2002. Space- and time-fractional diffusion and wave equations, fractional Fokker–Plank equations, and physical motivation. Chem. Phys. 284, 67–90. Preisler, H.K., Ager, A.A., Johnson, B.K., Kie, J.G., 2004. Modelling animal movement using stochastic differential equations. Environmetrics 15, 643–657. Ramos-Fernandez, G., Mateos, J.L., Miramontes, O., Cocho, G., Larralde, H., Ayala-Orozco, B., 2004. Levy walk patterns in the foraging movements of spider monkeys (Ateles geoffroyi). Behav. Ecol. Sociobiol. 55, 223–230. Rangarajan, G., Ding, M. (Eds.), 2003. Processes with Long Range Correlations: Theory and Applications. Springer Lecture Notes in Physics 621. Springer, New York, 392pp. Raviart, P.A., Thomas, J.M., 1983. Introduction a l’analyse nume´rique des e´quations aux de´rive´es partielles. Masson, Paris, 224pp. Revuz, D., Yor, M., 1999. Continuous martingales and Brownian motion, third ed. Springer, Berlin, 602pp. Samko, S.G., Kilbas, A.A., Marichev, O.I., 1993. Fractional Integral and Derivatives. Theory and Applications. Gordon and Breach, New York, 976pp. Shapiro, D.I., Glazer, I., 1996. Comparison of entomopathogenic nematode dispersal from infected hosts versus aqueous suspension. Environ. Entomol. 25, 1455–1461. Stamps, J.A., Buechner, M., Krishnan, V.V., 1987. The effects of edge permeability and habitat geometry on immigration from patches of habitat. Am. Nat. 129, 533–552. Viswanathan, G.M., Afanasyev, V., Buldyrev, S.V., Murphy, E.J., Prince, P.A., Stanley, H.E., 1996. Levy flight search patterns of wandering albatrosses. Nature 381, 413–415. Wiktorsson, M., Ryden, T., Nilsson, E., Bengtsson, G., 2004. Modelling the movement of a soil insect. J. Theor. Biol. 231, 497–513. Wilson, M.J., Glen, D.M., George, S.K., 1993. The Rhabditid nematode Phasmarhabditis hermaphrodita as a potential biological control agent for slugs. Biocontrol Sci. Technol. 3, 513–521. Wu, H., Li, B.-L., Springer, T.A., Neill, W.H., 2000. Modelling animal movement as a persistent random walk in two dimensions: expected magnitude of net displacement. Ecol. Modelling 132, 115–124. Zhang, X., Johnson, S.N., Crawford, J.W., Gregory, P.J., Young, I.M., 2007. A general random walk model for leptokurtic distribution of organisms movement: theory and application. Ecol. Modelling 200, 79–88.