## Abstract

The characterization of the dispersal of populations of non-identical individuals is relevant to most ecological and epidemiological processes. In practice, the movement is quantified by observing relatively few individuals, and averaging to estimate the rate of dispersal of the population as a whole. Here, we show that this can lead to serious errors in the predicted movement of the population if the individuals disperse at different rates. We develop a stochastic model for the diffusion of heterogeneous populations, inspired by the movement of the parasitic nematode *Phasmarhabditis hermaphrodita*. Direct observations of this nematode in homogeneous and heterogeneous environments reveal a large variation in individual behaviour within the population as reflected initially in the speed of the movement. Further statistical analysis shows that the movement is characterized by temporal correlations and in a heterogeneously structured environment the correlations that occur are of shorter range compared with those in a homogeneous environment. Therefore, by using the first-order correlated random walk techniques, we derive an effective diffusion coefficient for each individual, and show that there is a significant variation in this parameter among the population that follows a gamma distribution. Based on these findings, we build a new dispersal model in which we maintain the classical assumption that individual movement can be described by normal diffusion, but due to the variability in individual dispersal rates, the diffusion coefficient is not constant at the population level and follows a continuous distribution. The conclusions and methodology presented are relevant to any heterogeneous population of individuals with widely different diffusion rates.

## 1. Introduction

The study of organism movement and dispersal has become a key element for understanding a series of ecological questions related to the spatio-temporal dynamics of populations (Levin 1974; Kareiva 1990), foraging strategies (Bell 1991) and the community dynamics of interacting species (Banks *et al*. 1987; Taylor 1990, 1991). Therefore, a description and analysis of the process of movement is necessary to help us adequately understand the consequences of movement behaviour. One of the approaches to create such a link is to analyse the movement of individuals and use this description to directly infer the characteristics of the spatial spread of the whole population (Blanche *et al*. 1996; Skalski & Gilliam 2000; Wiktorsson *et al*. 2004). Unfortunately, in many cases, detailed observations of individuals and their behaviour are very difficult to obtain. In these situations, ecologists use classical diffusion models to quantify organisms' dispersal in different environments.

A common starting point of these models is to assume that the movement of the organism is random (Brownian motion; Croll & Sukhdeo 1981; Anderson *et al*. 1997; Turchin 1998; Okubo & Levin 2002). However, experimental observations over the past two decades have shown that the spatial spread of some organisms follows the leptokurtic distribution, where the probability per unit time of individuals moving short and long distances is higher compared with the Gaussian distribution predicted by the Brownian motion (Kot *et al*. 1996; Cain *et al*. 1998; Nathan 2001). Different assumptions have been made to explain such behaviour. One of the most frequent is that the individuals do not move according to Brownian motion and therefore cannot be described in terms of classical diffusion. Direct observations suggested that many entities, including animals, exhibit correlated movement, i.e. subsequent direction of the movement depends on the previous directions (Hall 1977; Kareiva & Shigesada 1983; Bergman *et al*. 2000). As a consequence, the turning angle between successive steps in a path is not uniformly distributed in the interval (−*π*,*π*) and characteristics of this distribution are very important when modelling the movement patterns by correlated random walks (CRWs; Kareiva & Shigesada 1983; Bovet & Benhamou 1988; Benhamou 2004, 2006; Codling *et al*. 2008). When there is a bias in either direction, the CRW can be used to model the loops observed in the movement of some organisms (Bell 1991; Blanche *et al*. 1996; Anderson *et al*. 1997; Conradt *et al*. 2000; Bengtsson *et al*. 2004), but sometimes this can lead to an oversimplification, since the resulting loops are orientated in one direction only. A solution to overcome this problem is to consider more memory to the CRW, by introducing a long-range correlation between the turning angles (Blanche *et al*. 1996; Wu *et al*. 2000; Preisler *et al*. 2004) or even by adding correlation between turning angles and the variable speed (Wiktorsson *et al*. 2004). These diffusion models can produce leptokurtic dispersal at small temporal scale; however, if there is only a short-range dependence in the CRW then at large times, the dispersal converges to normal diffusion (Tchen 1952; Bovet & Benhamou 1988; Codling & Hill 2005). A different approach, inspired again by observed data (Lewis *et al*. 1992, 2006; Ramos-Fernandez *et al*. 2004), is to divide individual behaviour into two events: periods of activity (jumps) and periods of inactivity in which the organisms are stationary (waiting times; Metzler & Klafter 2000; Wiktorsson *et al*. 2004). Derived originally from physics studies, when simulating the transport of particles, this description is the basis of the continuous-time random-walk (CTRW) model, in which the length of a given jump, as well as the waiting time between two successive jumps, follows a certain distribution. Depending on the form of this distribution, the CTRW model can generate different types of diffusion such as sub-diffusion, Lévy flight or Lévy walk (Metzler & Klafter 2000, 2004). These models have found application in many biological systems (Levandowsky *et al*. 1997; Atkinson *et al*. 2002; Ramos-Fernandez *et al*. 2004; Zhang *et al*. 2007), although more recent studies have suggested that they may not be so reliable due to difficulties in estimating model parameters from observed movement data (Benichou *et al*. 2006; Benhamou 2007; Edwards *et al*. 2007; Plank & James 2008).

