## Abstract

We investigated the phenomenon of population outbreaks in a spatial predator–prey model, and we found that pattern formation and outbreaks occur if the predators have a limited neighbourhood of interaction with the preys. The outbreaks can display a scale-invariant power-law tail, indicating self-organized criticality. We have also studied the system from an evolutionary point of view, where the predator home range is a hereditary trait subjected to mutations. We found that mutation drives the predator home range area to an optimal value where pattern formation and outbreaks are still present, but the latter are much less frequent. We developed analytical approximations using mean field and pair correlation techniques that indicate that the predation strategy is crucial for existence of this optimal home range area.

## 1. Introduction

Reaction–diffusion processes have attracted increasing attention over the last few decades, owing to their fundamental importance to physical, chemical and biological phenomena (Marro & Dickman 1999). In the last decade, the formalisms developed in the study of two-species reaction–diffusion processes have found application in the quantitative modelling of social and ecological systems (Marro & Dickman 1999; Cantrell & Cosner 2003). A common approach to modelling these systems invokes a Markov process with a suitable set of local transition rules. For specific dynamical rules, long-range order can emerge in reaction–diffusion systems. At high dimensions larger than the upper critical dimension, the behaviour of such systems is well described by mean field approximations (Araujo & de Aguiar 2007). This agreement with the mean field description arises because the number of interaction paths between any pair of particles becomes very large. So most particles interact (directly or indirectly) with most other particles. Since the mean field dominates, fluctuations away from the mean do not matter. Pattern formation, intermittency, avalanches and cascades thus do not play a significant role. However, in low-dimensional systems, the number of interaction paths between any pair of particles is reduced. No longer does the mean field dominate the behaviour. For this reason, low-dimensional systems have been studied using a variety of other approaches, e.g. Monte Carlo methods, field theoretical methods and cellular automata (CA). Here, we study pattern formation and self-organized criticality (SOC) in a low-dimensional predator–prey system using a CA model.

In 1987 Bak and co-workers made a seminal contribution towards understanding how complexity arises in CA models (Bak *et al.* 1987; Bak 1996). The pioneering studies of CA by Ulam and co-workers had already made it clear that complex behaviour could arise in low-dimensional CA with local interaction rules (for historical details, see Wolfram (2002)). Bak and co-workers showed that, for suitably chosen interaction rules, long-range correlations and scale-free power-law behaviour could arise via a mechanism that they named SOC. Standard critical behaviour in second-order phase transitions (e.g. in the two-dimensional Ising model) arises when a tunable parameter (e.g. the temperature in ferromagnets) equals the critical temperature. At this special critical point, two characteristic scales in the system become equal in magnitude. Typically, some important quantity that grows exponentially with distance is approximately counter balanced by another quantity, which decays exponentially and with a decay constant equal to the growth constant of the former quantity. This cancellation of two exponentials leads to power-law behaviour, resulting in a diverging correlation length (Stanley 1971). However, in SOC, the dynamical rules in CA with SOC drive the system towards the critical state. A suitably chosen dynamical rule thus plays the role of the tunable parameter (Bak 1996).

A poorly understood issue in the context of CA and SOC relates to the role of the interaction distance. In ecological systems, predators and prey cannot have an infinite area of interaction. Predators (preys) can only detect a prey (predator) that comes within some well-defined radius of vision or smell (Viswanathan *et al.* 1999). Moreover, predators often have a home location and forage only in a finite neighbourhood of that location, defining the range of interaction for the two reacting species. This finite range is a key ingredient underlying the rich and complex behaviour of low-dimensional systems, for the following reason. If the interaction range were infinite, then all prey and all predators would interact with each other, so we would recover the mean field behaviour. By contrast, the finite range of interaction allows fluctuations away from the mean field to play an important role in pattern formation, intermittency, etc. Here, the finite interaction distance is introduced to the predator–prey model by imposing a home range area where the predators can prey.

