Predator–prey, host–parasitoid and plant–herbivore models

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 245

Basic structure of predator–prey models

The most elementary predator–prey model is the Lotka–Volterra model:

Here, N is the population size of the prey and P is that of the predator. The prey are assumed to have an intrinsic growth rate a per unit of time, whereas the predator population decays at a rate b if there are no prey present. The death rate of each prey individual from predation is assumed to be αP, and the per capita growth rate of the predator population is assumed to be linearly dependent on prey density via the proportionality constant β. This model has the virtue of extreme simplicity, with only four parameters. The model is neutrally stable, meaning that it oscillates indefinitely with an amplitude determined solely by the starting conditions. Largely because of this property, described as ‘pathological’ by May (1976), the model is of little practical use in describing real populations, and has never been parameterized successfully for any particular interaction.

The following is the simplest plausible predator–prey model:

Here, rN is the intrinsic rate of increase of the prey population, f (N,P) is the predator functional response, K is the carrying capacity of the prey population, and g(N,P) is the numerical response of predator population.

Functional responses

Functional responses are generally classified as belonging to one of three types (Holling, 1959), as shown in Fig. 9.1. A type I functional response is a straight line, with prey consumption per predator increasing linearly with prey density. If predator and prey encounter each other at random (the ‘blundering idiot’ searching strategy) a type I response will occur. One parameter only is necessary to describe the response. Type I functional responses may occur in filter feeders, and are often assumed in host–parasite models. A type II functional response is far more common in host–parasitoid or predator–prey

248 C H A P T E R 9

where rP is the intrinsic rate of increase of the predator population, and κ is a ratio for converting prey numbers into the number of predators that they can support.

The model based on substituting eqn (9.6) into eqn (9.3) and eqn (9.4) is sometimes known as the Rosenzweig–MacArthur (Rosenzweig & MacArthur, 1963) model, and has been used quite extensively for building parameterized models (e.g. Hanski & Korpimäki, 1995).

Total responses

Some authors, particularly those studying mammals, work with the death rate from predation per prey individual, as a function of prey density (e.g. Messier, 1994; Pech et al., 1995). This is called the total response. Using a total response means that the predator dynamics are not considered explicitly. If predator dynamics occur on a much faster time scale than those of the prey (an unlikely scenario, unless the numerical response occurs through highly mobile predators aggregating on prey patches), then the shape of the total response is a product of the functional and numerical response. Alternatively, if the numerical response is slow or absent, the total response has the same shape as the functional response.

Ratio-dependent predation

A recent, and controversial, suggestion has been that predator–prey interactions are best understood in terms of the ratio of predators to prey, rather than in terms of the prey numbers themselves. The basic argument (see Berryman, 1992) is that it is unreasonable to assume that the predation rate per predator is a function only of prey density, as is assumed by the functional responses in Table 9.1. Rather, at high predator densities, there should be fewer prey per predator, and the attack rate should go down. Further, it is argued that the Holling functional response measures the attack rate of predators on prey on a very short behavioural time scale, but that eqns (9.3) and (9.4) represent population behaviour on the longer population dynamic time scale. The consequence of these assumptions, according to the proponents of the ratiodependent view, is that the predator isocline is vertical (see Fig. 9.2), meaning that predator populations will decrease below a certain threshold prey density, and increase above it, irrespective of the nature of the intrinsic density dependence in the prey species. In practice, this leads to the ‘paradox of enrichment’, in which increasing the food supply of the prey does not lead to any increase whatsoever in the prey population, but merely to an increase in the predator population. Further, it means that it is very difficult to generate a model in which the prey are regulated to a very low level by the predators,

Fig. 9.2 Graphical analysis of some common predator–prey models. In each of the figures, the zero net growth isocline for the prey N is shown as a solid line, and the zero net growth isocline for predators P is shown dashed. The arrows show the trajectory of the system around the equilibrium point, where the isoclines intersect. In all cases, the trajectories show a tendency to cycle. It can be shown that the equilibrium point is stable if the predator isocline intersects a descending part of the prey isocline, but is unstable if the intersection is on an ascending part of the isocline (see, for example, McCallum, 1987). These isoclines are very easy to construct and explore on a spreadsheet. (a) The standard Lotka–Volterra model (eqns (9.1) and (9.2)). Parameters are a = 1.0, α = 0.02, b = 0.05,

β = 0.0005. (b) Logistic prey density dependence (eqn (9.3)), a Holling type II functional response, and a numerical response of the form of eqn (9.5). Parameters are: rN = 1.0,

K = 250, α = 0.02, β = 2.0, b = 0.2, c = 0.2. The position of the vertical predator isocline does not depend on either K or rN. This means that increasing the prey carrying capacity or growth rate does not affect the equilibrium prey population size. (c) Logistic prey density dependence (eqn (9.3)), a Holling type II functional response, and a numerical response of the form of eqn (9.6). Parameters are: rN = 1.0, K = 250, α = 0.02, β = 2.0, κ = 0.3 (The position of the isoclines does not depend on rP.) (d) Logistic prey density dependence

(eqn (9.3)), a ratio-dependent type II functional response (eqn (9.7)), and a numerical response of the form of eqn (9.5). Parameters are: rN = 1.0, K = 250, m = 4.0, w = 2.4, b = 0.5, c = 0.2.

