Beyond resource selection: emergent spatio-temporal distributions from animal movements and stigmergent interactions

A principal concern of ecological research is to unveil the causes behind observed spatio-temporal distributions of species. A key tactic is to correlate observed locations with environmental features, in the form of resource selection functions or other correlative species distribution models. In reality, however, the distribution of any population both affects and is affected by those surrounding it, creating a complex network of feedbacks causing emergent spatio-temporal features that may not correlate with any particular aspect of the underlying environment. Here, we study the way in which the movements of populations in response to one another can affect the spatio-temporal distributions of ecosystems. We construct a stochastic individual-based modelling (IBM) framework, based on stigmergent interactions (i.e. organisms leave marks which cause others to alter their movements) between and within populations. We show how to gain insight into this IBM via mathematical analysis of a partial differential equation (PDE) system given by a continuum limit. We show how the combination of stochastic simulations of the IBM and mathematical analysis of PDEs can be used to categorise emergent patterns into homogeneous vs. heterogeneous, stationary vs. perpetually-fluctuating, and aggregation vs. segregation. In doing so, we develop techniques for understanding spatial bifurcations in stochastic IBMs, grounded in mathematical analysis. Finally, we demonstrate through a simple example how the interplay between environmental features and between-population stigmergent interactions can give rise to predicted spatial distributions that are quite different to those predicted purely by accounting for environmental covariates.


Introduction
The left-hand side shows the three possible pairwise interactions between two populations. On the right is an example network built from these interactions. One might imagine that A and B are competing prey being predated by mutualistic predators C and D.
For this, we assume that, as individuals move, they leave a trace of where they have been on of one population to be in the presence of another population. Since animals look at their 162 surroundings at a distance to make movement decisions, our model allows for individuals to 163 respond to the local average density of nearby marks. 164 Mathematically this situation can be described by writing down the probability f (x, t + 165 τ |x , t) of moving from lattice site x to x in a timestep of length τ . This function f is known 166 as a movement kernel. To construct our movement kernel, we use a generalised linear model 167 to describe the attraction to, or repulsion from, the local average density of nearby marks. A 168 second equation is then required to describe how marks are averaged over space. Finally, the 169 deposition and decay of marks over time is given by a third equation. We now give precise 170 functional forms of these three equations in turn. location x at time t, the movement kernel is given by otherwise. (1) Here K x = x f (x, t + τ |x , t) is a normalising constant ensuring that f (x, t + τ |x , t) is a 176 well-defined probability distribution; if a ij > 0 (resp. a ij < 0) then |a ij | is the strength of 177 population i's attraction to (resp. repulsion from) population j; andm δ j (x, t) represents the 178 average mark density over a radius of δ. Note that Equation (1) fits into the broad category of 179 functions that can be parametrised by integrated step selection analysis (Avgar et al., 2016).

180
The equation for average mark density is where S δ = {z ∈ Λ : |z| < δ} is the set of lattice sites that are within a distance of δ from 184 0 and |S δ | is the number of lattice sites in S δ . Note that Equation (2) 192 193 where N i (x, t) is the number of individuals at location x in population i at time t, µ τ is the 194 amount by which marks decay in a time step of length τ , and ρ τ is the amount of marking 195 made by a single animal in a single time step.

196
Equations (1-3) are not the only available functional forms to describe our stigmergent 197 process. However, the specific form for Equation (1) is advantageous because it arrives in 198 the form of a step selection function (Fortin et al., 2005;Avgar et al., 2016). It thus has the 199 potential to be parametrised by the methods of Schlägel et al. (2019), which deals with step 200 selection for interacting individuals (although here we focus on analysing the emergent features 201 of the model in Equation (3) rather than the question of fitting this model to data.) Equation

202
(2) assumes that marks are averaged over a fixed disc around the individual and was chosen 203 for simplicity, but other options, such as exponentially decaying averaging kernels, are also 204 possible. Equation (3) was, likewise, chosen for simplicity.

205
One drawback is that there is, in theory, no limit on the amount of marks in one location. If