We focus on predator–prey systems in which movement is restricted, i.e. confined. Random walkers can be confined via attractive potentials in Fokker–Planck equations. The idea that many animals have their movement restricted to a specific area was first introduced by Darwin (1859). Since then, many studies have been carried out in order to understand the mechanisms and consequences of this limited movement (Charnov 1976; Giuggioli *et al.* 2006; Araujo & de Aguiar 2007, 2008; Mitchell & Powell 2007; Borger *et al.* 2008; Goodnight *et al*. 2008; Zhang *et al.* 2009). *Home range* and *territoriality* are two different concepts related to this behaviour. Burt (1943, p. 351) proposed the following definitions: ‘Home range then is the area, usually around a home site, over which the animal normally travels in search of food. Territory is the protected part of the home range, be it the entire home range or only the nest’. Here we consider predators whose movement is restricted by a natural home range, but which are not territorial, in the sense that home ranges of different predators may overlap.

In what follows we address two questions: (i) How does the interaction range affect the behaviour of the predator–prey systems with home ranges? (ii) What effect do mutations and natural selection (possibly at the group level; Wolfram 2002; Goodnight *et al*. 2008; Wade *et al*. 2010) have on the system? In §2, we describe the CA model. Sections 3 and 4 present results for the system without and with mutations, respectively. For fixed home range sizes, we observed that population outbreaks can occur in the system and that their intensity is higher for smaller home range areas. Moreover, the outbreaks are self-organized and their distribution has a power-law decay. When the size of the home range is treated as an evolutionary characteristic, the predator population evolves to a preferential home range size that is of intermediate value. In §5, we develop analytical approximations that allow us to understand this behaviour in terms of the predator strategy. Finally, we present our concluding remarks in §6.

## 2. Model description

We describe the dynamics of spatially distributed predators and preys by a cellular automaton, where space is modelled by a two-dimensional lattice with *N*×*N* sites labelled by their location (*i*,*j*), with *i*,*j*=1,…,*N*. Each site can be empty or occupied by one or more individuals of the same type. The predators are characterized by the size of their home range, which is a circular area of radius *R* where the predator can look for preys. We call *R* the ‘predation radius’ and *v* the number of sites in the home range, which is approximately *πR*^{2}. In a first scenario, all predators will be assigned the same *R*, which will be held fixed throughout the simulation. Next, we will allow changes in *R*, similar to mutations, to study the evolution of this characteristic and its effect on the population outbreaks. Home range size will be treated as a social, or group, trait. Individuals born within a group adopt the home range of the group. Changes in home range size will be allowed only when a new group is founded, as explained below.

The state of site (*i*,*j*) at time *n* is indicated by , where 0 represents an empty site, *X*_{m} the presence of *m* preys, and the presence of *m* predators with predation radius *R*. Time is discrete and, given an initial configuration, the system is updated by choosing a random site and changing it according to the following rules:

If the site is empty, each prey in the nearest neighbourhood of the site has a probability

*h*_{x}of leaving one offspring in that site. If an attempt is successful the process stops, so that at most one prey occupies the site.If the site is occupied by preys, all predators that have the site in their home range can prey there, with probability

*h*_{y}. The predators that are nearest to the site are chosen to prey first; if the predation does not occur, the second nearest predators have the chance of preying and so on. This hierarchy is based on the fact that a predator spends more time in the centre of the home range, where its nest is located (Samuel*et al.*1985). When predation is successful, the prey population decreases by one individual and the predator population at the corresponding site increases by one individual. If the prey population decreases to a single individual and a further predation success occurs, the predator so generated migrates to the site (*i*,*j*), founding a new predator group. The offspring has a probability*μ*of having its predation radius increased or decreased in relation to its parent. The final state becomes or , where*R*is the parent’s predation radius. If, on the other hand, at the end of the predation process, there are still preys on the focal site, each remaining prey has a chance of*h*_{x}of having an offspring in the same site. In this case, the final state becomes , where 1≤*m*′≤2*m*. Notice that this rule updates not only the randomly selected site, but also the sites in its neighbourhood. Figure 1 illustrates the rule.If the site is occupied by predators, , it can be updated to or , considering that each predator has a probability

*d*of dying.