In another explanation of the origin of the anomalous diffusion described by the leptokurtic distribution, the foraging behaviour of individuals is assumed to be Brownian and the leptokurtic movement is caused by heterogeneities. Many studies consider spatio-temporal heterogeneities in which the organism moves fast in some regions and slow in others; such movement being modelled by diffusion equations with a diffusion coefficient that depends on space (Johnson *et al*. 1992*a*,*b*; Turchin & Thoeny 1993; Feltham *et al*. 2002) or by considering temporal changes in individual movement (Skalski & Gilliam 2003; Morales *et al*. 2004; Yamamura 2004; Yamamura *et al*. 2007). On the other hand, heterogeneity among individuals in a population has proven to be another important factor that can influence population dispersal. To incorporate population heterogeneity into dispersal rate, some studies assume that the movement of individuals is Brownian, but the travelling durations of organisms vary following an inverse-gamma distribution (Clark *et al*. 1999) or a gamma distribution (Yamamura 2002) leading to fat-tailed dispersal kernels. Another way to integrate population heterogeneity into dispersal models is to assume that populations consist of several subgroups with different diffusion coefficients (Skalski & Gilliam 2000, 2003; Okubo & Levin 2002). Inside each subgroup, the dispersal is described by a normal diffusion process with a constant diffusion coefficient that varies from one group to another. Owing to this variability, the overall dispersal cannot be described by normal diffusion and assumes a leptokurtic form. The model presented by Skalski & Gilliam (2000), based on two or more subgroups of stream fishes, showed that an increase in the number of subgroups not only significantly improves the model fit but also increases the number of parameters involved in the model. It is therefore reasonable to generalize this idea by considering dispersing populations as having diffusion coefficients distributed according to a continuous probability distribution function.

The present work uses digital recordings of the movement of the slug parasitic nematode *Phasmarhabditis hermaphrodita* as the basis of a theoretical model that quantifies the organism's movement over time and space. *Phasmarhabditis hermaphrodita* is a nematode capable of controlling slug damage in a wide variety of agricultural crops, and therefore a proper understanding of their dispersal abilities is very important for the optimization of this biocontrol system (Wilson *et al*. 1993; Tan & Grewal 2001; Hapca *et al*. 2007*a*). There are very few models of nematode movement in the literature. Earlier work assumes a simple random walk model (Croll & Sukhdeo 1981; Hunt *et al*. 2001) with subsequent effort incorporating a biased random walk (Anderson *et al*. 1997), where an advection term was added to the classical diffusion equation, to simulate the loops in the movement patterns of *Caenorhabditis elegans*. Recently, a two-dimensional correlated random model has been considered to simulate the movement of *P. hermaphrodita* in homogeneous environments (Hapca *et al*. 2007*b*). It is an individual-based model that incorporates a long-range correlation in the step length recorded per unit time and the direction, so that each new step length and each new direction depends on previous directions and step lengths.

Most models have been derived from movement data on plates of agar because, due to their microscopic size, observations on movement of nematodes in soil are difficult to obtain. Considering their role in regulating soil biophysical properties, it is of practical importance to understand how these organisms behave in a structurally heterogeneous environment. In this study, we attempt to understand the behaviour of *P. hermaphrodita* in soil, by undertaking a comparative study between nematode movement in a homogeneous environment (agar plates) and a simulated heterogeneous structure obtained by adding sand particles to the surface of the agar.

## 2. Materials and experimental methods

### 2.1 Experimental treatments, image capture and analysis

