A predictive model and a field study on heterogeneous slug distribution in arable fields arising from density dependent movement

Patch formation

The individual-based simulations produce, as an immediate result, the position of all slugs at any given time t; an example is given in Fig. 2a in the case of ‘small’ field (10times 10) m. Analysing information describing the system state when it is presented in this format can be difficult. It is a particular problem when comparing it with field data where the system state is quantified by the slug trap counts (regarded as a proxy for the local slug density) at selected spatial locations, e.g. in a rectangular grid (cf. Fig. 1a). We therefore split the spatial domain into 100 bins by dividing its linear size both in x and y into 10 equal intervals and calculate the number of slugs in each bin, i.e. the population density at the location of each bin. In order to make the results more accessible for the visual perception, we then show the binned numbers as a continuous plot of the population density. As an example, Fig. 2b shows the distribution of the population density corresponding to the simulation data shown in Fig. 2a.

We readily observe that the simulated spatial distribution of slugs is apparently heterogeneous. Two questions immediately arise here as to (i) whether this heterogeneity is different from the heterogeneity of a purely statistical origin and (ii) how the density distribution evolves in time. To answer these questions, Fig. 3b,c show the simulation results after a series of increasing time periods using the same initial condition (Fig. 3a) used for Fig. 2. It is readily seen that the distribution evolves with time and the degree of heterogeneity (e.g. as described by the difference between the smallest and the largest values of the population density, inferred from the numbers on the colour bar) tends to increase over the course of time resulting in the formation of high density patches (shown by the yellow colour). Also, the size of individual patches changes, with a tendency to increase until it reaches a certain value (see the next section for details). For comparison, Fig. 3e,f show the spatial distribution obtained at the same moments of time as in Fig. 3b,c respectively, but in case of purely random density-independent movement; no high density patches emerge in that case.

Simulations show that the emerging high density patches are dynamic rather than stationary (even in the large-time limit, for more details, see Appendix A.2 in online Supplement Information). No stationary distribution emerges in the course of time. However, inspection of the results shown in Figs. 2b and 3b,c (as well as results obtained in other simulation runs, not shown here for the sake of brevity) reveals that the patch dynamics is rather slow, so that some of the patches roughly preserve their size and location on the timescale of (t=10^4), i.e. about 3 weeks in dimensional units, which is in a good agreement with the field data, see Section 2.

Note that, since our model is inherently stochastic, the emerging spatial distribution will differ in the precise shape and position of the patches between model runs. However, the formation of a distinct patchy structure is a generic property of the system. In this sense, the patterns shown in Fig. 3b,c are typical for the system’s dynamics. Moreover, the formation of the patchy structure appears to be robust and does not depend on the initial conditions. For example, in the case where the initial condition is chosen as a dense release (i.e., all animals are initially inside a single patch), over the course of time the initial patch eventually splits into a number of smaller patches resulting, for the same parameter values as in Fig. 3, in a spatial distribution qualitatively similar to those shown in Fig. 3; see³⁴ for all simulation details.

We now recall that, while the movement parameters in Eqs. (1–5) are determined from the field data with a sufficient accuracy, the value of threshold density d where the movement type switches is only roughly estimated. Therefore, the next step is to investigate whether the formation of a patchy spatial distribution is sensitive to the threshold density. For this purpose, simulations were run with a different value of d. The results are shown in Fig. 4. We observe that the variation in d will not eliminate the heterogeneity, a distinctly patchy spatial distribution develops for all values of d used in Fig. 4. The shape and size of the patches (as is readily seen from the visual inspection of the spatial distributions) as well as the difference between the maximum and minimum values of the population density varies slightly for different d without showing a clear tendency. Patchiness appears to be robust to the value of density threshold in a broad range of d (see also Fig. 9 below), unless the average slug density is much smaller than the density threshold; in this case, slug movement is always density independent and distinct patches never form (apart from purely stochastic fluctuations of a small magnitude, cf. Fig. 3e,f).