For a predator whose home range has *v* sites, the total chance of success in preying is roughly proportional to *vh*_{y}. If *h*_{y} is kept constant for all predators, the dynamics will clearly favour larger home ranges. However, there is a cost in keeping a large home range and preying far from home. To take this into account we set the probability of success per site as *h*_{y}=*g*_{y}/*v*, where *g*_{y} is a constant parameter.

One relevant characteristic of this model is the absence of a carrying capacity. If we restrict the number of individuals per site to just one, outbreaks do not occur. We have checked that the results reported here are maintained if the carrying capacity is set to 40 or more individuals per site. Another important feature of the model is the hierarchical predation rule that gives preference to predators that are nearest a given prey. This turns out to be responsible for the existence of an optimal predation radius in the evolutionary scenario.

As a first example, we consider the dynamics without mutation, *μ*=0. In this case, the size of the home range *R* is fixed and equal for all predators. We start the simulations considering a random initial condition but ensuring that 50 per cent of the sites are empty, 25 per cent are occupied by one prey and the remaining 25 per cent are occupied by one predator. The probability values associated with rules 1, 2 and 3 are fixed to *h*_{x}=0.4, *g*_{y}=0.4 and *d*=0.2.

Figure 2 shows the time evolution of the average number of preys and predators per site and the spatial distribution at a fixed time. In these plots, the predation radius is *R*=1, so that the number of sites in the home range is *v*=4 and *h*_{y}=*g*_{y}/*v*=0.1. Time is measured in units of lattice size, *t*=*n*/*N*^{2}, and we used *N*=100. The average populations display a series of peaks, or outbreaks, where the population in a single site can reach up to 10^{5} individuals. The spatial distribution shows the populations at *t*=1.5×10^{4}. Notice that, although the average populations do not show any significant peak at this time, we still find up to 1000 individuals in a single site. We tested different grid sizes, *N*={20,50,100,150}, and found the same behaviour, except for a change in the scale of the fluctuations of the average populations.

## 3. Fixed home ranges

One of our main interests was to quantify the outbreaks and to understand the factors that drive the dynamics to such large fluctuations. Our simulations show that the outbreaks are particularly sensitive to changes in the predation radius. Figure 3 shows the probability distribution *p*(*n*) of the number of individuals in a single site for different values of *R*. The distribution *p*(*n*) was generated from the dynamics in the interval 5×10^{3}<*t*≤5×10^{4}.

When *R*=1, the home range has *v*=4 sites. As the radius increases, the number of sites in the home range increases discontinuously at the critical values where *n* and *m* are integers. At , for instance, the number of sites jumps from four to eight. Table 1 displays the first 13 and the 43rd critical values of *R* and the corresponding number of sites in the home range.

In all simulations, the grid size was kept constant at *N*=100. We tested different grid sizes and we observed that *N*=100 is large enough not to change the conclusions for *R* up to 10. The tail of distributions shown in figure 3 can be fitted by a power-law function, *p*(*n*)∝*n*^{α}, if *R* is sufficiently small. For *R*=1 the slope is maximum, *α*=−2.26. For *R* larger than 2.0, the distribution changes gradually from a power law. First, the modulus of the power-law exponent changes from less than 3 to more than 3. If the exponent is larger than 3, then the distribution has a finite variance. So, upon rescaling, the fluctuations will appear to follow Gaussian statistics since the central limit theorem holds when the variance is finite. Finally, the shape of the distribution undergoes a distortion and no longer appears to be a power law. Instead, the distribution seems to decay more rapidly than as a power law. For *R*=10, for example, a good fit for the prey population is provided by the stretched exponential .

Figure 3 also shows the kurtosis (Mandel & Wolf 1995), defined as the ratio of the fourth to the square of the second momentum, , for the same time interval 5×10^{3}<*t*≤5×10^{4}. We observed that the kurtosis seems to become larger by several orders of magnitude for *R*={1,1.42,2}. A diverging kurtosis is a characteristic of power-law distributions when −5<*α*≤0.

This analysis indicates that population outbreaks become less frequent as the radius increases. So fluctuations become less important as the predation radius increases. This result is consistent with what one would expect theoretically, since a larger radius leads to many more interactions such that the mean field approximation becomes better.

## 4. Variable home ranges