250 C H A P T E R 9

although this does in fact appear to be the outcome of some successful biological control programs.

The simplest ratio-dependent functional response can be obtained by substituting the ratio N/P for N in the Holling type II functional response (Berryman, 1992):

where m and w are constants. A disadvantage of eqn (9.7) is that it does not rest on a straightforward biological argument in the same way as the Holling type II response does.

Arditi and Saïah (1992) suggested a more flexible predator-dependent functional response:

where α and β have the same meaning as in the usual type II functional response in Table 9.1, and m is a parameter measuring the extent of dependence on predator population size, ranging from m = 0, corresponding to the usual prey-dependent model, to m = 1, corresponding to a ratio-dependent model.

The arguments in favour of ratio dependence have been trenchantly rejected by some authors, particularly Abrams (1994) – but see also the responses by Berryman et al. (1995) and Akçakaya et al. (1995). The problem with the vertical predator isocline only occurs if the assumption is made that the numerical response of the predator is some constant times the functional response. There is no reason why this should necessarily be the case, as there is certainly a difference in time scale between functional and numerical responses. A numerical response of the form of eqn (9.6) does not result in a vertical isocline, although Berryman (1992) argues that, as it involves the ratio P/N, it is, in fact, a ratio-dependent model. The time scale problem, as applied to the prey dynamics, is spurious, as Abrams points out. A differential equation model is concerned with instantaneous rates, and the functional response enters into the prey equation as a component of the instantaneous death rate of the prey. The translation to the longer population dynamic time scale is accomplished by integration of the equations. The incompatibility of the behavioural and population dynamic time scales occurs in the predator dynamics, and then only if the numerical response is assumed to be proportional to the functional response. Finally, it is certainly the case that high predator densities may influence the searching efficiency of individual predators. Prey depletion is taken care of by the prey equation, but mutual interference, or alternatively group hunting, may cause the searching efficiency of predators to be a function of their density. Another criticism of ratio-

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 251

dependent models is that they assume that the amount of prey available per predator increases to infinity as predator density declines to zero (Gleeson, 1994).

Ultimately, whether ratio dependence is a more meaningful base assumption from which to build predator–prey models than prey dependence needs to be decided on the basis of usefulness: what is the relative performance of the two approaches in practice? Models of the form of eqns (9.3) and (9.4) are very abstract indeed. In particular, they omit spatial heterogeneity. Appeals to the logic of model construction in this very simple model will not necessarily ensure that the optimal model is obtained to approximate the behaviour of more complex systems. The issues of spatial heterogeneity and spatial scale are dealt with later in the chapter.

Whether ratio-dependent models are more often appropriate in practice than prey-dependent models remains an open question. Proponents of each approach can cite examples in which their favoured approach has worked better than the alternative. My advice is to stick to prey-dependent functional responses, unless there is strong empirical evidence to the contrary, because the parameters have a simpler mechanistic interpretation than do their ratiodependent counterparts. However, as predator reproduction and prey consumption operate on very different time scales, a numerical response of the form of eqn (9.6) should generally be used in preference to one that links explicitly to the functional response such as eqn (9.5).

Host–parasitoid models

Parasitoids are insects that lay their eggs in the immature developmental stages of other insects. Successful development of the parasitoids then inevitably kills the host into which they have been laid. More attention has probably been directed at the modelling of host–parasitoid interactions than at any other form of interspecific interaction. There are several reasons for this interest. First, host–parasitoid interactions are extremely common – there are about 68 000 described species of parasitoid, but there may be up to 2 million species in total (Godfray, 1994). Second, parasitoids are often significant natural enemies of pest insect species, and are often investigated as agents of biological control. Finally, the biological properties of host–parasitoid interactions make them relatively straightforward to model. Generations of host and parasitoid usually run in step, there is often a one-to-one replacement of the developing host by the parasitoid, and, in temperate regions at least, there is often one single, non-overlapping generation of parasitoid and host each year, lending the system to a simple difference equation (see Chapter 1) modelling approach.

A major objective in modelling host–parasitoid interactions is to develop models to predict the likely success of biological control (May & Hassell,

252 C H A P T E R 9

1988; Godfray & Waage, 1991), although the extent to which many models have performed this task successfully is debated (Murdoch et al., 1984; 1985; Murdoch & Briggs 1996).

Hassell (1978) is the classic reference on host–parasitoid models. For a recent review see Mills and Getz (1996). The simplest possible host–parasitoid model is the Nicholson–Bailey model:

Here Nt is the number of hosts in generation t, and Pt is the number of parasitoids. The number of hosts in the subsequent generation is the overall reproductive rate λ multiplied by the current population size and the proportion of hosts that escape parasitism, exp(−aPt). The number of parasitoids in the next generation is simply equal to the number of hosts that do not escape times q, the number of female parasitoids emerging per parasitized host. The exponential term for the proportion of hosts attacked is obtained from two straightforward assumptions: first, that the number of encounters that parasitoids have with hosts is proportional to the densities of the two multiplied together; and second, that these encounters are randomly distributed amongst hosts. Together, these assumptions mean that the probability of a host escaping from parasitism is the zero term of a Poisson distribution with mean aPt. The model is thus assuming a type I functional response, with the parameter a measuring the search efficiency of the parasitoids.

