## Abstract

Predators often have type II functional responses and live in environments where their life history traits as well as those of their prey vary from patch to patch. To understand how spatial heterogeneity and predator handling times influence the coevolution of patch preferences and ecological stability, we perform an ecological and evolutionary analysis of a Nicholson–Bailey type model. We prove that coevolutionarily stable prey and searching predators prefer patches that in isolation support higher prey and searching predator densities, respectively. Using this fact, we determine how environmental variation and predator handling times influence the spatial patterns of patch preferences, population abundances and per-capita predation rates. In particular, long predator handling times are shown to result in the coevolution of predator and prey aggregation. An analytic expression characterizing ecological stability of the coevolved populations is derived. This expression implies that contrary to traditional theoretical expectations, predator handling time can stabilize predator–prey interactions through its coevolutionary influence on patch preferences. These results are shown to have important implications for classical biological control.

## 1. Introduction

Most species live in an inherently heterogeneous world in which their growth rates vary spatially. For prey species, the variability in growth and reproduction may stem from the variation in the nutritional quality of their resources, climate and abundance of predators. For predator species, variability in growth and reproduction may reflect variation in the abundance, size or quality of their prey. Understanding the coevolution of predators and prey patch selection strategies in heterogeneous environments and the population dynamic implications of this coevolution have been the focus of several theoretical endeavours (Holt 1984; van Baalen & Sabelis 1993, 1999; Hochberg & Holt 1995; Schreiber *et al*. 2000, 2002; Cressman *et al*. 2004). Under the assumption that only the life history traits of the prey vary among patches, van Baalen & Sabelis (1993) found that coevolution of patch selection inevitably results in congruent patch choices, where predators and prey prefer the same patches, and promotes ecological stability under restricted circumstances. However, if life history traits of the predator as well as the prey vary amongst patches, then coevolved predator–prey populations can exhibit contrary patch choices with prey and predators preferentially selecting different patches. Moreover, these contrary choices can have a strong stabilizing influence on predator–prey interactions (Schreiber *et al*. 2000, 2002).

In the aforementioned analyses, the underlying assumption is that the rate at which predators encounter and attack their prey is limited by their ability to find the prey. However, at higher prey densities, one would expect that the encounter rate of predators with prey saturates due to handling time, i.e. the time between encountering a prey and searching for the next prey item. Handling time leads to a type II functional response and has a destabilizing effect on predator–prey dynamics (Holling 1959; Oaten & Murdoch 1975). For instance, adding handling time to the classical Lotka–Volterra predator–prey model switches a neutrally stable equilibrium to a globally unstable equilibrium. Alternatively, for Nicholson–Bailey predator–prey dynamics, there is a stability trade-off between the degree of aggregation of predator attack and the predator handling time. Higher levels of predator aggregation (i.e. *k*<1 in the negative binomial escape function, where *k* is the negative binomial shape parameter; May 1978) require shorter predator handling times to stabilize predator–prey dynamics (Hochberg & Holt 1995; Getz & Mills 1996). The prevalence of the type II functional response in the empirical literature resulted in Getz and Mills questioning ‘the generality of a purely search-limited encounter rate in population models that underlie more sophisticated behavioural and evolutionary analysis’.

To address this concern, we examine the effect of handling time on the coevolution of predator and prey patch preferences and the implications for ecological stability. Our analysis yields explicit expressions for the patch preferences and criteria for ecological stability of coevolutionary stable populations. Using these results we examine how environmental variation and predator handling time influence spatial patterns of per-capita predation rates and population abundance, under what conditions predator handling time coupled with coevolution of patch preferences can stabilize predator–prey interactions, and discuss the implications for classical biological control.

## 2. The model