Dauer larvae of *P. hermaphrodita* were obtained from Becker Underwood (Littlehampton, UK) in a formulated product (Nemaslug) consisting of partially dehydrated dauer larvae bounded 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 were 1 mm in length and 0.03 mm in width.

In all treatments, the experimental arena comprised a 9 cm diameter Petri dish containing 1.2 per cent technical agar (Oxoid Ltd, Hampshire, UK). To analyse the effects of structural heterogeneity on nematode movement, two experimental treatments were selected: a nematode on (i) homogeneous agar (control, homogeneous) and (ii) agar with a monolayer of sand grains (test plate, heterogeneous). One active nematode was picked from the water where it was rehydrated, and transferred, by means of an eyelash mounted on a syringe pick, to the centre of the agar plate (Rodger *et al*. 2003). The same procedure was used for the second treatment where the heterogeneous structure was introduced by adding a monolayer of acid washed sand grains (0.1–0.3 mm particle size) uniformly dispersed on the plate by means of a very fine sieve until approximately 30 per cent of the surface area was covered.

An Axio MRc Zeiss camera was attached to a Leica microscope MZ16 and connected to the computer monitored movement of nematodes with an image field size of 3×4 cm. It was not possible to consider a larger image field and still be able to resolve the nematode. The position of each nematode was recorded at 8 s intervals over a period of 15 min, which was equivalent to 100 frames. The experiments were done under a constant temperature of 19°C and 25 replicates were prepared for each treatment.

Images were recorded with Axio Vision v. 3.1 (Carl Zeiss Inc., USA), then analysed using an image analysis software incorporating a tracking algorithm (Image-Pro Plus v. 5, Media Cybernetics Inc., MD, USA), in order to obtain the *x* and *y* coordinates of the nematode at each time point. Matlab v. 7.3 (The Mathworks Inc., USA) was used to plot the final digitized nematode trail, to compute the distance travelled between frames (step length) and the value of the turning angle between successive steps of the movement.

### 2.2 Data analysis

Direct observations of nematode movement on agar plates suggested a large variation in individual behaviour among the population. To measure this variation, we computed the individual's mean speed (mm s^{−1}) by calculating the average distance between successive positions recorded in an 8 s time interval, chosen on the basis that in 8 s the nematodes travel an average distance that is approximately equal to their body length (1 mm). Shorter time scales would have introduced artefacts due to changes in body shape during movement, and longer time scales would have reduced the amount of data available for analysis.

The values of the step length were then used to check if nematode activity on the agar plates was constant during the recording time. EViews v. 4 software was used to perform the Dickey–Fuller unit root test (Gujarati 1995) for each nematode trail in order to decide if the time series formed by the nematode step length was stationary over time. When performing the test for stationarity, a minimum of 50 points in the time series are required. In our experiments, due to the small image field size, some nematodes had moved out of the image view before the 15 min recording time. Therefore only 18 and 17 trails from the first (homogeneous) and second (heterogeneous) treatments, respectively, provided us with the 50 time points required and only those ones were checked for stationarity.

To determine the degree of correlation between steps in the individuals' movement we first used MINITAB v. 15 to compute the step length autocorrelation *ρ*_{ΔR}(*k*) as a function of lag *k* between steps. If there is no seasonality trend in the time series, in general, the autocorrelation function tends to decrease with the increase of the lags, becoming close to zero after a certain number of lags (Hanke & Wichern 2007). For each individual trail, we determined the order of autocorrelation in the step length time series as given by the largest lag *k*≥0 for which the corresponding coefficient of correlation was significantly different from zero. Then based on this autocorrelation order, the trails were grouped into two categories: trails displaying at most a first-order autocorrelation and trails with at least the second-order autocorrelation coefficient significantly different from zero. A 2×2 contingency table was produced by recording the number of trails in each category and for each treatment, and Fisher's exact probability test (Zar 1999) was used to determine significant differences in the autocorrelation function between treatments. Second, we generated nematode turning angle distributions by computing the angle between successive steps of the trails. The angular range chosen was from −*π* radians to *π* radians, and all angles were computed with respect to the previous direction. We adopted the convention that negative values correspond to the right turns and positive values to left turns. Circular statistics (Batschelet 1981; Mardia & Jupp 1999; Zar 1999) were applied to the values of the turning angle in order to compute the mean turning angle *a* and the mean vector length *r* for each nematode trail. The vector length *r* takes values between 0 and 1, and a value of *r* close to 0 would indicate that there is no preferential turning direction in the movement; on the contrary, if *r* becomes close to 1 then the turning angles are concentrated around a mean value and we deduce that the individual has a preferential turning direction given by the mean turning angle *a*. Next, the individuals' mean angles *a* and mean vectors *r* were used to perform a second-order analysis of angles in order to obtain the population mean angle and vector length for each treatment. Then Hotelling's procedure for parametric one-sample second-order analysis of angles (Zar 1999) was used for each treatment separately to test whether there was a mean population direction in the movement. As previously, in order to obtain reliable results for the step length autocorrelation function and the turning angle distribution, we considered those trails containing at least 50 time steps for the analysis, i.e. 18 and 17 trails for the homogeneous and heterogeneous treatments, respectively.