As does the Lotka–Volterra model of predation (to which it is closely related) the Nicholson–Bailey model generates unrealistic dynamic behaviour (it is unstable, leading to diverging oscillations) and has rarely been fitted to a real system. An exception is the laboratory interaction between a greenhouse whitefly and a chalcid wasp described by Burnett (1958), which appears in a number of textbooks (for example, Hassell, 1978; Begon et al., 1996b).

More realistic host–parasitoid models are discussed in detail by May (1978), Hassell (1978), Hassell and Godfray (1992) and Mills and Getz (1996). They may involve substituting a type II or type III functional response for the aPt term in eqn (9.9) (using one of the functional forms in Table 9.1); adding density dependence to the prey population (using a functional form from Table 6.1); nonrandom parasitoid searching; or including age or spatial structure in the host population.

Aggregation and heterogeneity in predator–prey and host–parasitoid models

It is very well established that heterogeneity in parasitoid or predator attack rates has major consequences for the dynamics of the corresponding interactions

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 253

– see Murdoch and Oaten (1975), or Mills and Getz (1996) for a recent review. Heterogeneity can come from a wide variety of sources. The simplest idea is that of spatial refuges: areas of the environment in which the prey are protected from the predators. Without going into the details of the models, the existence of such refuges stops the predators reducing the prey to low levels, which would in turn lead to very low predator populations, and hence exaggerated coupled cycles in abundance. A similar ‘refuge effect’ may also result from the existence of a prey age class that is invulnerable to attack (Murdoch et al., 1987).

In general, spatial variation in the attack rate of parasitoids, without explicit refuges, also leads to stability, whether or not the attack rate of the parasitoids depends on host density (Pacala et al., 1990; Hassell et al., 1991b; Pacala & Hassell, 1991; Godfray & Pacala, 1992). Using a model based on differential equations, Murdoch and Stewart-Oaten (1989) suggest that aggregation will not necessarily lead to stability if generations are overlapping, as is typical for predator–prey interactions. The situation is further complicated if functional responses are of type II or III rather than type I (Ives, 1992).

Given the importance of aggregation to the qualitative dynamics of host–parasitoid models in particular, it is unlikely that a parameterized model that does not take aggregation into account will describe an interaction adequately. From the point of view of parameter estimation, however, a major problem is that many of the models that have been developed are very abstract. The parameters they contain are compound entities that cannot easily be reduced to items that can reliably be estimated in the field or laboratory.

Pacala et al. (1990) have shown that, for a wide variety of host–parasitoid models, the interaction is stable provided the coefficient of variation squared (σ 2/μ2) of the parasitoid attack rate in the vicinity of individual hosts exceeds 1. This ‘CV 2 > 1 rule’ suggests that considerable progress is possible in predicting the outcome of a host–parasitoid interaction, provided that σ 2/μ2 can be estimated from (s2/y−2). They modified eqns (9.9) and (9.10) by replacing the proportion of hosts that escape parasitism, exp(−aPt), by a term exp(−ap), where a is the searching efficiency of the parasitoids within a patch and p is the local (within-patch) density of searching parasitoids. The quantity that needs to be estimated is the coefficient of variation of p. The estimation process is unfortunately not entirely straightforward, as data on the searching behaviour of individual parasitoids are rarely obtainable. Pacala and Hassell (1991) provide a maximum likelihood estimation procedure to obtain the CV 2 and related data from the frequency distribution of the proportion of parasitized hosts per patch, and host density per patch. These are data that are often relatively easy to obtain. Jones et al. (1993) describe an application of the approach in practice. These papers should be consulted for further details.

254 C H A P T E R 9

A second approach to handling spatial heterogeneity is to build a model that considers the distribution in space of victims and exploiters explicitly, and models interactions on a local scale. This approach is demanding both of data and computing power, but examples are beginning to appear in the ecological literature (e.g. DeAngelis et al., 1998).

Experimental approaches

There are many examples of experiments measuring functional responses in experimental arenas, using different densities of hosts or prey, and measuring the number of individuals attacked over a certain time interval (Hassell, 1978; Fan & Petitt, 1994). In fitting functional responses there are various statistical considerations that have caused some debate, particularly in the entomological literature. One common problem is that such experiments are typically run over time periods long enough that the number of prey individuals consumed by the end of the experiment is a substantial proportion of those initially present, and prey are not replaced during the experiment. If this is the case, the functional forms in Table 9.1 cannot be used to describe the data directly. Instead, the functional response must be integrated over the course of the experiment. The number of prey still present at time t, N(t), will be a solution of the equation

where f(N,P) is the appropriate functional response from Table 9.1. For most functional responses, the solution to eqn (9.11) does not have a simple closed form, substantially complicating the estimation process. For example, the integrated form of the Holling disc equation (no. 2 in Table 9.1) can be written as