We consider a system consisting of a population of predator and prey that disperse among *n* patches. To model this system, we generalize a model of Hassell & May (1973) which assumes that prey and predator populations have discrete and synchronized generations. Every generation the predators and prey select patches according to fixed behavioural strategies. The prey lay a fixed proportion *α*_{i} of their eggs in patch *i*, where *i*=1, …, *n* and . The predators spend a fixed proportion *β*_{i} of their searching time in patch *i*, where *i*=1, …, *n* and . We assume that the probability of a prey in patch *i* escaping predation in generation *t* is exp(−*a*_{i}*β*_{i}*S*_{t}), where *a*_{i} is the searching efficiency of the predator in patch *i* and *S*_{t} is the total number of searching predators in the environment in generation *t*. We assume that intrinsic rates of growth *λ*_{i} of a prey in patch *i* are greater than one and that the expected numbers *θ*_{i} of predators produced per attacked prey in patch *i* are positive. If *N*_{t} and *P*_{t} denote the total number of prey and predators, respectively, in generation *t*, then we arrive at the following *n*-patch model(2.1)To complete the formulation of this model, we need to specify *S*_{t} in terms of *N*_{t} and *P*_{t}. In the electronic supplementary material, appendix A, we show that if *b*_{i} is the mean handling time in patch *i*, then(2.2)Equations (2.1) and (2.2) generalize a model of Hassell (2000) by allowing all of the traits of the predator and prey to vary amongst patches and generalize a model of Holt & Hassell (1993) by including predator handling times. In the special case of one patch, equations (2.1) and (2.2) correspond to a Nicholson–Bailey model with a type II functional response (Getz & Mills 1996).

## 3. Coevolutionary stable strategies

To understand how the predator and the prey may coevolve their patterns of patch selection, we employ the dynamical theory of evolutionary stable strategies. Following this approach, we consider a population of predator and prey, the *resident population*, that play the patch selection strategies *α*=(*α*_{1}, …, *α*_{n}) and *β*=(*β*_{1}, …, *β*_{n}). Assume a subpopulation of *mutant prey* with density *M* that play the mutant patch selection strategy . If the resident and mutant population only differ in their patch selection strategy and have identical life history traits otherwise, then the resident–mutant dynamics arewhereIf the resident population is at the equilibrium , then the *mutant prey invasion rate* is given bywhereis the density of searching resident predators at equilibrium. Alternatively, if a subpopulation of *mutant predators* with density *Q* plays the mutant patch selection strategy , then the resident–mutant population dynamics becomewhereandare the densities of searching resident and searching mutant predators, respectively. If the resident population is at the equilibrium , then the *mutant predator invasion rate* is given bywhere is the density of searching resident predators at equilibrium. If (respectively, ), then the mutant prey (respectively, predator) population can invade, otherwise the mutant invasion fails. Unlike *I*_{N}, which is linear in , *I*_{P} is nonlinear in .

Following Cressman and co-workers (Cressman 2003; Cressman *et al*. 2004), a *coevolutionarily stable strategy* (coESS) is a strategy that, when adopted by most of the populations, resists invasion attempts by mutants playing any other strategy and is stable with respect to spatial perturbations. In other words,(3.1)for all , and at least one of the inequalities(3.2)is true for each perturbed distribution , close to (*α*, *β*). In the electronic supplementary material, appendix B proves that there is a unique strategy (*α*, *β*) satisfying equation (3.1). Moreover, numerical simulations and analysis of the two-patch case suggest that this strategy satisfies equation (3.2). This unique coESS is given by(3.3)where correspond to the equilibrium densities attained by the prey and the searching predators when restricted to patch *i*. More precisely,(3.4)From these expressions, it follows that the coESS corresponds to the prey preferring patches that intrinsically support higher prey numbers and searching predators prefer patches that intrinsically support higher searching predator numbers. When the predator and prey are playing the coESS, the equilibrium population density is given by(3.5)where(3.6)is the equilibrium density of the predator if both species were confined to patch *i*.