Finally, nematode step lengths and features of the turning angle distribution such as mean cosine *c*, mean sine *s*, mean vector length *r* and mean angle *a* were used for each trail separately to determine the individual's diffusion coefficient. The formula (A 5) presented in appendix A, to calculate the coefficient of diffusion, applies for CRW in which the turning angles are statistically independent; it was first derived by Kareiva & Shigesada (1983), then developed by Marsh & Jones (1988) and Benhamou (2004, 2006). From the previous analysis, the experimental results suggested that the movement in heterogeneous structure was closer to the CRW, unlike the homogeneous treatment where longer range dependence was found in the steps of the movement. Therefore, in the first instance, the diffusion coefficient was calculated only for those trails and then used to validate the diffusion model developed in the study. A gamma distribution was fitted to the distribution of the individuals' diffusion coefficients using maximum-likelihood methods for parameter estimates (Cramer 1986) and the Kolmogorov–Smirnov test was used for goodness-of-fit test (Statistics Toolbox, Matlab v. 7.3).

## 3. Experimental results

Examples of nematode trails produced in each of the two experimental treatments are presented in figure 1. A large variation in individual behaviour is reflected in their mean speed (figure 2), ranging from nematodes moving at the slowest speeds (figure 1*a*) to those moving fastest (figure 1(iv)).

Next, the Dickey–Fuller test for stationarity showed that, in both environments, the statistical properties of the movement of the majority of nematodes remained constant during the recording time, with only three trails in the homogeneous treatments and one trail in the heterogeneous treatment not satisfying the stationarity condition.

At the individual level, the step length autocorrelation function in the heterogeneous structure was significantly different from the control. Table 1 shows that, in the presence of sand particles, the autocorrelation in step lengths for the majority of nematodes (11 out of 17) is not significantly different from zero for lags greater than one, while in the homogeneous environment 13 out of 18 nematodes were displaying at least a second-order coefficient of correlation significantly different from zero. The results of Fisher's exact probability test showed that the difference between the two treatments was significant (*p*=0.0313).

The distribution of the turning angle corresponding to the two treatments is presented in figure 3. In the homogeneous environment, the distribution is peaked around the value 0 corresponding to a forward persistent movement. By contrast, in the heterogeneous environment it becomes more flat around this value, with an increase in frequency around −*π* and *π* corresponding to backward movement. The second-order analysis of the turning angles showed that in the homogeneous environment the forward persistence in the movement was highly significant (, *F*=6.60, d.f.=2,16, *p*=0.002), while in the heterogeneous environment this tendency was less significant (, *F*=4.13, d.f.=2,15, *p*=0.037).

The previous statistical analysis shows that the individuals' steps are correlated with both environments; however, the movement becomes less correlated with the presence of heterogeneous structure. Therefore, we assume that in the heterogeneous environment the movement can be described by a first-order CRW. It has been proved that after a sufficiently large number of steps, the CRW converges to normal diffusion (Tchen 1952; Bovet & Benhamou 1988; Codling & Hill 2005), with a diffusion coefficient derived in appendix A. The resulting distribution of the diffusion coefficients in the heterogeneous structure is presented in figure 4*a*. A Kolmogorov–Smirnov test showed a good fit (figure 4*b*, *p*=0.89) to this distribution by a gamma distribution with maximum-likelihood estimates *ν*=0.57±0.13, *λ*=0.21±0.07.

## 4. Model description

### 4.1 A stochastic process in two dimensions associated with nematode movement