where Ne is the number of prey eaten during the course of the experiment, N0 is the number of prey present initially, and T is the duration of the experiment. Equation (9.12) is an implicit equation, meaning that it cannot be rearranged so that Ne appears on the left-hand side only. This equation is known as the random predator equation (Royama, 1971; Rogers, 1972). The problems in estimating its parameters from data have been studied in some detail by Juliano and Williams (1987). In particular, they compared a regression approach using a linearization of eqn (9.12) developed by Rogers (1972),

where Y = ln((N0 − Ne)/N0), with a nonlinear least squares approach first developed by Cock (1977). The linearization approach is easier in practice, but

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 255

it has theoretical problems, in that Ne is present on both sides of the equation. This will cause standard regression approaches to be unreliable. Using simulated data where the parameters were known, Juliano and Williams (1987) also found that it returned very biased estimates for the parameters.

The nonlinear approach requires numerical solutions of eqn (9.12), with values of the parameters α and β being selected that minimize the residual sum of squares. This process can be accomplished fairly easily using the nonlinear regression procedures built into most modern statistics packages. In general, Juliano and Williams (1987) found that the nonlinear approach performed well on their simulated data.

More recently, Fan and Petitt (1994) used simulation to compare nonlinear least squares parameter estimation for the simple Holling disc equation with linearized versions. Their conclusion, in contrast to Juliano and Williams (1987), was that the linearization

performed better than a nonlinear fit, particularly if errors were multiplicative. As Williams and Juliano (1996) point out, this conclusion is valid only if the effects of prey depletion can be neglected. Presumably, nonlinear fitting after transformation (probably logarithmic) to achieve a normal error distribution would perform better that either of the above alternatives (see Chapter 2). Equation (9.14) is probably an adequate and simple approach if highly precise answers are not required, and prey depletion can be neglected.

More fundamental than the technical details of curve fitting, however, is the question whether a functional response measured in the laboratory bears any resemblance to the functional response between the same pair of species in the field. There are two basic issues here, which are closely related.

First, there is the question of scale. Unless the system being parameterized is a laboratory system, or conceivably a small greenhouse system, the objective will be to apply parameter estimates obtained in small arenas to a much larger ‘natural’ environment. A dimensional analysis of any functional response can provide some information on the appropriate way to apply parameter estimates obtained on a small scale to larger arenas. For example, consider the Holling type II equation, with prey measured as a density per unit area:

D is the density of prey per unit area. f(D) has units of prey consumed per unit area per unit time. α and β thus both have units of time−1. This suggests, to a first approximation, that it should be possible to use estimates of both parameters measured in small arenas directly in larger-scale situations, provided prey are

256 C H A P T E R 9

measured in consistent units of density per unit area. Similar arguments would also apply to predators foraging in three dimensions, except that density would need to be measured per unit volume. However, the foraging behaviour of a predator or parasitoid in a 1 m2 enclosure, surrounded by barriers, is likely to be quite different from that of a predator foraging within 1 m2 of an infinite plane with the same prey density. Simple application of parameter estimates from the small scale to the large scale is likely to be unreliable, even in the improbable event that the prey distribution can be considered to be homogeneous.

Second, there is the question of habitat and environmental complexity. It is almost inevitable that artificial experimental arenas will be less spatially heterogeneous than the natural environment, and that they will have fewer species of alternative prey or competing exploiters present. It is very well established that the shape of a functional response depends on the spatial arrangement of the prey (Murdoch & Oaten, 1975; Murdoch & StewartOaten, 1989; Godfray & Pacala, 1992; Mills & Getz, 1996). Type III functional responses frequently occur as a result of predators switching from alternative prey, as the species under consideration becomes more common. Hence, the abundance and spatial arrangement of alternative prey may have a major influence on the qualitative shape of a functional response.

The net result of the above discussion is that it is most unlikely that a functional response measured in a small, artificial arena will be able to provide even qualitative insight into the dynamics of the interaction between the same predator and prey species in a natural environment. It is therefore essential that models be parameterized using field data.

Using field data

Functional responses

Measuring the functional response itself in the field requires recording the predation rate, per predator, as a function of prey density. In principle, this could be done using either a prey-centred approach or a predator-centred approach.

A prey-centred approach involves recording, for a sample of prey individuals, the death rate from predation, by the predator in question, over a range of prey densities. In addition, the population size of the predator in question must also be known or estimated. Otherwise, it is a total response, rather than a functional response, which is being estimated. Ideally, one would also want the prey density to be manipulated by the observer, rather than natural spatial or temporal variation being used. This would eliminate the possibility of confounding variables affecting both the predation rate and prey density. Furthermore, in this ideal experiment, one would want a range of predator

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 257

densities to be explored, to examine the possibility of complications due to varying predator density.

A predator-centred approach involves recording the prey consumption rate of a sample of the predators in question, for a variety of prey densities. As this method records the predation rate per predator directly, it estimates a functional and not a total response. Ideally, as in the prey-centred approach, one would prefer the prey and predator densities to be manipulated, rather than merely observed.

We do not, of course, live in an ideal world, so most real attempts to measure functional responses fall short of these ideals. The best strategy will often be determined by particular idiosyncrasies of the biology of predator and prey in the particular environment being studied, or by logistic considerations. I will therefore review briefly some case studies in which functional responses have been measured, in the hope that readers may find an approach that can be modified for their situation, rather than attempting to make general recommendations.