A similar question arises about the effect of the perception radius, which value is only roughly estimated. Figure 5 shows the results obtained for different R. We observe that a distinct patchy structure emerges for values of R over a broad range, which includes the range estimated from the field data. However, contrary to the density threshold, the degree of spatial heterogeneity clearly depends on R. The typical size of the patches tends to increase with R while the number of the patches decreases accordingly. For a sufficiently large R, a single high density patch is formed, cf. Fig. 5c. This effect of the increase in R can be explained as follows. The perception radius is, by its definition, the distance within which slugs react to each other by slowing down their movement. It is a characteristic length of the population distribution, with the meaning similar to the correlation length. Slowing down of slug movement eventually leads to their numbers building up at the scale consistent with that characteristic length. Note that the average radius of the single patch shown in Fig. 5c is about 2–3 m, and this is consistent with the used value (R=3).

Patchiness quantification

We now complement the visual inspection of the patchy pattern with a more quantitative assessment. There are several measures or indices that are used in statistical ecology for this purpose³⁵. In particular, the Morisita index³⁶ (I_M) provides a measure of how likely two individuals randomly selected from a given spatial domain are found within the same bin (e.g., quadrat) compared to that of a random distribution. It can be shown³⁷ that (I_M=1) if the individuals are distributed randomly (with a constant probability density) and (I_M>1) if the individuals are aggregated for reasons other than purely statistical ones. The Morisita index has been widely used to quantify the heterogeneity of the spatial distribution^38,39,40. It can be calculated as follows:

where (n_k) is the number of individuals in the kth bin, Q is the total number of bins (quadrats) and N is the total number of individuals.

Figure 6 shows the Morisita index calculated at each time step for a few cases with a different total number of slugs. Note that, since the movement of any individual slug is a stochastic process, (I_M) is a stochastic quantity. In order to make sure that any tendency in (I_M) to change is not obscured by random fluctuations (which can be of considerable amplitude), Fig. 6 shows (I_M) averaged over ten simulation runs.

It is readily observed that, in each case shown in Fig. 6, starting from (I_M=1) at (t=0) (which corresponds to our choice of a uniform random initial distribution), the Morisita index then shows a clear tendency to increase on average (apart from the random fluctuations) until approximately (t=8000) when it stabilizes at a certain value (I_M^*>1). We therefore conclude that (i) the spatial patterns obtained in our model are self-organized, i.e. caused by interactions between the individual slugs and not by purely random, statistical effects, and (ii) in the course of time, the system reaches a dynamical equilibrium so that the Morisita index stops growing. For comparison, the red curves in Fig. 6 show the Morisita index obtained in the corresponding cases of purely random density independent movement when no high density patches are formed (cf. Fig. 3e,f).

Note that the Morisita index tends to increase slightly with an increase in the average slug density (i.e., for a fixed size of the spatial domain, with an increase in the total number of slugs N). Indeed, the value (I_M^*) at which the patchiness stabilizes after a long time period rises somewhat with an increase in N, cf. Fig. 6a–c.

In order to provide an overall description of the emerging heterogeneous distribution, with a focus on the density dependence of the properties of the emerging patchy pattern, we use the Taylor’s Power Law aggregation index²⁵. It is well known⁴¹ that, for populations of many different species, the mean (m) and the variance (v) of population numbers in a sample are not independent but related by a power law:

where (alpha) is a coefficient and exponent (beta) is called the aggregation index. In the case where a species has a uniform spatial distribution, (beta) tends towards zero; for a purely random distribution (e.g., described by Poisson distribution), (beta =1). Values (beta >1) reflect progressively greater aggregation, i.e. formation of patches in the field resulting from self-organized, density dependent dynamics of the system.

By writing Eq. (9) on the logarithmic scale:

values of (alpha) and (beta) can be established by fitting (10) to relevant data; in particular, the aggregation index (beta) is defined as the slope of the regression line.

Figure 7 shows the aggregation index calculated for the patchy spatial patterns obtained in simulations. Recall that, when starting with a random uniform distribution, it takes a certain time for the patchy structure to develop. Correspondingly, in each simulation, the system was allowed to evolve over a certain time (t^*) before the population was binned and the mean and the variance of the distribution were calculated. The (a) and (b) panels in Fig. 7 are obtained for (t^*=10^3) and (t^*=10^4), respectively. We readily see that in both cases (beta >1) ((beta =1.066) and (beta =1.173), respectively) confirming the self-organised, inherent nature of the spatial patterns. Note that (beta) is larger in Fig. 7b, that is for a larger (t^*), which is consistent with an earlier observation that the patchy structure becomes fully developed by (tsim 8000).