Nematode movement is described in terms of a two-dimensional stochastic process (*X*_{t}, *t*≥0) defined on the space *Ω*, representing the population of nematodes. Thus, for each nematode *ω*∈*Ω*, corresponds to the position of the nematode at a certain time *t*≥0. We aim to determine the probability density function (PDF), , associated with this stochastic process that will be used to quantify nematode dispersal.

Based on the experimental results, we assume that the process (*X*_{t}, *t*∈(0, *T*)) has stationary increments up to a certain time, *T*, and that the individuals' movement is governed by normal diffusion with a diffusion coefficient, *D*, that varies among the individuals. If we denote by the Brownian motion associated with normal diffusion, then the corresponding PDF is described by the Gaussian kernel(4.1)

We develop a heterogeneous population model that considers the contribution of each individual to the movement of the population as a whole. The case where *D* follows a discrete distribution has already been investigated by Skalski & Gilliam (2000). When *D* follows a continuous distribution with a PDF, *f*_{D} on (0,+∞) the stochastic process (*X*_{t}, *t*∈(0, *T*)) can be defined by its finite-dimensional distributions as follows:(4.2)for any *n*∈*N*^{*} and . The finite-dimensional distributions above satisfy the Kolmogorov's consistency criterion and therefore the stochastic process (*X*_{t}, *t*∈(0, *T*)) is well defined (Applebaum 2004).

Prompted by the experimental results, we present the case when the diffusion coefficient *D* follows a gamma distribution, *γ*(*λ*, *ν*). Replacing *f*_{D} in (4.2) by the gamma distribution we obtain(4.3)where *K*_{ν} is the modified Bessel function of the second kind (Watson 1966).

Some important properties of this stochastic process are detailed in appendix B. It can be easily deduced that, in common with Brownian motion, the process (*X*_{t}, *t*∈(0, *T*)) preserves the property of self-similarity with a Hurst parameter *H*=1/2, and the mean square displacement is linear in time. Next, the increments are stationary in time, but unlike the Brownian motion, they are not independent (B 2). It can be shown that the stochastic process described by (4.3) can be characterized by a leptokurtic distribution. This can be deduced from the asymptotical behaviour and the behaviour at origin of the PDF (4.3). Using the asymptotic properties of the Bessel functions (Watson 1966) one can obtain that for large we havewhich shows a longer tail compared with the normal distribution, while close to the origin, i.e. ,for 0<*ν*<1 and for *ν*≥1. In the particular case, *ν*=3/2, we rediscover the exponential kernel . Alternatively, by computing the kurtosis of the marginal distributions of the stochastic process we obtainwhich again is a characteristic of the leptokurtic distribution.

Finally, for a time step Δ*t*>0, the step lengths form a long memory process in the sense defined by Beran (1994) (see appendix B for further details).

## 5. Model application and results

The theoretical model presented above was first validated against the data obtained from nematode movement in the heterogeneous structure. The PDF of nematode dispersal was derived from (4.3) with gamma distribution parameters *ν*=0.57 and *λ*=0.21, obtained previously. The resulting profile is plotted in figure 5 for three different times, 5, 10 and 15 min, and then compared to the normal diffusion with diffusion coefficient *D*=0.118, computed by averaging over all the individual diffusion coefficients. Figure 6 shows clearly the difference between the two models, and how the model developed in this paper is characterized by a leptokurtic distribution.

In order to validate the model, for each trail, we computed the net displacement at 5, 10 and 15 min, where 15 min represented the total recording time. As some of the nematodes were out of the image field before the total recording time, we could not include measures of the net displacement for them. As an alternative, we divided the image area into four sections corresponding to less than 5 mm, 5–10 mm, 10–15 mm and greater than 15 mm distance from the nematode starting point and counted how many nematodes were in each of the four sections after 5, 10 and 15 min from the time they were released. Nematodes that were out of the image field by the end of the recording time were counted in the last section.

The resulting counts are presented in table 2 and compared with the numbers predicted by the new model developed in this paper. The results of the chi-squared goodness-of-fit test showed that for each of the three time intervals, the observed counts were not significantly different from the predicted ones (*p*>0.425). Also table 2 shows that this new diffusion model gives a better fit than a model assuming normal diffusion with diffusion coefficient *D*=0.181 derived from the data (*p*<0.005).