Functional responses of predators

Boutin (1995) reviewed attempts to measure functional responses of predators preying on small mammals with fluctuating populations. These studies looked at prey remains, either in scats, stomach contents or raptor pellets. For example, Erlinge et al. (1983) attempted to determine the functional response of a variety of mammalian and avian predators feeding on small rodents in Sweden. The information came from the proportions of various prey species identified in the scats or pellets of the predators. Together with data on the biomass of the prey, the food requirements per day of each predator, and the abundance of the predators, the number of prey individuals consumed was determined. For technical information on appropriate correction factors for diet analysis, see Simmons et al. (1991) (raptors) or Reynolds and Aebischer (1991) (red fox).

Boutin (1995) points out that such methods are likely to be unreliable, particularly at high and low prey densities. At high densities, surplus killing, caching and partial prey consumption are likely to lead to an underestimate of the actual number of prey consumed, and at low densities, if food is in short supply, fewer scats or pellets may be produced, leading to a decrease in prey consumption not fully being reflected in the reduction of the proportion of remains per scat or pellet. Furthermore, the overall food requirements of the predators usually are obtained from a diversity of sources (as was the case in the study of Erlinge et al., 1983), some of doubtful accuracy, such as laboratory feeding trials or allometric calculations (Box 9.1 summarizes some allometric relationships for calculating food consumption rate).

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 259

Box 9.1 contd

Relationships for nontetrapod consumers or poikilothermic herbivores do not seem to be as firmly established.

To use either of the above relationships in practice, it is necessary to translate the ingestion rate in watts to a rate of biomass per unit time, or number of individuals per unit of time. For animal prey, Peters (1983) suggests that 1 W = 0.012 kg fresh mass per day. For herbivores, given the very different energetic content of different sorts of vegetable matter, calorimetry would be needed.

Using further allometric relationships between predator and prey size, Peters gives the following relationships between numbers of prey killed per day Kr and predator body mass W (kg) for large prey eaters (e.g. raptors or carnivores) and small prey eaters (e.g. insectivores): for large preyeating homoiotherms, Kr = 3.00 W − 0.47; small prey eating homoiotherms, Kr = 137 W − 0.49; large prey eating poikilotherms, Kr = 0.22 W − 0.34; small prey eating poikilotherms, Kr = 10.1 W − 0.36. As the average mass of prey items can usually be estimated fairly easily, better accuracy can be obtained by using the relationship given above between biomass consumed per day and watts, together with the relationships between ingestion rate and body mass.

In some cases, it may be possible to do better than this. Sinclair et al. (1990) estimated the functional responses of raptors feeding on house mice by examining pellets and counting the number of right mandibles, or by halving the number of upper or lower incisors. As these raptors had been found to produce very close to one pellet per day, the counts of bone fragments could be translated into prey consumed per day, per predator, without needing to estimate prey mass or predator food requirements. Unfortunately, the resulting functional responses (see Fig. 9.3(a)) contain so much variation that it is hard to determine even if they are of type I, II or III.

As an alternative to examining scats or pellets, several workers have attempted to determine prey consumption rates directly from stomach contents. The problem is to convert field data on the average number of prey items, or amount of food, per stomach into a prey consumption rate per unit of time. This requires a model of the stomach evacuation rate, together with some assumptions about the feeding rate of the predators in question. For example, Pech et al. (1992) estimated the functional response of foxes preying on rabbits in Australia, using an approach based on the proportion of fox stomachs containing rabbits, the average amount of rabbit per stomach, an estimate of the stomach emptying time, and an assumption that the feeding rate was uniform over the nocturnal foraging period (Fig. 9.3(b)).

Fig. 9.3 Functional responses estimated from field data. (a) The functional response of diurnal raptors preying on mice in Australia (Sinclair et al., 1990). The response is based on the number of mouse skeletal elements per pellet (see text), and the prey density is

estimated from mark–recapture. Data points are based on sampling through time. The solid dots and solid line represent the functional response in a non-outbreak phase in the mouse population, and the open circles and dashed line represent the functional response in the outbreak phase. (b) The functional response of foxes preying on rabbits in semi-arid Australia (Pech et al., 1992). The response is based on the weight of rabbit in fox stomachs sampled by shooting, and the prey density is based on spotlight counts from vehicle transects. Data are based on samples through time: triangles are from the spring–summer period, and circles are from the autumn–winter period. The solid square is an anomalous point. (c) The functional response of wolves preying on moose in North America (Messier, 1994). The killing rate is the number of moose killed per wolf per 100 days. Data are based on samples through both time and space. The solid line is a type II functional response, with parameters shown in the body of the figure. (d) The functional response of birds (species pooled) preying on spruce budworms (Crawford & Jennings, 1989). Data are based on replication through space. The fitted curve is a quadratic, and the concave form is described as evidence for a type III functional response. The large influence of the outlying data point at a larval density of 1.1 × 106 ha−1 means that any tests of significance of the quadratic term should be interpreted with caution.

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 261