The analysis of the previous section revealed that very large outbreaks can occur if the predation home range is small; as *R* increases, the outbreak intensity decreases and becomes independent of *R*. In this section, we test whether the system can evolve towards a preferential home range size. To do that, we treat *R* as a variable parameter that is transmitted through generations but is also subjected to small random changes. As described in §2, when a predator with home range *R* migrates to a new site by consuming the last prey on that site, its home range can remain equal to *R* with probability 1−*μ* or change to *R*±*δ* with probability *μ*/2. The value of *R* is treated as a social feature, which is kept constant in a group of predators inhabiting a site but may change when a new group is founded in a newly invaded site. In what follows, we have fixed the variation in the radius at *δ*=0.1 and the ‘mutation rate’ at *μ*=0.1.

We start the simulations assigning the same predation radius *R*_{0} to all predators. Figure 4 shows the average number of preys and predators as a function of time for two different initial values, *R*_{0}=1 and *R*_{0}=5. For *R*_{0}=1, the density of individuals increases initially because of the frequent outbreaks that occur for small *R*. However, the density decreases again once the predators have adjusted their home range to an evolutionary stable value. Figure 4 also shows a typical spatial pattern of the population and the distribution of home range sizes over space. The local patterns of home range sizes change in time very slowly when compared with the frequency of outbreaks. While the local outbreaks oscillate with period of about 10 time steps, regions dominated by a given value of *R* take about 100 iterations to evolve to a different value. Comparing with the case without mutation (figure 2) reveals that the outbreaks are greatly suppressed by mutation.

Figure 5 shows the evolution of *R* for the two initial conditions. In both cases, the average of *R* converges to *R*≈2.8 with significant occurrences in the interval 2.24<*R*<3.16. The distribution of the number of individuals in a single site is also plotted in figure 3 and is very similar to those of fixed *R* for 2.24≤*R*≥2.83. In conclusion, if the radius is allowed to vary, it evolves to a preferential value around *R*=2.8 where the outbreaks are reduced but are not completely absent.

## 5. Analytical approximation

In order to better understand the existence of a preferential home range size, we consider here a simplified version of the model that is amenable to analytical treatment. We introduce two simplifications: the first consists in restricting the number of individuals per site to just one. This ensures that the dynamical rules change the state of only one site per time step. The second simplification is the restriction of *R* to two values *R*_{1} and *R*_{2} only, where *R*_{2}=*R*_{1}+*r*. With these modifications outbreaks are no longer possible, since the number of individuals is limited to *N*^{2}. However, the new model still allows the study of preference for a home range size.

The strategy is to construct a master equation for this simplified CA and to derive mean field approximations for the average density of prey and predators. As we shall see, the plain mean field approximation is not accurate enough to capture the behaviour of the automaton and we will need to include pair correlations. We will follow closely the methodology presented in Satulovsky & Tome (1994) and de Aguiar *et al.* (2003, 2004).

The simplified dynamics of the system is as follows: there are four possible states for each site in the lattice, empty (0), occupied by a prey (*X*), occupied by a predator with predation radius *R*_{1} (*Y*_{1}) or occupied by a predator with predation radius *R*_{2}=*R*_{1}+*r* (*Y*_{2}). At each time step, preys can reproduce with probability *h*_{x} in a nearby empty site. Predators *Y*_{1} (*Y*_{2}) can prey with probability *h*_{1} (*h*_{2}) in the sites that are inside their respective home ranges. We set *h*_{1}*v*_{1}=*h*_{2}*v*_{2}, where *v*_{k} is the number of sites in the home range *R*_{k}. When predation occurs in a given site, a predator is born in that site and its predation radius has a probability *μ*/2 of mutating. Each predator has a probability *d* of dying, leaving its site empty. Therefore, the updates follow the cyclic transitions . We set the mutation probability to *μ*/2 because there are only two possible values of *R*, in contrast with the general case of many values where the mutation probability is *μ*/2 for changing to both larger and smaller values of *R*.

In this section, to simplify the notation, we relabel the sites (*i*,*j*) with a single index *i*. We call *σ*=(*σ*_{1},*σ*_{2},…,*σ*_{N2}) the state of the whole lattice and *P*(*σ*,*n*) the probability of finding the system in the state *σ* at iteration *n*. The value of *P*(*σ*,*n*+1) can be calculated as
5.1where and ^{i}*σ*_{w} are states that are equal to *σ* at all sites except at the site *i*, whose state is *w*. The sum over *w* is necessary because two different site states can be generated from the same state: observe that *Y*_{1} and *Y*_{2} can evolve to the same empty state, and they both can be originated from the same *X*. and are the probabilities that the transition indicated occurs in one time step. According to the dynamical rules, they are given by
5.2and
5.3where
5.4