However, these results do not provide any evidence that the new model produces a better prediction than the one developed by Skalski & Gilliam (2000), where the diffusion coefficient was characterized by a discrete distribution. In order to address this, we arranged the values of the individuals' diffusion coefficient in ascending order and divided the population of nematodes into four subgroups, the first subgroup containing four individuals corresponding to the smallest four values in the dataset. The other three subgroups were constructed in a similar way, with each containing six individuals. A diffusion coefficient was calculated for each subgroup, by averaging over the values of the diffusion coefficient of the individuals within that subgroup. This resulted in the following four values: *D*_{1}=0.026, *D*_{2}=0.147, *D*_{3}=0.433, *D*_{4}=1.647. These values were then used to compute the predicted numbers of nematodes in each of the four sections and for the three different instant times, 5, 10 and 15 min. Table 2 shows that, although we have divided the population into four discrete subgroups, this model is inferior to the new model.

Finally, in order to complete the analysis we attempted to validate the model with the data obtained from nematode movement in homogeneous environment. Despite the fact that the statistical analysis presented previously showed that movement on homogeneous plates is characterized by a longer range correlation compared with heterogeneous structure, we employed again the CRW techniques described in appendix A, to derive the values of the individuals' diffusion coefficients. Kolmogorov–Smirnov test showed a good fit (*p*=0.92) of the distribution of individuals' diffusion coefficients to a gamma distribution with maximum-likelihood estimates *ν*=0.84±0.21 and *λ*=0.19±0.061, values that were next implemented into the model. The model was validated as previously by counting the number of nematodes in each of the four sections and for the three different times and then compared with the ones predicted by the model. The results of the chi-squared goodness-of-fit test (*p*=0.522, 0.402 and 0.902 at *t*=5, 10 and 15 min, respectively) showed a good fit of the model to the data, proving that the first-order CRW techniques (appendix A) were reliable when applied to random walks characterized by a longer range correlation.

## 6. Summary and conclusions

In this paper, we develop a diffusion model for the dispersal of a heterogeneous population that can be used to predict nematode movement in soil. The findings obtained from the experiments on agar plates can be extrapolated to the soil environment, where irrespective of how complex the structure is, the population will preserve the variation in individual behaviour observed in the experiments. We do not discount the possibility that the nematode will exhibit a correlated movement in soil. However, we have seen that in the presence of heterogeneous structure, this is a short-range dependence and therefore at large times it converges to normal diffusion. Therefore, the model assumes normal diffusion at the individual level but shows that, due to the heterogeneity observed in the experiments, the diffusion of the population becomes leptokurtic.

While the general finding in relation to variation between individuals is not in itself novel (Clark *et al*. 1999; Skalski & Gilliam 2000, 2003; Yamamura 2002), the diffusion model presented here is new and of significant relevance for the ecological sciences. It is an obvious extension of the diffusion model that can be used to predict spatio-temporal distribution of moving organisms, where classical models fail as they do not consider that individuals within a population have different movement tendencies. Therefore, the techniques are relevant for such heterogeneous populations where individuals disperse with different diffusion coefficients that follow a continuous distribution. When describing heterogeneity effects on dispersal by using continuous distributions, most authors use the inverse gamma distributions or gamma distribution, for describing the fluctuation of dispersal parameter ‘*α*=*Dt*’ (Skellam 1951*a*,*b*; Clark *et al*. 1999), but without specifying whether fluctuation in the diffusion coefficient *D* or the travelling duration *t* is more important. Yamamura (2002) and later Tufto *et al*. (2005) considered the fluctuation of the travelling time to be more important, especially for predicting long-term dispersal. In both studies the travelling time (time to settling) follows a gamma distribution, the resulting dispersal model being cited in the literature as the gamma model. The PDF describing the gamma model is similar to equation (4.3) developed in this paper; however, the essential difference is that in Yamamura (2002), the model predicts only the spatial distribution due to the fact that an integration over time is made. Yamamura (2002) mentioned in the discussions that a diffusion coefficient variable within the population can be assumed as well, but might not be very realistic when modelling long-term dispersal, especially when the individual diffusion coefficient depends on the environmental conditions that are likely to vary at large temporal scales. Therefore, in a later study Yamamura (2004) and Yamamura *et al*. (2007) considered a fluctuating diffusion coefficient over time, and therefore again, the time does not appear explicitly, but is hidden in the parameter *ν* that consequently increases with the increase of time, hence a PDF different from (4.3). Supported by experiments, here we adopted the assumption that the individual diffusion coefficient is constant over time but varies among the population and so we developed a spatio-temporal model that predicts how the distribution evolves over space and time simultaneously, and this arises from the fact that in (4.3) we integrate over the diffusion coefficient and not over time. We think that it is reasonable to consider both temporal changes in individual behaviour as well as differences among individuals as realistic assumptions for a long-term dispersal model, and this will be carried out in a future study.