Crawford and Jennings (1989) estimated the functional response of forest birds preying on spruce budworm (Choristoneura fumiferana). They used natural spatial variation in budworm density, and used shot samples to collect songbirds in two-hour time blocks from dawn onwards. The number of individual budworms in each bird gut was determined by counting pairs of left and right mandibles. Given previous work, which indicated that 35% of budworm remains could be identified 2.5 h after feeding, the total number of budworms consumed was determined by dividing the observed count by 0.35. This correction is not entirely satisfactory, but will not influence the shape of the functional response. The results indicated that the gradient of the functional response increased with prey density (Fig. 9.3(d)), consistent with the lower portion of a type III functional response.

The problems of estimating prey consumption rates from stomach contents have been explored in some detail in the fisheries literature. For most fish (Elliott & Persson, 1978) and most mammals and birds (Warner, 1981), the rate of stomach emptying in the absence of feeding is exponential. Thus, the change in S(t), the amount of food in the stomach of an animal at time t, can be described by the following equation:

where α(t) is the rate of ingestion of food (α(t) may be a function of time), and c is the exponential rate of stomach emptying. This basic equation can often be solved to predict daily food consumption rates from observed data on stomach contents.

The simplest possible case is one in which predators forage for a period τ each day, and the rate of prey consumption during this foraging period is a constant α. This would mean that the daily food intake or ration R was simply.

If the rate of passage through the stomach is reasonably quick, as it is for most mammalian or avian predators (Warner, 1981), then the stomach can be assumed to be empty at the start of foraging. Given these assumptions, the stomach fullness at time t after the onset of foraging, of an average predator, can be found by solving eqn (9.16), with an initial condition of S(0) = 0, yielding

If the time from the start of the foraging period when each predator was sampled had been recorded, eqn (9.18) could be fitted directly to a set of stomach fullnesses and sampling times, and both α and c and hence R could be estimated.

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 263

estimate functional response parameters is its cost, both in maintaining enough radio-collared prey to obtain worthwhile results, and in sufficient personnel to follow up dead animals rapidly enough to determine the cause of death.

Some of the best vertebrate functional response data have been obtained from long-term studies of wolves (Canis lupus) preying on moose (Alces alces) in North America (Messier and Crête, 1985; Petersen and Page, 1988; Messier, 1994). Both predators and prey are large, predation is relatively infrequent (less than 1.1 moose per 100 wolf days) and wolves remain near a carcass for a considerable time (at least 8 days for adult moose; Messier and Crête, 1985). It is therefore possible to use a predator-centred approach, in which individual wolf packs are intensively tracked and all successful predation events for an individual pack are recorded. Some typical results are shown in Fig. 9.3(c).

Distinguishing between type II and type III functional responses

Theory suggests that the qualitative effects of predation depend critically on whether the functional response is of type II or type III. A type II functional response is destabilizing, as it has an inversely density-dependent effect at all prey densities, whereas a type III functional response is stabilizing over a range of prey densities, as it is directly density-dependent at low prey densities. Distinguishing between these two types of response is therefore important. Unfortunately, it is not straightforward.

One apparently attractive method that does not work is to compare the fit of a cubic regression with that of a quadratic regression. It is possible to test a null hypothesis that a response is type I versus an alternative that it is nonlinear by testing whether a term which is quadratic in prey density improves the fit of a regression of predation rate per predator versus prey density. This approach is valid because the type I functional response is indeed linear. A type II functional response may, in some cases, be adequately fitted with a quadratic regression, and cubic regressions will often be an adequate fit to curve with an inflection point (which a type III response has). However, none of the suggested functional forms suggested for a type II functional response in Table 9.1 is, in fact, a quadratic. Cubic and even higher-order terms will, in general, improve the fit of a polynomial to a convex curve that reaches an asymptote as does a type II functional response. A cubic term improving the fit over a quadratic is quite consistent with a type II functional response, and conversely, the failure of a cubic term to improve the fit over a quadratic is not evidence that the relationship is no more complex than a type II response.

When nonlinear relationships are being compared, the most appropriate way to determine whether a more complex model fits a given set of data better than a simpler one is to choose models so that they are nested, meaning that the simpler model can be obtained from the more complex one by setting

268 C H A P T E R 9

Predator numerical responses in the field

Almost all published studies of predator numerical response have attempted to estimate the strength of the response by simply regressing predator abundance against prey abundance. Observations have been based on existing spatial variation in prey abundance (e.g. Crawford & Jennings, 1989), temporal variation in prey density (Newsome et al., 1983; Rohner, 1995) or sometimes a combination of both (Messier, 1994). Such correlative approaches are not ideal, as a third factor (e.g. habitat quality) could affect both variables. The time delays that are a typical feature of predator–prey interactions are also not taken into account by a simple regression of predator abundance versus prey abundance, although lagged responses can easily be included if the variation in prey density is through time (e.g. Rohner, 1995). It is generally better to link prey abundance to predator reproduction (e.g. Shine & Madsen, 1997) or survival, rather than to predator abundance. If this is done, the mechanistic link is clearer, although not necessarily unequivocally demonstrated.

Plant–herbivore models