206
it is necessary to account for such a limit, one might exchange the ρ τ N i (x, t) term for something Finally, there is an analogy between marks and resource depletion that enables our mod-214 elling framework to be used in situations where animals both deplete resources and move up 215 resource gradients. The idea is to view the total number of marks in a location, from all the 216 populations, as the extent of depletion of a resource. In this case, each population would avoid 217 'marks' left by either population, as animals will tend to avoid depleted resources. We do 218 not explicitly examine this situation here, but it is a possibility for future investigations and 219 expands the potential applicability of our work.
the maximum local population density across space minus the minimum. We want to find out 234 whether the amplitude ever becomes higher than would be expected from individuals moving 235 as independent random walkers (i.e. when a ij = 0 for all i, j in Equation 1), assuming that the 236 individuals are initially distributed uniformly at random on the lattice. For this, we calculate 237 A i,d (t) in the case a ij = 0 for all i, j (i.e. no mark deposition so no interactions between walkers) 238 and take the average over a sufficiently long time period to calculate the mean to a given degree 239 of accuracy (i.e. so that the standard deviation of the mean is below a pre-defined threshold, 240 determined by the needs of the simulation experiment). We call this mean amplitude A rw (for 241 'random walk'). Then the extent to which the patterns are heterogenous can be determined 242 with reference to this base-line value. In all panels, two populations of 100 individuals each were simulated on a 25×25 lattice, with initial locations distributed uniformly at random on the landscape. Also µ = 0.001 and ρ = 0.01 for all panels (Equation  3). Panel (a) shows a system where two populations form a single, stable aggregation. Here, a 11 = a 22 = 0, a 12 = a 21 = 2, δ = 10l (Equation 1). In panel (b) the populations segregate into distinct parts of space. Here, a 11 = a 22 = 0, a 12 = a 21 = −2, and δ = 5l. In both Panels (a) and (b) the snapshot is taken at time t = 5000τ . Panels (c) and (d) show a situation where one population (blue) chases other (red) around the landscape in perpetuity, with snapshots at two different times. Here, a 11 = a 22 = 1, a 12 = 10, a 21 = −10, and δ = 10l.
sensitive to stochastic fluctuations. We therefore introduce a corrected circular mean which 248 accounts for this, and denote it by c i (t) (notice that this is a location in two dimensions, for we say that a system has become (R,

266
The separation index is a simple metric that is quick to calculate for multiple simulation 267 analysis. However, one could also use more sophisticated measures of range overlap, such as 268 the Bhattacharyya's Affinity (Fieberg & Kochanny, 2005) between kernel density estimators 269 (Worton, 1989;Fleming et al., 2015). Here, though, we will keep things simple, to enable 270 analysis of a broader range of simulation scenarios in a realistic time-frame. simplicity (which is also called Turing pattern analysis, after Turing (1952)).

284
To arrive at a PDE system, we take a continuous limit in both space and time, sending 285 l and τ to 0 such that l 2 /τ tends to a finite constant, D > 0. This is sometimes called the 286 diffusion limit, as D is a diffusion constant, but is also referred to as the parabolic limit (Hillen 287 & Painter, 2013). If we take this limit, and also assume that infinitesimal moments beyond the 288 second are negligible, we arrive at the following system of PDEs (see Supplementary Appendix the average of q(x, t) over a ball of radius δ. Here, we assume that animals move at the same 296 rate, so D is independent of i. It is possible to drop this assumption, and we discuss the effect  It is sometimes helpful to simplify calculations by assuming that q i equilibrates much faster 300 than u i , so that the scent mark is in quasi-equilibrium ( ∂q i ∂t = 0), leading to the following 15 equation for each i = 1, ..., N 302 This assumption says, in effect, that the distribution of marks accurately reflects the space 305 use distribution of the population. When terrain marking is deliberate, its usual purpose is 306 precisely to advertise space use. Therefore this quasi-equilibrium assumption is likely to be 307 biologically reasonable in many realistic situations.

308
The LSA technique enables us to use Equations (4-5) to construct the pattern formation distributed on the landscape), and also whether these patterns begin to oscillate as they emerge.

313
This technique dates back to Turing (1952) and is essentially an extension to PDEs of stability 314 analysis for ordinary differential equations (May, 2019).

315
The emergence of heterogeneous patterns is expected whenever there is an eigenvalue whose has positive real part and a non-zero imaginary part then small perturbations of the homoge-319 neous system will oscillate as they grow, at least at small times. Often (but not always) these 320 oscillations will persist for all times, so give an indication of the likely answer to Question (II).

321
We stress that this is just an indication, though, and that discrepancies may exist between the 322 answer to (II) and whether or not the dominant eigenvalue of M is real. Full analysis of when 323 to expect non-constant stationary patterns in Equation (4-5), or when to expect perpetually 324 changing patterns, requires more sophisticated techniques.

Simulation experiments 326
To give some insight into the sort of patterns that can emerge from our model (Equations 1-3), 327 we perform a variety of simulations in the simple case of two populations (N = 2). Throughout,

331
First, we examine the situation where populations have a symmetric response to one an-332 other, so that a 12 = a 21 = a. For simplicity, we set a 11 = a 22 = 0. In this case the continuum 333 limit PDE system (Equations 4-5) has the following pattern formation matrix (derived in Sup- 336 337 Here, κ is the wavenumber of the patterns that may emerge at small times, if there is an 338 eigenvalue of M with positive real part (i.e. the wavelength of these patterns would be 2π/κ).

339
For our simulation experiments, we fix the scent-marking rate ρ = 0.01 to be a low number and 340 vary the decay rate µ. We consider two different values of a: either a = 2, which corresponds 341 to populations having a mutual attraction, or a = −2, corresponding to mutual avoidance.

342
In either case, the dominant eigenvalue of M is always real (Supplementary Appendix C).

343
Furthermore, it is positive if and only if µ < 0.0064. In other words, this mathematical 344 analysis predicts that the system will bifurcate at µ = 0.0064 from homogeneous patterns 345 (µ > 0.0064) to heterogeneous patterns (µ < 0.0064). This means that if marks remain long 346 enough on the landscape, they will affect movement to such an extent that the overall space 347 use patterns change from being homogeneous (so indistinguishable from independent random 348 walkers) to heterogeneous. This hetergogeneity will be either aggregative, if a = 2, analogous to the example in Figure 2a or segregative, if a = −2, like Figure 2b.

350
To test whether we see a similar change in stability in simulations, we start by simulating  This means that it can be readily used to incorporate the effect on movement of environ-379 mental or landscape features. Suppose that we have n such features, denoted by functions 380 Z 1 (x), . . . , Z n (x). For each k = 1, . . . , n, denote by β k the relative effect of Z k (x) on move-381 ment. Then, to incorporate these into the movement kernel, we modify Equation (1) as follows Red dots denote simulation runs that were not (R, T )-stable (for R/l = 1, T /τ = 7500), whereas those on the purple-to-brown spectrum were (R, T )-stable. This colour spectrum corresponds to the separation index, from aggregative to segregative. Linear pattern formation analysis of the PDEs in Equations (4)-(5) predicts stationary (resp. non-stationary) patterns to emerge in the top-right and bottom-left (resp. top-left and bottom-right) quadrants, which corresponds well with the dot colours. Notice that the top-right (resp. bottom-left) quadrant corresponds to mutual attraction (resp. avoidance) and, likewise, the dot colours indicate aggregation (resp. segregation) patterns. Panel (b) gives a schematic of the between-population movement responses corresponding to the four quadrants in panel (a). An arrow from u i to u j represents attraction of population i towards population j. An arrow pointing out of u i away from u j represents u i avoiding u j .
We use this to investigate the effect on space use of interactions both between populations 385 and with the environment, by considering a simple toy scenario, but one that is based on a 386 particular empirical situation. Specifically, we consider two populations competing for the same 387 heterogeneously-distributed resource, Z 1 (x) (here, n = 1). One population is assumed to be 388 a weaker competitor, so avoids the stronger competitor, whilst the movements of the stronger 389 are not affected by the weaker. Both have a tendency to move towards areas of higher-density 390 resources.

391
In our simulations, each population consists of 100 individuals. We examine three cases.  Table 1. Extent to which analytic predictions agree with our simulation analysis for different choices of R and T . The third column gives the percentage of the simulations from Fig. 4 for which the analytic prediction for stability agrees with that measured from stochastic simulations using our method. The fourth (resp. fifth) gives the percentage for which the stochastic simulations were deemed unstable (resp. stable), for the given values of R and T , but the analytic prediction is stable (resp. unstable), denoted as SU/AS (rep. SS/AU).

R/l T /τ
The first is where the effect of the stronger competitor on the weaker is ignored (so animals 393 are assumed to act independently, which mirrors many basic resource/step selection studies).

394
The second incorporates the effect of the stronger on the weaker's movements, but treats each system. However, the latter occurs at a slightly lower value of µ than for the IBM indicating 418 that a slightly lower decay rate of marks is necessary to overcome the stochastic effects and 419 allow patterns to form. In other words, the stochasticity has a mild homogenising effect. (rather than decreasing as before). The red dots in Fig. 3b show that there is indeed hysteresis 429 in the IBM system, whereby the system is bistable for 0.006 µ 0.0065, a phenomenon 430 that has also been observed in single population aggregation models with a slightly different 431 class of differential equation models (Potts & Painter, 2021). This means that if a population This resource is shown in shades of yellow-green, with darker (resp. lighter) green denoting higher (resp. lower) density of resources. Magenta (resp. blue) dots denote the stronger (resp. weaker) competitor. In Panel (a) the individuals do not alter their movement in response to the presence of others, and we simply see a preference for higher quality resources. In Panel (b), as well as attraction to better resources, the weaker (blue) population has a tendency to move away from the stronger (magenta) population. In addition to this avoidance mechanism, in Panel (c) the magenta population is strongly territorial. This leads to the emergence of interstitial regions where the blue population can access resources that may be quite high quality. we found from the ones tested, inasmuch as the results corresponded to the pattern formation 438 analysis in the highest proportion of cases, N% (Table 1). Notice too that the mutually-439 avoiding populations (with a 12 , a 21 < 0) tend to have much higher separation indices, s * 12 , than 440 the mutually attracting populations, as one would expect.  (Fig. 5a). When we account for the avoidance of the weaker population by the 445 stronger then the model predicts that the stronger population will live where the resources are 446 better, driving the weaker to resource-poor areas (Fig. 5b). This, of course, may ultimately A principal technical obstacle was to seperate-out random noise from actual spatial patterns, 500 be they stationary or fluctuating. The fact that our techniques agreed well with the analogous 501 PDE analysis provides a validation and ground-truthing of the methods, suggesting they are 502 capturing the key features of patterning with good accuracy.

503
Furthermore, even in the relatively simple example situations studied here, our IBM anal-504 ysis revealed some interesting theoretical insights. It appears that segregation patterns emerge 505 in a continuous fashion as a parameter value moves past the bifurcation point (Figure 3a).

506
However, when aggregations emerge, they appear suddenly (Figure 3b), with a small change in 507 parameter value causing a sudden jump from homogeneous patterns to clearly-defined aggrega-508 tions. Moreover, this is accompanied by a hysteresis effect, meaning that identical underlying 509 processes can give rise to either aggregation or homogeneity, depending on the history of the 510 system.

511
As well as using our methods to analyse IBMs, it is conceivable that the same methods 512 may be valuable for analysing pattern formation in empirical data. One would, admittedly, 513 need some rather high quality data: large quantities of co-tagged animals for sufficiently long 514 time periods to observe changes in space use patterns. However, in the present 'golden age' of 515 animal movement data (Wilmers et al., 2015), with ongoing rapid increases in the magnitude 516 and quality of datasets (Williams et al., 2020), it is good idea to ensure the methodological 517 and theoretical tools exist to deal with such data as it emerges. We have not focused on data 518 analysis here, but we encourage researchers to test this idea in future studies if they have such 519 data.

520
On the more ecological side, we have shown how accounting for feedbacks between the 521 movement mechanisms of constituent populations may help explain the emergence of intersti-522 tial regions in territorial animals that could provide safe-havens for weaker competitors. Such 523 patterns have been observed in coexistent wolf and coyote populations in the Greater Yellow-524 stone Ecosystem. There, these interstitial regions have also been observed as refuges for deer 525 (Lewis & Murray, 1993). Although we did not consider the mobility of prey resources in our 526 simple example, one could add extra complexity by considering the attempts of mobile prey to 527 avoid predators, and observe how this affects the spatial patterns. However, for the purposes 528 of our simple illustration, this level of modelling complexity was not required.

529
An important thing to note is that emergent patterns from interacting populations cannot 530 be revealed by correlative models alone. To take the example from Figure Schlägel (2020). This would lead to precisely the sort of IBM studied here, which enables 541 analysis of predicted space use patterns in novel environments.

542
In general, our approach is valuable for predicting the distribution of populations whenever 543 the locations of two or more populations affect the movements of each other (Schlägel et al.,544 2020). This has been observed in a variety of situations. We have already mentioned com-545 petition between carnivores, and indeed the movements of coexistent carnivore populations in 546 response to the presence of others has been measured in various studies (Vanak et al., 2013;547 Swanson et al., 2016). Also the effect of predator movement on prey locations (sometimes 548 called prey-taxis), and vice versa (the 'landscape of fear'), has been documented in a variety of 549 scenarios (Kareiva & Odell, 1987;Laundré et al., 2010;Latombe et al., 2014;Gallagher et al., 550 2017). Despite this, the predominant species distributions models used in ecology tend to 551 not to account for the underlying between-population movement processes and the emergent 552 features that they engender, even in cases where they model species jointly (Ovaskainen & Abrego, 2020). Explicit modelling of the underlying movement mechanisms, as we have done 554 here, would help plug this gap and lead to more accurate description and forecasting of species 555 distributions.