Abstract
We propose a minimal model of predator–swarm interactions which captures many of the essential dynamics observed in nature. Different outcomes are observed depending on the predator strength. For a ‘weak’ predator, the swarm is able to escape the predator completely. As the strength is increased, the predator is able to catch up with the swarm as a whole, but the individual prey is able to escape by ‘confusing’ the predator: the prey forms a ring with the predator at the centre. For higher predator strength, complex chasing dynamics are observed which can become chaotic. For even higher strength, the predator is able to successfully capture the prey. Our model is simple enough to be amenable to a full mathematical analysis, which is used to predict the shape of the swarm as well as the resulting predator–prey dynamics as a function of model parameters. We show that, as the predator strength is increased, there is a transition (owing to a Hopf bifurcation) from confusion state to chasing dynamics, and we compute the threshold analytically. Our analysis indicates that the swarming behaviour is not helpful in avoiding the predator, suggesting that there are other reasons why the species may swarm. The complex shape of the swarm in our model during the chasing dynamics is similar to the shape of a flock of sheep avoiding a shepherd.
1. Introduction
Many species in nature form cohesive groups. Some of the more striking examples are schools of fish and flocks of birds, but various forms of collective behaviour occur at all levels of living organisms, from bacterial colonies to human cities. It has been postulated that swarming behaviour is an evolutionary adaptation that confers certain benefits on the individuals or group as a whole [1–5]. These benefits may include more efficient food gathering [6], predator avoidance in fish shoals [7] or zebra [4] and heat preservation in penguin huddles [8]. An example is the defensive tactics used by a zebra herd against hyaenas or lions [4]. These defence mechanisms may include evasive manoeuvres, confusing the predator, safety in numbers and increased vigilance [4,9,10]. On the other hand, a countervailing view is that swarming can also be detrimental to prey, as it makes it easier for the predator to spot and attack the group as a whole [1].
Figure 1 gives some idea of the variety and complexity of predator–swarm interactions that occur in nature. A common characteristic is the formation of empty space surrounding the ‘predator’ (or a human shepherd as in figure 1a). There is also a presence of a relatively sharp boundary of the swarm.
In this paper, we investigate a very simple particlebased model of predator–prey interactions which captures several distinct behaviours that are observed in nature. There are several wellknown mechanisms whereby the prey tries to avoid the predator. One wellstudied example is predator confusion, which occurs when the predator is ‘confused’ about which individual to pursue. Predator confusion decreases the predators' ability to hunt their prey. To quote Krause & Ruxton [3, p. 19], ‘predator confusion effect describes the reduced attacktokill ratio experienced by a predator resulting from an inability to single out and attack individual prey’. Bazazi et al. [12] studied marching insects and demonstrated that their collective behaviour functions partly as an antipredator strategy. During hunting, predators become confused when confronted with their prey swarm [13] and predator confusion was observed in 64% of the predator–prey systems studied in [14]. Predator–prey dynamics were also studied using computer models (e.g. [5,15–17]). Zheng et al. [15] studied a mathematical model of schools of fish, which demonstrates that collective evasion reduces the predator's success by confusing it. Olson et al. [16] used simulated coevolution of predators and prey to demonstrate that predator confusion gives a sufficient selective advantage for swarming prey. Similar preference for swarming in the presence of a confused predator was investigated in [5].
While there are many models in the literature that demonstrate complex predator–prey dynamics, most of these models are too complex to study except through numerical simulations. The goal of this paper is to present a minimal mathematical model which is carefully chosen so that (i) it is amenable to mathematical analysis and (ii) it captures the essential features of predator–prey interactions. A commonly used approach to swarm dynamics is to represent each prey by a particle that moves based on its interactions with other prey and its interaction with the predator. There is a large literature on particle models in biology, where they have been used to model biological aggregation in general [1,18–22] and locusts [21] or fish populations [15,23–27] in particular. This is the approach that we take in this paper as well.
We now introduce the model that we study in this paper. We assume that there are N prey whose positions j = 1 … N follow Newton's law so that Here, F_{j}_{,prey−prey} + F_{j}_{,prey−predator} is the total force acting on the jth particle, μ is the strength of ‘friction’ force and m is its mass. We make a further simplification that the mass m is negligible compared with the friction force μ. After rescaling to set μ = 1, the model is then simply so that the prey moves in the direction of the total force. This reduces the secondorder model to a firstorder model, which makes it easier to analyse mathematically. Similar reduction was used, for example, in the analysis of locust populations [28] and other biological models [19,29]. Various forms can be considered for prey–prey interactions. To keep cohesiveness of the swarm, we consider the interactions which exhibit pairwise shortrange repulsion and longrange attraction, averaged over all of the particles. For concreteness, we consider the endogenous prey–prey interaction of the form The term represents Newtoniantype shortrange repulsion that acts in the direction from x_{j} to x_{k}, whereas −a(x_{j} − x_{k}) is a linear longrange attraction in the same direction. While more general attraction–repulsion dynamics can be considered, we concentrate on this specific form because more explicit results are possible. In particular, in the absence of exogenous prey–predator force, this particular interaction has been shown to result in uniform swarms [30,31]. In general, the distribution inside the swarm can vary and have fluctuations; however, uniform density of a swarm is often a good firstorder approximation for many swarms. For example, Miller & Stephen [32] found that the flocks of sandhill cranes feeding in cultivated fields had distribution close to uniform, regardless of flock size. See [19, pp. 537–538], and references therein for further examples and discussion of prevalence of nearly uniform distribution of flocks in nature.
The prey–predator interactions are modelled in a similar fashion: again for concreteness assume that there is a single predator whose position we denote as . Assuming that the predator acts as a repulsive particle on the prey, we take with b being the strength of the repulsion. Finally, we model the predator–prey interactions as an attractive force in a similar way, . We consider the simplest scenario where F_{predator−prey} is the average over all predator–prey interactions and each individual interaction is a power law, which decays at large distances; the prey then moves in the direction of the average force. These assumptions result in the following system: 1.1and 1.2
To illustrate the results and motivate the analysis in this paper, consider the numerical simulations of the particle model (1.1) and (1.2) shown in figure 2. We use the strength c of the predator–prey attraction as the control parameter, with other parameters as given in the figure. In the second row with c = 0.4, random initial conditions for prey and predator positions are taken inside a unit square. The swarm forms a ‘ring’ of constant density with a predator at the centre of the ring. Our first result is to fully characterize this ring in the limit of large swarms; see result 2.1. Our main result characterizes the stability of this ring. In result 3.1, we show that the ring is stable whenever 2 < p < 4 and 1.3
With parameters as chosen in figure 2 this corresponds to 0.2190 < c < 0.7557. When c is decreased below 0.2910 (row 1), the ring becomes unstable and the predator is ‘expelled’ out of the ring; the swarm escapes completely. A very different instability appears if c is increased above 0.7557 (row 3). In this case, we show that the ring also becomes unstable owing to the presence of oscillatory instabilities, whereby the predator ‘oscillates’ around the ‘centre’ of the swarm. After some transients, the system settles into a ‘rotating pattern’ where the predator is continually chasing after its prey, without being able to fully catch up to it. As c is further increased (row 4), the motion becomes progressively chaotic until the predator is finally able to catch the prey (row 5).
Our approach is to take the continuumlimit N→∞ of (1.1) and (1.2), which results in the nonlocal integrodifferential equation model [19–22] 1.4 1.5 1.6
Here, ρ(x, t) denotes the density distribution of the prey swarm at position so that and v(x, t) is the swarm's velocity field. The system (1.4)–(1.6) is obtained by choosing the initial density to be where δ is the delta function. Equation (1.4) simply reflects the conservation of mass of the original prey system (1.1) (as no prey particles are created or destroyed); with the mass normalized so that ρ(x, t) represents a probability distribution. By taking different pairwise endogenous forces, the steady state to (1.1) and (1.2) with no exogenous force (b = 0) presents a wide variety of patterns [33–35]. Similar equations have been used to model animal aggregation in [21,28,36–39]. The classical Keller–Segel model for chemotaxis also contains a Newtonian intraspecies interaction [40,41]. Aggregation models also appear in material science [42–44], vortex motion [45–48] where Newtonian potential arises for vortex density evolution and granular flow [49,50].
We now summarize the paper. In §2, we construct the steadystate solution consisting of a ring of prey particles of uniform density that surround the predator at the centre. In §3, we study its stability. We conclude with some extensions of the model and discussion of some open problems in §4.
2. ‘Confused’ predator ring equilibrium state
We start by constructing the ‘ring’ steady state of the model (1.4)–(1.6), as shown in the last picture of the second row of figure 2. Consider a steady state for which the predator is at the centre of the swarm, surrounded by the prey particles. The predator is ‘trapped’ at the centre of the prey swarm while the prey forms a concentric annulus where the repulsion exerted by the predator cancels out owing to the symmetry. We state the main result as follows.
Result 2.1.
Define 2.1
The system (1.4)–(1.6) admits a steady state for which z = 0, ρ is a positive constant inside an annulus and is zero otherwise.
Figure 3a illustrates this result. For parameters as shown in the figure, the discrete model (1.1) and (1.2) generates a stable ring steady state, which is shown with dots. Solid curves show the continuum result (2.1), in excellent agreement with the discrete model (1.1) and (1.2).
The fact that the density is constant inside a swarm is a result of the careful choice of the forces in (1.1): namely, the nonlinearities are both Newtonian. The proof of result 2.1 follows closely [30,31] and uses the method of characteristics, a common technique to find steady states in the aggregation model.
Derivation of result 2.1. Define the characteristic curves X(X_{0}, t) which start from X_{0} at t = 0 2.2Using (1.4), along the characteristic curves x = X(X_{0}, t), ρ(x, t) satisfies 2.3
Note that so that from (1.5) we obtain 2.4where is conserved. Then (2.3) becomes 2.5which has a solution ρ(X(X_{0}, t), t) approaching aM/π as t → ∞ and independent of the location, as long as ρ(X_{0}, 0) > 0.
Next, we seek a steady state such that ρ is constant inside , ρ zero outside , where is an annulus with R_{1}, R_{2} and z to be determined. Using the identity 2.6and for , we compute 2.7The assumption of the steady state implies that (2.7) is zero for all , which in turn implies that z = 0, and so that and Conversely, with this choice of R_{1} and R_{2}, v = 0 whenever ρ ≠ 0. Moreover by symmetry, dz/dt = 0 so that v, ρ, z as in result 2.1 constitute a true steady state of (1.4)–(1.6).
3. Transition to chasing dynamics
As illustrated in figure 2, the ring steadystate configuration can transition to a moving configuration in two ways: if the predator strength c is sufficiently decreased, the swarm will escape the predator. If c is increased past another threshold, the predator becomes more ‘focused’ and less ‘confused’, resulting in ‘chasing dynamics’ which can lead to very complex periodic or chaotic behaviour. Similar dynamics can be observed in nature, as figure 4 illustrates. The onset of these dynamics can be understood as a transition from stability to an instability (i.e. bifurcation) of the ring steady state. The destabilizing perturbation corresponds to the translational motion of the predator as well as the inner or outer boundary of the ring.^{1}
To understand these bifurcations, we consider the perturbations of the inner boundary, outer boundary, as well as the predator itself. These perturbations are of the form 3.1 3.2 3.3where Note that this form of perturbation preserves the total mass which is an invariant of the model. In appendix A, we show that λ satisfies the eigenvalue problem 3.4where
The eigenvalues of are given by λ = 0 and which satisfy where and The eigenvalues are stable (i.e. ) if and only if B > 0 and C > 0. Note that, when c = 0, we get B = 1, C < 0 so that λ_{−} < 0 < λ_{+} and the ring is unstable. As c is increased, either λ_{+} or λ_{−} cross zero. This occurs precisely when c = c_{0}, where 3.5with C > 0 if and only if c > c_{0}. If p ≤ 2, then B > 0 for all c > c_{0} so that . If 2 < p, a Hopf bifurcation occurs when B = 0 with C > 0; i.e. when c = c_{hopf} > c_{0}, where 3.6Note 0 < c_{0} < c_{hopf} if and only if 2 < p < 4 (with c_{0} > c_{hopf} if p > 4, c_{hopf} = ∞ if p = 2 and c_{hopf} < 0 if p < 2). Therefore, if and only if one of the following holds: (i) p ≤ 2 and c > c_{0}; (ii) 2 < p < 4 and c_{0} < c < c_{hopf}.
We summarize as follows.
Result 3.1.
Consider the ring steady state of (1.4)–(1.6) given in result 2.1. Let c_{0}, c_{hopf} be as defined by (3.5) and (3.6), respectively. The ring stability with respect to translational perturbations is characterized as follows:

— If p ≤ 2: the ring is translationally stable if c_{0} < c, and unstable if c < c_{0}.

— If 2 < p < 4: the ring is translationally stable if c_{0} < c < c_{hopf}. It is unstable owing to the presence of a negative real eigenvalue if c < c_{0}. As c is increased past c_{hopf}, the ring is destabilized owing to a Hopf bifurcation.

— If p > 4: the ring is unstable for all positive c.
This analysis reveals that there are three distinct regimes, which depend on the power exponent p of the predator–prey attraction. If p < 2, then at close range the prey moves faster than the predator and can always escape. As a result, the predator can never catch the prey no matter how large c is. The most interesting regime is 2 < p < 4. As c is increased just past c_{hopf}, complex periodic or chaotic chasing dynamics result, but the predator is still unable to catch the prey. The shape of the perturbation is reflected in the actual dynamics when c is close to c_{hopf} (such as in figure 2, row 3); however, as c is further increased, nonlinear effects start to dominate and linear theory is insufficient to describe the resulting dynamics (see figure 2, rows 4 and 5). For even larger c, the predator finally ‘catches’ the prey; this is illustrated in figure 2, row 5; see §4 for further discussion of this.
Note that c_{0} = c_{hopf} when p = 4, in which case the stable band disappears. If p > 4, then c_{hopf} < c_{0} and the ring configuration is unstable for any c. In this case, the swarm escapes completely if c < c_{hopf} but chasing dynamics and catching of the prey can still be observed if c > c_{hopf}.
4. Discussion and extensions
The minimal model (1.1) and (1.2) supports a surprising variety of predator–swarm dynamics, including predator confusion, predator evasion and chasing dynamics (with rectilinear, periodic or chaotic motion).
Biologically, our model is useful in two ways. First, despite its simplicity our model has an uncanny ability to reproduce the complex shapes of a swarm in predator–prey systems. This is illustrated in figure 4. Second, the mathematical analysis of this model provides some rudimentary biological insight into general forces at play, which we now discuss.
Formula (3.6) shows that the prey–prey attraction that is responsible for prey aggregation, controlled by parameter a in (1.1) and (1.2), is detrimental to prey: c_{hopf} is a decreasing function of a so that increasing a makes it easier for the predator to catch the prey. This is also in agreement with several other studies. For example, Fertl and coworkers [51,52] observed groups of about 20–30 dolphins surrounding a school of fish and blowing bubbles underneath it in an apparent effort to keep the school from dispersing, while other members of the dolphin group swam through the resulting ball of fish to feed. In a survey [1], the authors suggest that factors other than predator avoidance, such as food gathering, ease of mating, energetic benefits or even constraints of physical environment, are responsible for prey aggregation. Our model also supports this conclusion.
The parameter b in the model (1.1) and (1.2) can be thought of as the strength of prey–predator repulsion. Formula (3.6) shows that c_{hopf} is an increasing function of b so that increasing b is beneficial to the prey.
The parameter p can be vaguely interpreted as the predator ‘sensitivity’ when the prey is close to the predator and can be thought of as a measure of how sensitive the predator is to a nearby prey. Simple calculus shows that c_{hopf} has a minimum at p = p_{optimal} given by 4.1provided that b > a (no optimal p exists otherwise with c_{hopf} → 0 as p → ∞; figure 5). From the point of view of the predator, this choice of sensitivity requires the least strength c for success. It is unclear however whether this optimal value has a true biological significance or is simply an artefact of the model.
So far, we have concentrated on the onset of chasing dynamics as c crosses c_{hopf}, as this value is computable analytically. This is a precursor to the predator catching the prey, but for values of c just above c_{hopf} the prey still escapes. Let us investigate further numerically what happens for larger values of c when the predator can actually ‘catch’ the prey. For concreteness, we say that the prey is caught if the distance between it and the predator falls below a certain kill radius, which we take to be 0.01 in our simulations (numerically, the problem becomes unstable when this distance becomes too small as the velocity of the prey and predator increases without bound). Whenever the prey is caught, we remove it from the simulation (and decrease N by 1 in (1.1) and (1.2)). Consider the parameters p = 3, a = 1, b = 0.2, c = 1.8 and suppose there are N = 200 prey initially. Figure 6 shows the number of prey as a function of time. It shows that the rate of consumption is higher with fewer individuals. The reader is invited to see the movie of these simulations.^{2}
Let c_{catch} be the smallest value of predator strength c for which the predator is able to catch the prey. We compute this value using full numerical simulations of (1.1) and (1.2) for several values of N, while fixing the other parameters to be p = 3, a = 1, b = 0.2. The results are summarized in the following table:
Note that c_{catch} is increasing with N, which is also consistent with figure 6 showing that the kill rate increases when there are fewer particles. This suggests that all else being equal, having more individuals is beneficial to prey, in that a higher predator strength c is required to catch the prey when N is increased. This may be owing to the fact that the predator becomes more ‘confused’ by the various individuals inside the swarm when there are more of them.
From a mathematical point of view, our analysis is rather nonstandard: the main result is obtained by doing a stability analysis on the entire swarm in the continuum limit, which can be thought of as an infinitedimensional dynamical system, or, alternatively, a nonlocal PDEODE system (1.4)–(1.6). Below we discuss several possible extensions of the model.
4.1. Nonuniform state
The first extension is to replace the prey–predator interaction in (1.1) by a more general power nonlinearity, for example 4.2with the equation for the predator unchanged; the original model corresponds to q = 2. As before, there is a steady state with the predator z = 0 at the centre with the swarm forming a ring around it. Unlike the q = 2 case, the density of the swarm is no longer uniform. Using a computation similar to the q = 2 case, we find that, in the continuum limit, the density is given by 4.3with R_{1}, R_{2} satisfying 4.4result 2.1 is recovered by choosing q = 2 in (4.4).
From (4.3) we note that for q < 2, the density is higher further away from the predator; conversely for q > 2 the density is higher closer to the predator. This compares favourably with full numerical simulations as shown in figure 3. However, the computation of stability for the nonconstant density state remains an open problem.
4.2. Multiple predators
It is easy to generalize (1.1) and (1.2) to include multiple predators. For example, replace (1.1) by 4.5and replace z by z_{k} in (1.2) (more complex predator–predator interactions can similarly be added). Even more complex dynamics can be observed. Multispecies interaction has been studied in several other contexts recently, including crowd dynamics and pedestrian traffic [53,54], decisionmaking in the group with strong leaders [55] and generalization of the Keller–Segel model to multispecies in chemotaxis [56–58].
Here, we briefly consider the possible steady states of the swarm in the presence of two stationary predators (i.e. c = 0). Consider two predators located symmetrically at z_{1} = d and z_{2} = −d. Figure 7 shows some of the possible steady states for various values of d. As d is decreased, the swarm splits into two. The swarm is symmetric with respect to x and yaxes but is not radially symmetric.
The solid curve in figure 7 shows the continuum limit of (4.5) which is obtained by computing the evolution of the boundary ∂D of the swarm, while assuming that swarm density is constant. Using the divergence theorem, the velocity can then be computed using only a onedimensional integration 4.6
where we assumed that the centre of mass of the swarm is at the origin, and where the area is also a onedimensional computation.
4.3. Acceleration and other effects
Introducing acceleration allows for a more realistic motion. A more general model is 4.7 4.8where μ_{j}, μ_{0} are friction coefficients of prey and predator, respectively, and m_{j}, M are their masses. Figure 8 illustrates some of the possible dynamics of these models. Even more complex models exist in the literature. For example, to obtain a more realistic motion for fish an alignment term is often included, which can lead to milling and flocking patterns even in the absence of predator [26,59,60].
Many models of collective animal behaviour found in the literature include terms such as zone of alignment, angle of vision, acceleration, etc. These terms may result in a more ‘realisticlooking’ motion, although it can be difficult in practice to actually measure precisely how ‘realistic’ it is (but see [11,61] for work in this direction). Moreover, the added complexity makes it very difficult to study the model except through numerical simulations. Our minimal model shows that these additional effects are not necessary to reproduce complex predator–prey interactions.
Funding statement
T.K. was supported by a grant from AARMS CRG in Dynamical Systems and NSERC grant no. 47050. Some of the research for this project was carried out while Y.C. and T.K. were supported by the California Research Training Program in Computational and Applied Mathematics (NSF grant DMS1045536).
Acknowledgements
We are grateful to Ryan Lukeman and Yann Arthus Bertrand for providing us with their photographs. We also thank the anonymous referees for their valuable comments, which improved the paper significantly.
Appendix A
In this appendix, we derive eigenvalue problem (3.4) for the perturbations of the form (3.1)–(3.3). Let . The velocity then becomes A1
Using (2.6) with we get A2At the steady state o_{i} = 0 and v = 0 so that (A 2) simplifies to A3
On the inner boundary, we have and linearizing we obtain
Evaluating the perpendicular component yields
We equate along the perpendicular component to finally obtain A4
The same computation along the outer boundary yields A5
Next, we linearize predator equation (1.6) around ring steady state (2.1). We estimate where h.o.t. denotes higher order terms that are quadratic in o_{i}. We then compute explicitly and approximate Linearizing predator equation (1.6) then yields A6
The three equations (A 6), (A 4) and (A 5) then yield a closed threedimensional eigenvalue problem A7
Problem (3.4) is obtained by substituting (2.1) into (A 7).
Footnotes
↵1 As discussed in the derivation of result 2.1, for large time, the density ρ(x,t) rapidly approaches a constant on its support, ρ → aM/π and the equation for ρ along characteristics is independent of the boundary shape or the form of predator–prey interactions (parameters c and p in (1.2)). As such, tracking the evolution of the boundary and the predator is sufficient to determine the stability of the ring state.
↵2 We created a website which contains the movies showing the simulations of predator–swarm interactions from this paper. These can be viewed by following the link: http://goo.gl/BC6pyC.
 Received December 29, 2013.
 Accepted February 6, 2014.
 © 2014 The Author(s) Published by the Royal Society. All rights reserved.