A predictive model and a field study on heterogeneous slug distribution in arable fields arising from density dependent movement
Patch formationThe 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.Figure 2Example of simulation results and their visualization obtained at (t=100) for (R=1) (meters) and the total number of slugs (N=10^4) in a (10times 10) m field. The chosen threshold density is equal to the average slug density, i.e. (d=100) (slugs/m(^{-2})). Other movement parameters are given in the text, see the lines after Eqs. (1), (3) and (5). (a) Positions of all individual slugs, (b) the corresponding density distribution reconstructed from the bin counts (see details in the text) using linear interpolation.Full size imageWe 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.Figure 3The spatial distribution of (N=10^4) slugs shown at different moments of time: (a,d) (t=0), (b,e) (t=10^3), (c,f) (t=10^4). Distributions in the upper row (a–c) are obtained using the density dependent movement model (as is described in “Methods” section). Parameter values are the same as in Fig. 2. For comparison, the lower row (d–f) shows the distributions obtained in the case of density independent movement. While patches of high density (shown by yellow colour) emerge in the course of time in the case of density dependent movement, they do not emerge in the case of purely random, density independent movement.Full size imageSimulations 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; see34 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).Figure 4The spatial distribution of (N=10^4) slugs at (t=10,000) simulated for different values of the threshold density: (a) (d=80), (b) (d=100) and (c) (d=120) (slugs/m(^{-2})). Other parameters are as in Fig. 2.Full size imageA 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).Figure 5The spatial distribution of (N=10^4) slugs at (t=10^4) simulated for different values of the perception radius: (a) (R=0.5), (b) (R=1) and (c) (R=3). Other parameters are as in Fig. 2.Full size imagePatchiness quantificationWe 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 purpose35. In particular, the Morisita index36 (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 shown37 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 distribution38,39,40. It can be calculated as follows:$$begin{aligned} I_M = frac{Q}{N(N-1)}sum _{k=1}^{Q}n_k(n_k-1), end{aligned}$$
(8)
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.Figure 6The mean Morisita Index from 10 simulation runs for different number of slugs: (a) (N=2.5cdot 10^3), (b) (N=5cdot 10^3) and (c) (N=10^4). Here (R=200) and (d=10), other parameters are as in Fig. 2. The red curves show the Morisita index obtained in the corresponding cases of a purely random individual movement, i.e., without any density dependence.Full size imageIt 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 index25. It is well known41 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:$$begin{aligned} v = alpha m^{beta }, end{aligned}$$
(9)
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:$$begin{aligned} log (v)=alpha + beta log (m), end{aligned}$$
(10)
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).Figure 7The variance of bin populations plotted against the mean bin population on a grid of 100 bins shown on a log-log scale. Each point is calculated from a single simulation and the total population is varied between simulations from (N=300) to (N=9900) in intervals of 300. The density dependent parameters are (d=1) (equal to the average slug density) and (R=2). (a) (t=1000), the regression equation (10) is (log (v)=1.066log (m)-0.1774), (r^2=0.9686), (b) (t=10,000), (log (v)=1.173log (m)-0.4551), (r^2=0.9427).Full size imageEvaluating trap countsIn 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 found18 to be in the range (1.09 More