Since *S*_{i} is independent of *b*_{i}, the handling time of the predator does not influence the coevolved patch preferences of the searching predators. However, handling time does influence the patch preferences of the prey as well as the coevolved distribution and abundance of the predator and the prey. To understand the role of predator handling time on these quantities, we assume that the handling time and searching efficiency of the predator do not vary from patch to patch (i.e. *b*_{1}=⋯=*b*_{n} for some *b* and *a*_{1}=⋯=*a*_{n}=*a* for some *a*). Moreover, without any loss of generality, we assume that the coevolved prey prefer the lower numbered patches (i.e. *α*_{1}>*α*_{2}>⋯>*α*_{n}) when *b*=0. In the electronic supplementary material, appendix B, we show is an increasing, convex function of *b* that has a vertical asymptote at *b*=*θ*_{i}(1−1/*λ*_{i})ln *λ*_{i} for all *i*=1,2, …, *n*−1. Since is increasing with *b*, predators with longer handling times result in the prey exhibiting greater preferences for the lower numbered patches. For predators whose handling times are close to the critical value *θ*_{1}(1−1/*λ*_{1})ln*λ*_{1}, prey preferences and densities are highly skewed toward patch one (figure 1). Equation (3.6) implies that the aggregation of prey in the lower numbered patches results in the predator densities accumulating primarily in these patches despite the fact that searching predators may not prefer these patches.

To extract some additional information about the coESS and how it depends on environmental heterogeneity, we consider two types of environmental gradients. First, as considered by van Baalen & Sabelis (1993), suppose that only the prey intrinsic rate of growth varies among the patches, i.e. *λ*_{1}>*λ*_{2}>⋯>*λ*_{n} and *θ*_{1}=⋯=*θ*_{n}. Then for all predator handling times, coevolution yields prey and predator patch preferences and densities decreasing along the environmental gradient, i.e. , , and . Moreover, since the per-capita predation rates decrease with prey density along the environmental gradient, the per-capita predation rates exhibit positive density dependence for all *b*≥0 (Hassell & May 1973, 1974; Pacala *et al*. 1990). However, as longer predator handling times result in the prey aggregating primarily in the first patch, longer handling times weaken the positive correlation between per-capita predation rates and prey density.

Second, suppose the differences in patches form a gradient for the prey and the predator, i.e. *λ*_{1}>*λ*_{2}>⋯>*λ*_{n} and *θ*_{i} is proportional to *λ*_{i}. In the electronic supplementary material, appendix B, we show that the prey densities and patch preferences increase along the environmental gradient while searching predator preferences decrease along the environmental gradient, i.e. and . Since the per-capita predation rates increase along the environmental gradient, the per-capita predation rates exhibit inverse density dependence (Hassell 1984; Chesson & Murdoch 1986; Pacala *et al*. 1990). Since long predator handling times result in the prey aggregating into lower numbered patches without affecting the per-capita predation rates, long predator handling times weaken the negative correlation between per-capita predation rates and prey density. Predator handling times play an important role in determining the pattern of predator abundance along the environmental gradient. For short handling times, the combined abundance of searching and handling predators decreases along the environmental gradient, i.e. . Alternatively, for long handling times, the combined abundance of searching and handling predators increases along the environmental gradient, i.e. . This latter pattern occurs as searching predators handle prey with greater likelihood in lower numbered patches.

## 4. Ecological stability

The coESS analysis was based on the assumption that the resident populations are at an equilibrium of equations (2.1) and (2.2). Therefore, to determine when the coevolutionary analysis is relevant, we need to determine under what conditions this equilibrium is stable. In the electronic supplementary material, appendix C proves that the coESS equilibrium is stable if, and only if,(4.1)This stability criterion implies that by promoting the coevolution of aggregation, predator handling time can stabilize predator–prey interactions (figure 2). To understand under what conditions predator handling time promotes ecological stability, we consider two scenarios involving an environment consisting of low and high-quality patches. The first scenario assumes that the prey performs poorly in low-quality patches and the predators perform equally well in all patches. Thus, the coevolved populations exhibit congruent choices and tend to aggregate in the high-quality patches. The second scenario assumes that both the predator and prey perform poorly in the low-quality patches. Thus, the coevolved populations exhibit contrary choices and tend to aggregate in the low-quality patches. Figure 3 illustrates that predators with longer handling times promote stability provided that there is a sufficiently small percentage of patches in which the populations aggregate. For congruent choices, a very small percentage of high-quality patches are necessary for stability (figure 3*a*). Alternatively, when the populations exhibit contrary choices, predator handling times can stabilize interactions even if more than 50% of the patches are low quality (figure 3*b*). Figure 3 also illustrates that the range of percentages of high-quality patches resulting in stability increases with predator handling time for populations exhibiting contrary choices but decreases for populations exhibiting congruent choices.

## 5. Discussion

The theory of coevolution proposes that the interactions between predator and prey undergo reciprocal evolution due to natural selection (van Baalen & Sabelis 1993, 1999; Thompson 1999; Abrams 2000; Schreiber *et al*. 2000, 2002; Cressman *et al*. 2004). In patchy landscapes where there is variation in environmental factors influencing fitness, reciprocal evolution of patch preferences can at first glance yield a paradox. Predators should prefer patches where prey are common and should avoid patches where predators are common. While each species may choose strategies that thwart the other, the coevolutionary resolution to this paradox may be not be obvious (Lima 2002). A powerful methodology to resolve this paradox is to combine the concepts of evolutionary game theory with the tools on nonlinear ecological dynamics.

Our coevolutionary analysis reveals that depending on the nature of environmental heterogeneity and the length of predator handling times, the coevolved populations exhibit different spatial patterns of abundance, preferences and per-capita predation rates. First and foremost, predators with long handling times promote the coevolution of aggregation. The coevolved prey exhibit strong preferences for patches that in isolation would support highest prey densities. By aggregating to the preferred patches, prey dilute predation risk and form a ‘selfish herd’ (Hamilton 1971; Rands *et al*. 2004). Moreover, independent of predator searching preferences, predator densities are highest in these patches where predators are most likely to be handling prey. Second, when there is heterogeneity in prey life history traits but no variation in predator life history traits, the patterns of abundance, patch preferences and per-capita predation rates are congruent. Mirroring previous work on related models (Holt 1984; van Baalen & Sabelis 1993, 1999; Schreiber *et al*. 2000, 2002; Cressman *et al*. 2004), we find that prey and searching predators prefer patches that in isolation would support higher prey densities, the abundance of the prey and predator are greatest in these preferred patches, and per-capita predation rates correlate positively with prey density (Pacala *et al*. 1990). The independence of these patterns from handling time is somewhat surprising in light of the ecological expectation that long handling times should yield inverse density dependence in per-capita predation rates (Hassell 1984, 2000).

When the life history traits of the prey and predator covary, we find that the patch preferences of the prey and the searching predators are contrary: prey prefer patches that the searching predators tend to avoid. Moreover, these patch preferences always result in per-capita predation rates being negatively correlated with prey density. These results are consistent with prior work (Schreiber *et al*. 2000, 2002) and occur as the predators are searching for high-quality prey and the prey are shifting to ‘enemy-free space’ (Jeffries & Lawton 1984; Holt & Lawton 1993). These contrary patch preferences have been observed experimentally in a system involving the moth-specialist solitary parasitoid (Fox & Eisenbach 1992). Unlike per-capita predation rates, the correlation between predator densities and prey densities depends critically on the length of the predator handling time. For shorter handling times, predator densities correlate negatively with the prey densities. For longer handling times, predator densities correlate positively with prey densities. Hence, for longer handling times, predator densities and per-capita predation rates are inversely related. This inverse relationship has been observed in studies of the egg parasitoid *Trichogramma evanescens* and its host, the stored product moth *Plodia interpunctella* (Hassell 1982) in which ‘parasitism (rates are) inversely density dependent, despite the tendency for the parasitoids to aggregate on the leaves with highest host densities’ (Hassell 2000, p. 547) .

Placing predator–prey dynamics in a coevolutionary context can alter our perception of mechanisms for ecological stability. Throughout the ecological literature, a type II functional response has been shown to destabilize predator–prey interactions (Rosenzweig 1971; Oaten & Murdoch 1975; Getz & Mills 1996). For example, the inclusion of a type II functional response into a Lotka–Volterra predator–prey model transforms a neutrally stable equilibrium to a globally unstable equilibrium where oscillations of ever increasing amplitudes result in the demise of one or both species. Alternatively, Getz & Mills (1996) found that long predator handling times disrupt the stabilizing effect of aggregation in Nicholson–Bailey models with a negative binomial escape function. In stark contrast to these ecologically derived predictions, our analysis shows that predator handling time can stabilize predator–prey interactions by altering the coevolutionary patch preferences of the prey. Intuitively, this stabilizing effect stems from handling times promoting the coevolution of aggregation. Aggregation has been promoted as a stabilizing factor in many theoretical studies of predator–prey dynamics (Hassell & May 1974; Chesson & Murdoch 1986; Murdoch & Stewart-Oaten 1989; Holt & Hassell 1993). Hence, through its coevolutionary influence on patch preferences, handling time creates a tension between two ecological forces: the stabilizing force of aggregation and the destabilizing force of predator saturation. We have illustrated two environments that result in the stabilizing force overcoming the destabilizing force: environments that promote contrary choices and that contain sufficiently many high-quality patches, and environments that promote congruent choices and that contain sufficiently few high-quality patches.

Our results are based on Nicholson–Bailey models that assume the populations have discrete generations. While there have been no studies of how handling time and coevolution effect ecological stability in continuous-time models, previous studies (Křivan 1997; Cressman *et al*. 2004) have considered the effect of ideal free dynamics for two patch models with classical Lotka–Volterra predator–prey interactions. For these models, ideal free distributions are coevolutionarily stable (Cressman *et al*. 2004) and stabilize the ecological dynamics by creating a global attractor consisting of a neutrally stable equilibrium surrounded by periodic motions (Křivan 1997). Despite the qualitative similarity to some of our results, ideal free behaviour in these continuous-time models has no effect on ecological stability at equilibrium. One possible explanation for this discrepancy is that the individuals in these continuous time models are ‘freer’ than the corresponding individuals in the discrete-time models. In the continuous-time models individuals move instantaneously from one patch to another patch, while in discrete-time models individuals require a generation to move between patches.

Our analyses have several potentially important implications for classical biological control in which predators are used to control agricultural pests. Successful biological control is often equated with a stable equilibrium at which the prey density is below an economic threshold. Key predatory traits that promote low prey abundance at equilibrium include high searching and conversion efficiencies. For continuous time models in which the predator has a type II functional response, high values of these traits destabilize predator–prey interactions. This trade-off between low prey equilibrium densities and stability for continuous time models yields ‘the paradox of biological control’ (Arditi & Berryman 1991). For coevolved populations in our model, this paradox need not arise. For example, when the predator searching efficiency does not vary between patches (i.e. *a*_{1}=*a*_{2}=…=*a*_{k}), the stability of the coevolved populations (cf. equation (4.1)) is independent of the predator searching efficiency. Hence, if environmental variation in other life history traits of the predator and prey promote the coevolution of ecological stability, then predators with high-searching efficiencies yield successful biological control in these environments. One particularly intriguing prediction concerns the stabilizing effect of marginal patches: the introduction of alternative host plants of marginal quality can stabilize otherwise unstable predator interactions without raising the equilibrium prey densities on the higher quality plants. This stabilizing effect may only be observed when there has been a sufficiently long period for coevolution to occur. In particular, for novel predator–prey associations or novel environments, unstable dynamics may be observed before an ecologically stable coESS is reached.

Biological control often involves the use of parasitoids whose ability to regulate prey may be limited by their egg load (Heimpel & Rosenheim 1998; Lane *et al*. 1999; Jervis *et al*. 2001). Previous theoretical studies (Getz & Mills 1996) suggest that parasitoid egg limitation can adversely effect biological control in two ways: egg limitation tends to destabilize host–parasitoid interactions and simultaneously raise the host equilibrium. Using the Biocat database, Lane *et al*. (1999) tested this theoretical prediction. They found that the hypothesis of egg limitation being negatively correlated with successful biological control was not supported. Lane *et al*. proposed several explanations including an evolutionary trade-off between parasitoid egg limitation and searching efficiency. As our model for one patch is mathematically equivalent to these earlier models, our analysis offers another possible explanation: egg limitation (i.e. handling time) via the coevolution aggregation can stabilize predator–prey interactions and thereby promote successful biological control.

## Acknowledgments

The authors thank Priyanga Amarasekare, Wayne Getz, Bob Holt, Volker Rudolf and two anonymous referees for encouraging comments and useful suggestions regarding this work. This research was supported in part by US National Science Foundation Grants DMS-0071994 and EF-0436318.

## Footnotes

The electronic supplementary material is available at http://dx.doi.org/10.1098/rspb.2005.3236 or via http://www.journals.royalsoc.ac.uk.

- Received April 13, 2005.
- Accepted June 25, 2005.

- © 2005 The Royal Society