In their basic structure, the most commonly used plant–herbivore models do not differ fundamentally from Lotka–Volterra type predator–prey models (Crawley, 1983). However, as herbivores often (although not invariably) consume biomass from plants rather than killing them outright, the variable describing the plant population is usually biomass, rather than numbers of individuals. A detailed discussion of plant population models themselves, and the best way to parameterize them, is beyond the scope of this book. Crawley (1983; 1997) or Harper (1977) are good places to start for information in this area.

For herbivores themselves, the parameterization problems are the same as those for predators. It is necessary to quantify the functional response, which is the way in which the grazing off-take depends on the amount of forage present, and the numerical response, which is the way in which the population growth rate of the herbivore depends on the amount of forage. In some unusual cases, such as herbivores feeding on planktonic algae, or seed predation, the methods needed to parameterize these responses are identical to those for predator–prey models. Normally, however, the immobile nature of plants and the likelihood of partial consumption mean that the solutions to parameterization problems are quite different from those needed for predator– prey or host–parasitoid models.

Models for vertebrate herbivores

Caughley and Gunn (1993) provide a framework for modelling the population

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 269

dynamics of mammalian herbivores in variable environments. There are four basic components:

1 The plant response function, which is the growth rate of plant biomass as a function of standing biomass and relevant environmental variables. In deserts, the particular focus of the article, the fundamental environmental variable driving the plant response is rainfall. The plant response is analogous to the rNN(1 − N/K ) component in eqn (9.3), with the environmental variable driving variability in rN.

2 The herbivore’s functional response, which is typically of type II.

3 The herbivore’s numerical response, which is the rate of increase r as a function of forage biomass.

4 Possible density dependence in the herbivore population that is not directly resource related (the ‘intrinsic response’).

This basic structure is broadly the same as that of eqns (9.3) and (9.4), with the addition of the climatic forcing, and the expectation that the model will be run with short discrete time steps, rather than as a set of continuous differential equations.

Figure 9.4 shows the first three of these components for a red kangaroo (Macropus rufus) model parameterized using data from Kinchega National Park, New South Wales.

The plant response model is based on data from Robertson (1987), which should be consulted for a detailed description of the methods. In brief, the data are based on 0.25 m2 caged plots, in which pasture growth over periods of three months was recorded as a function of initial pasture biomass and rainfall over the period.

The functional response is based on data from Short (1987), and is based on a ‘graze-down’ experiment, in which kangaroos were penned in small enclosures. The food intake per animal per day was then recorded through time as a function of the decreasing vegetation cover, until the pens were essentially devegetated. As Short himself discusses, this is not the ideal way to obtain a functional response of a free-ranging grazer in a heterogeneous environment. Spalinger and Hobbs (1992) suggest a series of mechanistic models of mammalian herbivore functional responses that may be useful alternative ways of parameterizing herbivore functional responses.

Finally, the numerical response was based on analysis by Bayliss (1987), modified by Caughley (1987). Over a four-year period, finite rates of increase for three-month segments were converted to exponential annual rates of increase r, and related to pasture biomass. The strongest relationship was observed with no lag between vegetation and rate of increase. The data were then fitted to an Ivlev (1961) curve of the form

(c)Biomass of grasses and forbs (kg per ha)

Fig. 9.4 Components of Caughley’s interactive model for kangaroos. (a) Vegetation response model. This shows the increment, over a three-month period, in biomass of grasses and forbs at Menindee, New South Wales, Australia, as a function of rainfall over that period and initial biomass. One would expect the instantaneous growth increment to have an approximately parabolic shape, as small starting biomass would lead to a small growth increment. The three-month time period over which growth is measured is relatively long for these fast-growing plants, and so only the descending arm of the growth curve is detected. (b) Herbivore functional response. The fitted food intake in kilogram dry weight per kangaroo over three months, as a function of availability of forage. (c) Herbivore numerical response. From Caughley and Gunn (1993).

where V is the vegetation biomass, and a, b and d are constants. In principle, eqn (9.23) could be fitted to observed data by nonlinear least squares or maximum likelihood, provided the problems of non-independence of successive error terms was ignored. However, Bayliss (1987) used an ad hoc procedure in which the asymptote was fitted first by inspection. Caughley (1987) further

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 271

adjusted the estimates in order to make them more biologically ‘sensible’. Such fitting procedures would offend statistical purists, but if the objective is to produce a well-functioning mechanistic model, rather than to estimate parameters as an end in themselves, a certain amount of adhocery can be defended.

The outcome of a sample run of this model is shown in Fig. 9.5. Perhaps the most interesting result is the way in which the dynamics of the plant– herbivore interaction transform the essentially uncorrelated white noise of the rainfall into long runs of population increase and decrease in the kangaroo population. Repeated runs of a model of this type can be used to generate frequency distributions of animal numbers, in a similar way to population viability analyses (e.g. Burgman et al., 1993). However, the correlated noise introduced by the dynamics in this model is likely to produce quite different (and probably more reliable) predictions than the simple white noise incorporated into most population viability analysis models.

It should be apparent, however, that to parameterize this model for this particular species and location, albeit in a not entirely satisfactory way, has required an extensive research effort.

Approximate methods when few data are available