In these equations *n*_{i} is the number of preys neighbouring the site *i*. The probability is the probability of predation times the conditional probability *Ω*_{k} that the new offspring is a specific predator with predation radius *R*_{k},
5.5and
5.6where *m*_{ki} is the number of predators *Y*_{k} that are a distance smaller than or equal to *R*_{k} from the site *i*. The term *p*_{k}=1−(1−*h*_{k})^{mki} is the probability that one of the *m*_{ki} predators will prey and the indices *k*≠*k*′ can be 1 or 2.

The mean field equations for the density of preys and predators can be derived with the help of a suitably chosen auxiliary function *f*(*σ*). Let be the ensemble average of *f*(*σ*) at iteration *n*. Using the master equation for *P*(*σ*,*n*), equation (5.1), and the cyclic property of the transitions we obtain
5.7Choosing *f*(*σ*)=*δ*(*σ*_{i},*α*), where *α*={0,*X*,*Y*_{1},*Y*_{2}}, we see that is the probability that a site is in the state *α* at the time *n*. Explicitly,
5.8

The transition probabilities in the averages above involve the number of preys *n*_{i} or predators *m*_{ki} in a given neighbourhood. Since we may write, for example, for *j* in the neighbourhood of *i*, it becomes evident that these averages will involve pair correlations of the form 〈*δ*(*σ*_{i},*α*)*δ*(*σ*_{j},*β*)〉 and higher order correlations as well. If we approximate all these correlations’ by-products of one-site averages, such as , we obtain the mean field approximation. Renaming the probabilities , and we find (see appendix A)
5.9where
5.10where *v*_{k} is the number of sites in the home range with radius *R*_{k} and *v*_{x}=4.

If we keep the pair correlations as independent variables and only approximate third and higher order correlations in terms of one- and two-site terms we obtain the so-called pair approximation. In this case, we need to write down another set of six equations for the independent probabilities of finding neighbouring sites in states *α* and *β*. We shall not write these equations down here and will only present numerical results. In §6, we describe in more detail the calculation of the mean field approximation.

Figure 6 compares the average number of preys and predators for three cases: simplified automaton version, mean field and pair approximation. In all plots, we consider *R*_{2}=*R*_{1}+*r*, where *r* is the minimum value that increases *R*_{2} to the next critical radius value (table 1). According to the automaton, the density of predators *Y*_{2} is greater than *Y*_{1} if *R*_{1}≤4. For *R*_{1}>4, the density of *Y*_{1} is slightly greater, indicating a preferential radius around *R*=4. The crossing of the curves at *R*=4 is not an effect of fluctuation owing to system size or poor statistics: it is robust against the increase in system size and the time interval used to generate the averages.

The mean field approximation is almost insensitive to variations in *R* but the pair approximation is quite accurate (except for *R*=1), demonstrating the importance of the spatial structure in the dynamics. However, even the pair approximation is not capable of determining the preferential radius, since the curves for *Y*_{1} and *Y*_{2} never really cross. Nevertheless, the two curves come very close to each other for *R*≈3.5 and remain so afterwards, indicating a change of behaviour at this value.