The choice of the gamma distribution in equation (4.3) was inspired by the experimental results on nematode movement, where the values of individuals' diffusion coefficients were fitting well this distribution. The advantages of this theoretical distribution lie in the fact that because two parameters are involved, it can be easily fitted to a large range of distributions obtained from experimental data. However, should the diffusion coefficients derived from experiments be consistent with a different distribution, the techniques that we have presented in this paper, in particular formula (4.2), can be easily applied to any continuous distribution. The only problem that can occur is that for a different distribution it may not be possible to deduce the PDF associated with the stochastic process explicitly. If this is the case, numerical approaches can be used to calculate the integral in (4.2). In this study, the gamma distribution leads to a convenient and interesting formula of the PDF described by (4.3) that can be used to further investigate the theoretical properties of this stochastic process. However, due to the fact that it is expressed in terms of Bessel functions, we must resort to numerical approximations to calculate the predicted probabilities, as was the case in table 2.

The model presented here has been successfully validated using the data obtained from nematode movement in heterogeneous structure and it proves that we can only obtain a reliable prediction by choosing a continuous distribution of the diffusion coefficient. The previous models in which the population was divided into two or more subgroups of individuals with different diffusion coefficients among the subgroups can be applied when we deal with different species from the same population and when the differences are very distinguishable and so it is easy to place the individual in one category or another. This was the case in the study done by Skalski & Gilliam (2000), where the model was validated on a population formed by two species of fishes. But when individuals with different behaviour come from the same species, as is the case here, it is much harder to group them into a finite number of subgroups. In this study, we showed that by grouping the individuals into four subgroups we obtain a model with four parameters that is inferior to the new model that involves only two parameters. While it is true that by increasing the number of subgroups we can improve the model fit, this is just a way of approximating a continuous distribution using a discrete one. However, the use of successive subdivisions introduces its own limitations as the number of required parameters becomes prohibitive. Therefore, the only way forward is to assume a continuous distribution of the individuals' diffusion coefficients.

In this study, the distribution of the diffusion coefficient was obtained by recording the individual movement on agar plates at the same time scale for all the nematodes. Related to this procedure, there are two points that need to be discussed here. First, when trying to derive the value of the diffusion coefficient assuming a CRW in the movement patterns, this value can be significantly affected by the scaling rate (Bovet & Benhamou 1988; Codling & Hill 2005). In particular, when individuals travel with different movement rates, choosing a right time scale for each individual based on its ability to move can give better estimates of the diffusion coefficient. Here, we have chosen to work with the same time scale for all the individuals, in order to be consistent with the statistical analysis of the nematodes' trail, including the analysis for stationarity and for the autocorrelation function of the step length. The values of the individual diffusion coefficients obtained by using the same time scale for all individuals were shown to be reliable as they provided good parameter estimators to validate the model. It is not the aim of this study to investigate how the use of different time scales for different individuals would affect the parameter estimates and whether these estimates provide a better fit to the dispersal data when implemented in the model. This analysis can lead to interesting results and it will be considered in a future study. On the other hand, the way the experiments were designed to obtain movement data and implicitly the model that we have produced based on these data does not take into account possible ‘alterations’ in individuals' movement caused by interactions between organisms, or the presence of chemical signals or other biological processes including death rate or reproduction in relation to environmental conditions. In its simple form, the study presented here presents a framework for a new diffusion model that successfully quantifies the dispersal of heterogeneous populations in the absence of any other physical or biological processes. These important aspects need to be addressed and added carefully in future studies if we want to have a detailed picture of population dynamic over time and space.

## Acknowledgments

We thank P. Budha for the image recordings of nematode trails and B. Sleeman for helpful discussions related to the theory of Bessel functions. We are also grateful to the anonymous referees for their constructive comments. This research was funded by the Biotechnology and Biological Sciences Research Council (grant no. BBSB01065).

## Footnotes

- Received June 17, 2008.
- Accepted July 28, 2008.

- © 2008 The Royal Society