In most cases where it might be helpful to apply a predator–prey model, we cannot conduct extensive experiments, nor are there long runs of field data at various predator and prey densities. A variety of ad hoc approaches has therefore been applied. Quite obviously, these approaches are not as satisfactory as a controlled experiment or the analysis of extensive data. However, a model based on approximate estimates is often better than no model at all, provided its limitations are recognized. By their very nature, it is difficult to generalize about ad hoc models. I will therefore review a particular example.

The interaction between voles and their predators in northern Europe has attracted a great deal of attention over many years. Vole population dynamics, particularly at high latitudes, are characterized by pronounced oscillations. One of the main hypotheses put forward to explain these cycles is that predators drive them. The following model is the most recent in a series of relatively simple predator–prey models (Hanski et al., 1991; 1993a; Hanski & Korpimäki, 1995; Turchin & Hanski, 1997) that have been parameterized to investigate the hypothesis that the cycles are driven by predators.

Two sources of predation are included. Generalist predators have a range of possible prey species. Their dynamics are therefore decoupled from the vole population dynamics, and are not explicitly included in the model. The vole death rate from generalist predators is determined by their functional response. Generalists switch on to voles when voles become sufficiently common, and will switch to alternative prey when voles are rare. They should

274 C H A P T E R 9

Table 9.2 Parameters of the vole–mustelid model. See the text for explanation of how these estimated values were chosen. Where more than one value is given in the table, Turchin and Hanski (1997) investigated results for each of the listed values. Asterisked parameters are those used in the scaled predator–prey equations

in Chapter 2. In this model, the scaling reduces by two the number of parameters that need to be examined in any sensitivity analysis, and shows that it is the ratios G/K, etc. that determine the dynamic behaviour of the model. Unless it is possible to estimate such ratios directly, however, all the original parameters still need to be estimated to relate the model to the real world.

Environmental fluctuations were included by assuming that each parameter changed in each year i, randomly and independently of other parameters, such that the parameter p was replaced by

where N(0,σ 2) is a random number drawn from a normal distribution with mean 0 and standard deviation σ. Equation (9.28) is one of the least satisfactory parts of this model. In general, one would expect that environmental variation (such as an unusually warm winter) would cause several parameters to co-vary, and there is also no reason to expect that the additive variance in parameters should be constant. The model is deliberately highly simplified, to keep the number of parameters that need to be estimated down, as far as is possible. Nevertheless, eight parameters must be estimated. These are summarized in Table 9.2.

I discuss methods for estimating aggregated demographic parameters such as intrinsic growth rates and carrying capacities in Chapters 5 and 6. The intrinsic rate of increase of voles, r, required for this model is a yearly average: the peak rate of increase in summer would be r(1 + e), and the minimum rate of increase would be r(1 − e) (see eqn (9.26)). Turchin and Hanski (1997)

P R E D A T O R – P R E Y , H O S T – P A R A S I T O I D M O D E L S 275

used an estimated r = 6.0 y−1. This was based on a detailed study (Turchin & Ostfeld, 1997) of another vole species (Microtus pennsylvanicus) in another habitat entirely (Millbrook, NY, USA). In that study, r was estimated from the intercept of a regression of the realized rate of increase (see eqn (5.5)) versus population density. Turchin and Hanski estimated the intrinsic rate of increase of the specialist predators, s, in an even more approximate fashion. Weasels produce one generation per year, in a maximum of two litters, with a mean litter size of five. Assuming a 1:1 sex ratio, and no mortality whatsoever, the annual rate of increase would be ln(10/2) = 1.6. Turchin and Hanski adjusted this downwards to allow for density-independent mortality, yielding s = 1.25. In their study of M. pennsylvanicus in New York State, Turchin and Ostfeld used a regression-based method to estimate the seasonality parameter e at 0.53. On the basis that the boreal environments they were modelling were more extreme, Turchin and Hanski used e = 1.0 in their study.

As I discuss in Chapter 6, the carrying capacity K for voles should be the equilibrium density that would occur in the absence of predation. An additional complication here is that the structure of eqn (9.26) means that the carrying capacity in summer should be 2K if e = 1. Turchin and Hanski estimate that the peak carrying capacity in summer is 300 ha−1, and hence K = 150. The prey–predator ratio Q, which determines the carrying capacity of the specialist predators, proved very difficult to estimate adequately from the data available. Turchin and Hanski suggest that Q should be in the range of 20–40, based on an argument that a weasel needs about one vole per day to maintain itself, and that about 10 female voles will produce an average of one offspring per day. Given that the sex ratio is 1:1, 20 voles in total will then have a productivity of one vole per day.

The maximum consumption rate of specialist predators was derived from calculations of food requirements of weasels per gram of body weight per day, with an ad hoc adjustment for surplus killing and higher predation on juvenile voles. Turchin and Hanski describe the estimate of C = 600 y−1 as conservative. The half-saturation constant D was based on the observation that there is a critical minimum vole density Ncrit, below which weasels cannot breed. Estimates from the literature put this at about 10 ha−1. At this prey density, the consumption rate of voles by an individual weasel would be Ccrit, which was estimated (see above) at about 1 vole per day = 365 y−1. Using the third term in eqn (9.24),

or D ≈ 6 voles ha−1.

Turchin and Hanski estimated the parameters for the generalist predators using literature data on avian predators feeding on voles in several locations in