The inability of the approximations to clearly point to the critical radius can be understood in terms of the predator strategy. In the automaton, the action of predators on a site containing preys happens in such a way that nearby predators have priority to prey first. If the first neighbours to the prey are unsuccessful, other predators, in order of distance from the prey, try their luck. Consider a site with preys and all the predators trying to feed on them. For those with small *R* the prey is likely to be outside their home range, whereas predators with large *R* have low probability of success in a single attempt, since *h*≃1/*v*≃1/*R*^{2}. Therefore, predators with intermediate values of *R* will be more effective.

The predators’ strategy of giving priority to the nearest sites cannot be included in the pair approximation and indeed is not reflected in its dynamics, confirming the interpretation above. However, when *R* is very large, the spatial fluctuations in the distribution of preys becomes irrelevant and the average number of preys inside the home range compensates the low probability of success per attempt, resulting in an overall probability of success that is independent of *R*, as shown in the figures.

## 6. Conclusion

In this paper, we have considered the population dynamics of spatially distributed preys and predators whose movements are restricted by a home range. We have shown that self-organized population outbreaks occur if the size *R* of the home range is sufficiently small. The distribution of the number of individuals per site follows a power law for small *R*. As *R* increases, the shape of the distribution changes and seems to decay faster than a power law. There are two important facts to note: (i) scale-free power-law behaviour is seen even though the model itself does not possess scale-free properties and (ii) the power-law behaviour is seen not for a special critical value of the predation radius but rather for any sufficiently small value of this radius. Hence, the behaviour is closer to SOC than to standard criticality. We interpret our results as evidence of SOC, such that, for low radii, the system is driven to scale-free behaviour by the dynamical rules. Empirical evidence of such outbreaks has been observed by Cooke & Lorenzetti (2006) and Cooke *et al.* (2009) and power-law distributions of populations have also been found in similar situations (Scanlon *et al.* 2007).

We have also studied the dynamics of the system, when the home range is allowed to mutate whenever a predator invades a site previously occupied by preys, founding a new group. In this case, the population evolves to a dynamical equilibrium where the predators have a distribution of home ranges. We find that the predators do not evolve to maximize their home range. Instead, the average value of the distribution converges to a finite value.

Our model suggests that the reason for this preferential home range size is the strategy of the predators to forage. If the predators forage randomly throughout their home range, no preferential home range emerges. However, if the predators forage preferentially around their home site, a trade-off between large and small home ranges shows up. Indeed, predators with small home ranges will be very effective in preying at the small number of sites in the range, but the chances of finding preys there are small. On the other hand, if the home range is too large the predator effectiveness per site is small, and the preys might be caught by other more effective predators. The predator strategy adopted here is backed up by previous work on animal foraging strategies (Charnov 1976; Pyke 1984; Samuel *et al.* 1985; Borger *et al.* 2008) and home range characterization (Samuel *et al.* 1985; Giuggioli *et al.* 2006; Mitchell & Powell 2007). These qualitative arguments were corroborated by a mean field analysis of the CA under certain simplifying conditions.

## Acknowledgements

This work was partially supported by CNPq, Fapesp and the Consortium of the Americas for Interdisciplinary Science.

## Footnotes

One contribution of 13 to a Theme Issue ‘Complex dynamics of life at different scales: from genomic to global environmental issues’.

- This journal is © 2010 The Royal Society

## Appendix A. Mean field approximation

In this appendix, we calculate the average terms appearing in equation (5.8). We start with , which can be calculated immediately: since we find . The remaining terms are calculated below.

Expanding , equation (5.2), we obtain
A 1where and *j* runs through the nearest neighbours of *i* only. This is the number of preys close to *i*. The terms of this series can be calculated using the mean field approximation. The first term gives
A 2where *v*_{x}=4 is the number of nearest neighbours.

The second term, equation (A 1), is calculated analogously,
A 3Separating the sums over *l* and *j* into *j*=*l* and *j*≠*l* a partial cancellation occurs and we obtain
A 4The other terms can also be calculated in the same way, resulting in
A 5where *f*_{x}=1−(1−*h*_{x}*P*_{X})^{vx}.

This time the calculation is slightly more involved, but proceeds in the same way. We first expand , equation (5.5), in powers of *h*_{1} and *h*_{2} as
A 6
where
and so on. The quantity is the number of predators *Y*_{k} that can prey in the site *i*; the index *j* runs in the corresponding home range centred on *i*. The contribution of the first term is
A 7where *v*_{k} is the number of sites in the home range of predator *Y*_{k}.

The remaining contributions are calculated analogously. When all the terms are put together, we obtain an expression that can be rearranged as
A 8where and . This series is the Taylor expansion of the function 1−(1−*B*)^{v2−v1}(1−*A*−*B*)^{v1}. Therefore,
A 9where
A 10

We start by rewriting
A 11where
A 12and
A 13Using the same procedure described above we can calculate the average, 〈*δ*(*σ*_{i},*X*)*q*_{k}〉, and find
A 14and
A 15where
A 16Now we use the fact that
A 17Substituting equations (A 9), (A 14) and (A 15) in equation (A 17) we find
A 18One possible solution for *χ*_{k} is
A 19and, therefore,
A 20