Evaluating trap counts

In the above, we have shown that our IBM model parameterized using field data on individual slug movement produces a distinctly heterogeneous, patchy spatial distribution. The degree of aggregation, both in terms of the Morisita index and Taylor’s aggregation index is higher than it would be due to purely stochastic reasons. The emerging patchy structure is self-organized in the sense that it emerges not due to the effect of external factors but due to an inherent property of the system such as the density dependent slug movement.

The simulated spatial patterns exhibit properties similar to the distribution of slugs in the field, in particular showing similar values of the aggregation index, which in the field data was found¹⁸ to be in the range (1.09<beta <1.52). There is, however, one difference: while the aggregation of slugs in the field is considered using trap counts, in our simulations we calculated the total size of the slug population in each bin. Obviously, the latter is considerably larger than the former, as in reality only a small fraction of the total population is caught by traps. This explains the difference (almost by an order of magnitude) between numbers shown on the colour bar in Fig. 1 (a typical example of field data) and the numbers shown on the colour bar in Figs. 2, 3, 4, 5 (simulation results).

In order to resolve this discrepancy, one has to evaluate the trap counts arising from a given population density distribution in the area around the trap. This is a difficult problem, as it depends on the trap design and also requires a good understanding of fine details of animal behaviour. Reliable models to calculate trap counts only exist for the basic case of a completely random Brownian movement (which is not the case that we consider in this study) and for pitfall-type traps^42,43,44, i.e. the kind of trap where, once an animal is caught, it is removed from the field. However, the design of refuge traps used in the related field study^18,24 to which we endeavour to link our results was different. Slugs typically remain in a refuge trap for a variable period of time, but they are not removed from the field and can leave the trap at any point. Detailed modelling of trap counts arising for this trap design (e.g. taking into account a variety of slug behavioural responses) requires a separate study and hence lies beyond the scope of this paper. For the purposes of this paper, we use a different, semi-empirical approach to evaluate trap counts. Namely, we take into account that the majority of animals caught by a trap over a given exposition time come from a certain catchment area^42,43,45, so that the closer is the position of a given animal to the trap, the larger is the probability for this animal to enter the trap. Note that, for a simpler (Brownian) movement pattern, the existence and effect of catchment area was earlier shown using rigorous mathematical arguments⁴³. We assume that, for a sufficiently small catchment area, all animals end up in the trap. In order to simulate the catchment area, each of the bins in (10times 10) grid is split to a number of ‘sub-bins’ or cells of a smaller size (see Fig. 8). We then interpret the total number of animals in the central cell as the trap count.

The spatial distribution of trap counts simulated as described above in the case of 25 cells for each bin in the (10times 10) grid (equivalent to the grid of traps used in the field) is shown in Fig. 9. Trap counts are calculated after the initial distribution was allowed to evolve over (t=10^4) to ensure that a distinct patchy structure emerges. It is readily seen that the simulated trap counts correspond well to those collected in the field; see Fig. 1 and also^18,19 for more examples of field data (in particular, the supplementary material in¹⁸ contains data on 167 field visits).

A question arises here regarding whether the aggregation index is sensitive to the characteristics of the proxy used for local slug abundance, i.e. using the total number of animals in each bin or only the ‘trap count’ in its central cell. The relation between the mean and the variance of trap counts obtained in the case of (5times 5) cells is shown in Fig. 10a. We obtain that, in this case, the aggregation index is approximately 1.09, a slightly lower value than in the case with no splitting (cf. Fig. 7), albeit clearly larger than one.

In order to check whether the results are sensitive to the size of the catchment area (i.e. to the number of cells in each bin), we consider an intermediate case of a larger trap catchment area (i.e. splitting to a smaller number of cells). Figure 10b shows the relation between the mean and the variance of trap counts obtained in the case of (3times 3) cells. It is readily seen that the results are similar between the two cases and similar to the baseline case with no splitting (cf. Fig. 7). The aggregation index is consistently larger than one, being approximately between 1.09 and 1.17 for all cases considered, hence indicating a self-organized patchiness (irreducible to purely statistical reasons) in a good agreement with the field data.

Source: Ecology - nature.com