Forgot your password?
Document Actions
A functional response model of a predator population foraging in a patchy habitat
1.
Functional response models (e.g. Holling’s disc equation) that do not take the spatial
distributions of prey and predators into account are likely to produce biased estimates
of predation rates.
2.
To investigate the consequences of ignoring prey distribution and predator aggregation,
a general analytical model of a predator population occupying a patchy environment
with a single species of prey is developed.
3.
The model includes the density and the spatial distribution of the prey population,
the aggregative response of the predators and their mutual interference.
4.
The model provides explicit solutions to a number of scenarios that can be independently
combined: the prey has an even, random or clumped distribution, and the
predators show a convex, sigmoid, linear or no aggregative response.
5.
The model is parameterized with data from an acarine predator–prey system consisting
of
Phytoseiulus persimis
and
Tetranychus urticae
inhabiting greenhouse cucumbers.
6.
The model fits empirical data quite well and much better than if prey and predators
were assumed to be evenly distributed among patches, or if the predators were distributed
independently of the prey.
7.
The analyses show that if the predators do not show an aggregative response it will
always be an advantage to the prey to adopt a patchy distribution. On the other hand,
if the predators are capable of responding to the distribution of prey, then it will be an
advantage to the prey to be evenly distributed when its density is low and switch to a
more patchy distribution when its density increases. The effect of mutual interference is
negligible unless predator density is very high.
8.
The model shows that prey patchiness and predator aggregation in combination can
change the functional response at the population level from type II to type III, indicating
that these factors may contribute to stabilization of predator–prey dynamics.
A functional response model of a predator population
Journal of Animal
Blackwell Publishing Ltd
Ecology 2006
foraging in a patchy habitat
75, 948–958
GÖSTA NACHMAN
Department of Population Biology, Institute of Biology, University of Copenhagen, Denmark, Universitetsparken 15,
DK 2100 Ø Copenhagen
Summary
1. Functional response models (e.g. Holling’s disc equation) that do not take the spatial
distributions of prey and predators into account are likely to produce biased estimates
of predation rates.
2. To investigate the consequences of ignoring prey distribution and predator aggre-
gation, a general analytical model of a predator population occupying a patchy envi-
ronment with a single species of prey is developed.
3. The model includes the density and the spatial distribution of the prey population,
the aggregative response of the predators and their mutual interference.
4. The model provides explicit solutions to a number of scenarios that can be inde-
pendently combined: the prey has an even, random or clumped distribution, and the
predators show a convex, sigmoid, linear or no aggregative response.
5. The model is parameterized with data from an acarine predator–prey system consist-
ing of Phytoseiulus persimis and Tetranychus urticae inhabiting greenhouse cucumbers.
6. The model fits empirical data quite well and much better than if prey and predators
were assumed to be evenly distributed among patches, or if the predators were distributed
independently of the prey.
7. The analyses show that if the predators do not show an aggregative response it will
always be an advantage to the prey to adopt a patchy distribution. On the other hand,
if the predators are capable of responding to the distribution of prey, then it will be an
advantage to the prey to be evenly distributed when its density is low and switch to a
more patchy distribution when its density increases. The effect of mutual interference is
negligible unless predator density is very high.
8. The model shows that prey patchiness and predator aggregation in combination can
change the functional response at the population level from type II to type III, indicat-
ing that these factors may contribute to stabilization of predator–prey dynamics.
Key-words: aggregative response, mutual interference, spatial distribution, stability
Journal of Animal Ecology (2006) 75, 948–958
doi: 10.1111/j.1365-2656.2006.01114.x
involved. The most important component in such
Introduction
models is the density of prey that determines the func-
The capacity of predators to find, kill and consume tional response, i.e. the ability of a predator individual
prey plays a fundamental role in shaping the trophic to adjust its feeding rate to changes in prey density
interactions of food webs (Begon, Harper & Townsend (Solomon 1949). Phenomenological functional response
1996). The success of an individual predator depends models like the ones by Holling (1959) and Ivlev (1961)
on a combination of prey and predator traits that need (see Jeschke, Kopp & Tollrian 2002 for a review) predict
to be incorporated in a predation model to understand the predation rate as a function of prey density only.
fully the temporal and spatial dynamics of the species This will be true only if the prey is evenly distributed in
space, which limits the general applicability of these
models as most prey populations occur in aggregated
© 2006 The Author. Correspondence: Department of Population Biology,
patterns (Turchin & Kareiva 1989). To improve the pre-
Journal compilation Institute of Biology, University of Copenhagen, Denmark,
cision and realism of functional response models it is
© 2006 British Universitetsparken 15, DK 2100 Ø Copenhagen.
Ecological Society thus necessary to include the spatial distributions of
E-mail: gnachman@bi.ku.dk
allows for a rather general analysis of the ecological
both prey and predators, especially the ability of the
949
and evolutionary consequences of nonhomogeneous
predator to aggregate in patches with abundant prey
Functional
spatial distributions of prey and predators.
(Murdoch & Stewart-Oaten 1989; Ives 1992; Ives et al.
response in a patchy
1999). Furthermore, as many predators tend to con-
environment
centrate their searching efforts to patches where prey is
The empirical background
plentiful, mutual interference between searching pred-
ators also needs to be included as this is likely to deter- Spider mites and predatory mites occurring on all leaves
mine the limit beyond which further aggregation no from six cucumber plants Cucumis sativus L. were
longer pays-off (Beddington 1975; Hassell & Rogers counted, amounting to 92 293 spider mites and 24 801
1972; Van der Meer & Ens 1997). predatory mites distributed over 440 leaves. For each
As pointed out by Ives et al. (1999), a functional plant consisting of n leaves, the instantaneous per capita
predation rate f was computed as
response depends on the spatial scale on which it is
measured. Thus, the average per capita predation rate n
1
∑ ∑ y jt f jt =
f obs =
of predators occupying a complex environment with a
Y j =1 t
nonhomogeneous distribution of prey is likely to differ
′
n
ast x js y jt / Aj
1
∑∑∑
from the per capita predation rate of a single predator
Y
individual living in a small homogeneous area, even
1 + ∑ast hst x js + ηt ( y jt − 1) Aj
j =1 t s
′
s
though the mean prey density is the same. Ives et al. (1999)
eqn 1
coined the former ‘the population functional response’
and the latter ‘the behavioural functional response’. It where Y denotes the total number of predators on a plant
is usually the behavioural response that is studied ex- and fjt is the per capita predation rate of a predator indi-
perimentally, while it is the population response that is of vidual of stage t staying on leaf j with area Aj. Leaf j is
inhabited by xjs prey of stage s (s = eggs, juveniles, adults)
ecological interest, because it affects the population
and yjt predators of stage t (t = eggs, juveniles, adults). ast ′
dynamics of predators and prey.
Recently, Williams & Martinez (2004) suggested a and hst denote the stage-specific parameters of Holling’s
generalized version of Holling’s (1959) disc equation, (1959) disc equation, i.e. the attack rate and the handling
which incorporates an extra parameter (q) that, depending time, respectively, of stage t predators attacking prey of
stage s. Finally, ηt is a constant expressing the time a pred-
on its value, can change the functional response from
being of type II (convex) to type III (sigmoid). The ator of stage t wastes per encounter with conspecifics
model shows that even small deviations from a type II due to mutual interference (Murdoch 1973; Beddington
response in the direction of a type III response can have 1975; Hassell 1978). If stage structure was omitted from
profound effects on food-chain dynamics. However, eqn 1 (by assuming a fixed stage distribution across
leaves), the values of f were rather close to those obtained
the problem is that q has no obvious biological inter-
pretation as it is unclear how it relates to prey distribu- from eqn 1. This indicates that stage structure can be
tion and/or predator behaviour. Thus, it seems likely modelled implicitly by using lumped parameters derived
that q is not merely a predator-specific constant but a by weighting the relative contributions of the different
variable that depends on how the prey is distributed. stages to the parameter values (Nachman 2006).
In order to investigate this in more detail, I present a
functional response model that explicitly includes the
The model
density and spatial distribution of prey, and the aggre-
gative response and mutual interference of predators. Omitting stage distributions from eqn 1, and replacing
n
the observed distributions of X prey ( X = ∑ j=1 x j ) and
Thus, the contribution of each component to the over-
n
Y predators (Y = ∑ j=1 y j) with probability distributions
all functional response can be studied separately and
its potentially stabilizing effect on the predator–prey lead to the generalized functional response model
interaction can be assessed. ∞ ∞
1
∑ p(x ) ∑ p( y | x ) yf (x, y)
f=
In another study (Nachman 2006), empirical data eqn 2
¥ x=0 y=0
were used to quantify the spatial distributions of two-
spotted spider mites Tetranychus urticae Koch and preda- where p(x) denotes the probability that a patch is
inhabited by exactly x(x = 0, 1, 2, ... , ∞) prey and p(y | x)
tory mites Phytoseiulus persimilis Athias-Henriot, and
to estimate the per capita predation rate of P. persimilis the conditional probability that patches with x prey are
occupied by y( y = 0, 1, 2, ... , ∞) predators. f (x, y) is the
taking the actual distributions of the two species on
greenhouse cucumbers into consideration. The results instantaneous predation rate of a predator individual
indicate that prey patchiness may be an advantage staying in a patch occupied by x prey and y predators.
© 2006 The Author. to the prey at high prey densities but an advantage to f(x, y) depends on the size of a patch (which for simplicity
Journal compilation
the predator at low prey densities, primarily because is assumed to be the same for all patches), the number
© 2006 British
the specialist predator does not search at random but of prey (the behavioural functional response), the
Ecological Society,
spends relatively more time in the most profitable patches. number of predators (the aggregative response), and the
Journal of Animal
Applying the functional response model developed in the interaction between predators (mutual interference).
Ecology, 75,
present paper to the Tetranychus–Phytoseiulus system, The challenge is to find mathematical expressions to
948–958
substitute p(x), p( y | x) and f (x, y) in eqn 2. These
950
expressions should be simple enough to provide explicit
G. Nachman
convergent solutions to the infinite summations and at
the same time be reasonably realistic.
Three different types of prey distribution to substitute
p(x) will be considered in the following: (1) a clumped
distribution described by the negative binomial distri-
bution (NBD) with clumping parameter k; (2) a random
distribution described by the Poisson distribution with
Fig. 1. The expected number ( ¥x) of P. persimilis on a leaf
parameter ≈ (the mean number of prey per patch); and
inhabited by x T. urticae for two different combinations of
(3) an even distribution where all patches are occupied
mean prey (≈) and predator ( ¥) density. (a) Full line: ≈ = 400 T.
by the same density of prey. The clumping parameter k urticae/leaf and ¥ = 120 P. persimilis/leaf; (b) broken line:
can be considered either as a constant or as a function ≈ = 200 T. urticae/leaf and ¥ = 60 P. persimilis/leaf. Curves are
of prey density (Nachman 2006). computed by means of eqn A7 with parameter values given in
Table 1, except that c in graph (b) is constrained to 1·298 to
The predator distribution has two components: an
prevent ¥x from becoming negative.
aggregative response determined by the local density of
prey and some amount of spatial variation that cannot
be attributed to prey density (Chesson & Murdoch 1986; response levels off at high prey densities, only cases in
which µ is negative will be considered, encompassing con-
Hassell & Wilson 1997). The prey density-dependent
vex (m = 0) and sigmoid aggregative responses (m = 1).
component expresses the expected value of y at a given
density of prey, i.e. E( y | x) = ¥x. The prey density-
independent component expresses the variance of y for
-
a fixed x and is described by the conditional probability
2
function p(y | x) with variance σ y | x.
Three specific distributions of p(y | x) are considered to
describe the prey-independent component of predator
aggregation: (1) a clumped distribution characterized
2
by σ 2 | x > ¥x ; (2) a random distribution ( σ y | x = ¥x ); and
The expected number of predators inhabiting a patch y
(3) an even distribution ( σ y | x = 0). Applying the
2
with x prey is assumed to increase with x according to
the general expression NBD to represent the clumped distribution yields
2 2
σ y | x = ¥x + ¥x /κ , where κ expresses the tendency of
d ¥x m µx
= cx e eqn 3 the predators to clump independently of the prey. The
dx
random (Poisson) distribution appears as a special case
where m, c and µ are constants (Nachman 2006). of the NBD as κ → ∞.
Five main types of aggregative response (see, e.g.
Van der Meer & Ens 1997) can be produced by eqn 3
depending on its parameter values: (i) c = 0: the preda-
tors do not show any aggregative response; (ii) c > 0, The behavioural functional response of a predator
m = 0, µ = 0: the aggregative response increases linearly individual staying in a patch with x prey is assumed
with prey density; (iii) c > 0, m = 0, µ > 0: the response to be convex (type II). A model that describes such
accelerates with prey density; (iv) c > 0, m = 0, µ < 0, response is given by Ivlev (1961) as
the response increases with decelerating slope and
− ψx/ A
approaches an upper asymptote: (v) c > 0, m = 1, µ < 0: f ( x ) = f m (1 − e ) eqn 4
the response is sigmoid. Type (ii), (iv) and (v) corre-
where fm is the maximal predation rate per individual, ψ
spond to what Gascoigne & Lipcius (2004) classify as
type I, II and III aggregative response, respectively. a positive constant expressing the efficiency of the
General solutions for ¥x are found by integration of eqn predators to find and attack prey and A the patch area.
3 (Appendix 1). The higher the value of c, the more will For simplicity, A is assumed to be the same for all
patches (and equal to the mean patch area A ).
the predators tend to aggregate in patches with abun-
dant prey, whereas ¥x will decrease in patches with few
© 2006 The Author. prey. As ¥x cannot be negative, it sets an upper limit to
Journal compilation
how large c can be (see Fig. 1).
© 2006 British
To obtain explicit solutions to the aggregative response
Ecological Society,
function, it is necessary to specify whether the prey dis- Mutual interference among the y predators in a patch
Journal of Animal
tribution p(x) is clumped, random or even (Appendix 2). with x prey may tend to reduce the per capita predation
Ecology, 75,
Furthermore, as it is most realistic that the aggregative rate. This effect is included in eqn 4 as (cf. Royama 1992)
948–958
¥x in eqn 7 depends on prey distribution and on the
951 − ψx/ A
− ε ( y − 1)/ A
f ( x, y ) = f me (1 − e ) eqn 5
shape of the aggregative response (Appendix 2). If
Functional
where ε is a positive constant expressing the intensity of the prey is either clumped or randomly distributed and
response in a patchy
mutual interference (ε = 0 implies no mutual interference). the predators show a convex aggregative response,
environment
substitution of ¥x in eqn 7 by eqn A6 gives
fm c µx
∞
− ψx/ Á c
∑
I p( x )(1 − e
f= ) ¥ Q0 (µ ) + e
−
µ µ
¥ x = 0
Replacing f(x,y) in eqn 2 with eqn 5 yields
which leads to
∞ ∞
f
∑ p(x )(1 − e )∑ p( y | x ) ye
− ε ( y − 1)/ A
− ψx/ A
f= m eqn 6
ψ
ψ c
¥ ψ
x=0
y=0
f = fm I 1 − Q0 − + Q0 (µ ) Q0 − − Q0 µ −
A µ ¥ A A
According to Appendix 4, eqn 6 reduces to eqn 8a
− ( κ + 1)
∞
¥ where the functions Q0(·) and Q1(·) depend on the den-
fm
∑ p(x)(1 − e − ψx/ A
) ¥x 1 + x (1 − e− /A )
ε
eqn 6a
f=
κ
¥ sity and distribution of the prey (Appendix 5).
x=0
If the predators show a sigmoid aggregative response,
if p(y | x) follows a NBD with parameters ( ¥x,κ), and to
and the prey is either clumped or randomly distributed,
¥x is replaced by eqn A7, yielding
∞
fm − ε /A
∑ ( p(x )(1 − e
−¥x (1−e
− ψx/ A )
f= ) ¥x e ) eqn 6b
¥ x=0
fm
∞
∑
I p( x )(1 − e ψ )
− x/ A
f=
if p(y | x) follows a Poisson distribution with mean ¥x. ¥ x=0
Finally, if the aggregative response has no prey density-
c µx
− ( k + 1)
µ µ
≈
2
independent component (i.e. σ y | x = 0), all patches with ¥ + xe − (≈e ) 1 + (1 − e )
µ
k
x prey are expected to be inhabited by ¥x predators, so
− ε ( ¥ −1)/ A − ε ¥ /A
∑∞=0 p( y | x ) ye− ε( y−1)/A = ¥x e x ∑∞=0 p ( y | x ) ≈ ¥x e x and
1 x
−k
≈
y y
− eµ − 1 + (1 − eµ )
eqn 6 therefore becomes µ
k
∞
c ψ
ψ ψ
fm
∑( p(x )(1 − e
− ε ¥x / A
− ψx/ A
f= ) ¥x e = fm I 1 − Q0 − + Q1(µ )Q0 − − Q1 µ −
) eqn 6c
A µ ¥ A A
¥
x=0
ψ
c ψ
+ 2 Q0 µ − − Q0 (µ )Q0 −
The three special cases of eqn 6 can be solved numeri-
µ¥ A A
cally, because the infinite sums converge as x → ∞.
eqn 8b
However, approximate analytical solutions are derived
by making some simplifying assumptions. The terms If the predators show no aggregative response and prey
− ε /A
(1 + (¥x /κ)(1 − e−ε/A))−(κ+1) in eqn 6a, e−¥x (1 − e ) in eqn 6b, distribution is either clumped or random, ¥x will be
and e−ε¥ /A in eqn 6c will be equal to 1 when ε = 0 and less equal to ¥ for all x. Setting c = 0 in eqn 8b yields
x
than 1 when ε > 0. Hence, these terms represent the
ψ
inhibiting effect of mutual interference on the predation
f = f mI 1 − Q0 − eqn 8c
A
rate. As long as ε/A is small, all three terms will be close
to unity, but decrease with an increase in ¥x, which in
Finally, if the prey is evenly distributed, so that prey
turn increases with ¥. Therefore, as an approximation, the
density is ≈/A in all patches, eqn 7 reduces to
overall inhibiting effect of mutual interference is assumed
to depend on the overall density of predators ( ¥) rather − ψ ≈ /A
f = f mI (1 − e ) eqn 8d
than on ¥x. Equation 6 can therefore be replaced by
∞
fm
I ∑( p( x )(1 − e
− ψx/ A
Analyses
f= ) ¥x ) eqn 7
¥ x=0
where I is a factor accounting for mutual interference
(0 < I ≤ 1), calculated as Iclumped = (1 + (¥/κ)(1 − e−ε/A))−(κ+1),
− ε /A
−¥ (1 − e
and Ieven = e−ε¥/A, depending on whe-
)
I random = e The models require three state variables (≈, ¥, and A)
and seven parameters ( fm, ψ, ε, c, µ, k, and κ) to com-
ther the predators for a given value of x are clumped,
© 2006 The Author. pute f. Furthermore, as k may be density-dependent it
randomly or evenly distributed, respectively. It appears
Journal compilation
that for ε > 0, Iclumped < Irandom < Ieven and that Iclumped could be replaced by a function of ≈, which would add
© 2006 British
declines with decreasing κ. This implies that scatter two extra parameters to the model (Nachman 1984).
Ecological Society,
Table 1 gives the estimated parameter values for the T.
(prey density-independent variation) in the aggregative
Journal of Animal
urticae–P. persimilis system if k is density-independent,
response will increase mutual interference and thereby
Ecology, 75,
otherwise see Nachman 2006).
reduce the overall predation rate.
948–958
tional response has the capacity to be stabilizing. It
952 1. Parameter values of the model. Parameters allowing for density dependence in
Table
G. Nachman Nachman (2006)
k are given in should be emphasized that even if this criterion is met,
it does not necessarily mean that the predator–prey system
Parameter Parameter description Value
will be stable, because other conditions have to be ful-
filled as well (see Discussion). However, applying the above
7·787 day−1
Maximal predation rate
fm
ψ 0·887 cm2 criterion to the model can reveal whether incorpora-
Attack efficiency
ε 0·0402 cm2
Mutual interference coefficient tion of spatial heterogeneity, mutual interference and
} 2·220
c
predator aggregation change an otherwise type II func-
Aggregative response of P. persimilis
µ − 0·0834
tional response to a type III.
Aggregation index of T. urticae 0·091
k
κ Aggregative response of P. persimilis conditioned on ¥x 0·450
Results
The model-predicted predation rates ( fpred) were com-
pared with the empirical values ( fobs), which were obtained
by means of eqn 1. It should be emphasized that observed
and predicted values of f are not completely independ- The correlation coefficient (R) between the observed pre-
dation rates ( fobs) and the model-predicted rates ( fpred)
ent, as data from the same plants were used to calculate
the values of fobs and to estimate the model’s parameters. expresses the qualitative fit of a model to data. R was
However, no attempts were made afterwards to calibrate high and almost the same irrespective of whether k was
the parameters to the observed predation rates so as to assumed to be density-dependent (using the relationship
β
−α ≈
(1 + x/k)−k = e , where α and β are positive constants
improve the agreement between fobs and fpred.
(cf. eqn A3 in Nachman 2006) or density-independent
provided the prey was clumped and the predators showed
an aggregative response. The correlation decreased when
the predators were assumed to lack an aggregative response
The functional response of a predator has a stabilizing or when the prey was assumed to be evenly distributed
effect on the dynamics of the prey if the predation risk (Table 2).
Table 2 also shows how much of the variation in fobs
per prey increases with increasing prey density (Murdoch
& Oaten 1975). Convex functional response curves (type that could be explained by the respective models. The best
quantitative agreement (88·9%) between fobs and fpred
II), like Holling’s disc equation or Ivlev’s predation model,
lead to a monotonous decline in the predation risk for was achieved by eqn 6 assuming density independence
a given predator density, and therefore do not posses the of k, although k was found to vary with prey density
potential of being stabilizing, whereas a sigmoid (type III) (Nachman 2006). However, density dependence of k will
response may be stabilizing up to a certain prey density. only be of importance when prey density is either very low
Murdoch (1977) suggests using ∂f/∂x > f/x as a crite- or very high, but not within the range of observed densi-
ties (0·002–3·758 T. urticae cm−2). Equation 8 explained
rion to identify the density interval for which the func-
Table 2. Comparisons between observed ( fobs) and predicted predation rates ( fpred). The models are compared with respect to the correlation between fobs
and fpred and the percentage variation in fobs explained by means of the model
Predicted predation rates ( fpred)
Prey clumped Prey even
Greenhouse data Predator aggregation No predator aggregation
k density-
dependent k density-independent
Average leaf
P. per./leaf ( ¥) fobs
size (cm2) T. urt./leaf (≈)
Plant Eqn 6* Eqn 6† Eqn 8‡ Eqn 8§ Eqn 8¶
1 166·8 416·3 73·52 3·932 4·676 4·569 5·006 1·871 6·558
2 168·3 17·53 7·096 3·318 2·772 2·508 2·574 0·475 0·680
3 160·6 0·316 4·776 0·096 0·186 0·075 0·078 0·013 0·014
4 234·6 58·33 99·43 1·471 2·217 2·260 2·547 0·780 1·460
5 191·1 0·262 2·857 0·014 0·185 0·068 0·070 0·009 0·009
6 2006 The175·5
© Author. 659·6 187·4 3·362 4·299 3·388 3·832 1·914 6·580
Journal compilation
Correlation between fobs and fpred (R) 0·959 0·949 0·937 0·828 0·771
© 2006 British fobs explained by fpred
% variation in 84·7 88·9 79·6 1·6 0
Ecological Society,
Journal of Animal
*Predictions based on eqn 6 and parameter values in Table 1, but with a density-dependent prey clumping parameter (k). †Predictions based on eqn 6 but
Ecology, 75,
with density-independent k. ‡Predictions based on eqn 8 but with the same assumptions as in the previous column. §Predators are assumed to search
948–958
independently of prey distribution. ¶Prey is assumed to be evenly distributed.
953
Functional
response in a patchy
environment
Fig. 3. The effect of prey distribution and predator density on
Fig. 2. The predicted per capita predation rate ( fpred) for the functional response predicted by means of eqn 6 and
parameter values given in Table 1.
different densities of T. urticae and P. persimilis based on
eqn 6 and parameter values given in Table 1.
79·6%, which indicates that the approximations made
to derive this equation from eqn 6 were acceptable. The last
two columns of Table 2 demonstrate the importance of
taking the aggregative response and distribution of prey
into consideration when predicting predation rates.
Without predator aggregation (i.e. c = 0), the predicted
predation rates become consistently lower than fobs,
which shows the advantage to the predators of being
able to adopt nonrandom search. On the other hand, if
the prey is assumed to be evenly distributed the pre-
dicted predation rates will be lower than fobs at low prey
densities but higher than fobs when prey density is high.
Fig. 4. The instantaneous risk of predation calculated from
Fig. 2, assuming a mean predator density of one individual
per cm2.
exceed c.0·02 prey cm−2. Nonrandom search for a patchily
distributed prey (k = 0·091) leads to higher predation
The parameters obtained from the P. persimilis–T. urticae
system were used as the standard case to which changes rates than random search, in particular at low predator
in the model assumptions could be compared. All densities (Fig. 5).
predictions of f were obtained from eqn 6 assuming
density independence of k.
Discussion
Figure 2 shows that the functional response of P.
persimilis is predicted to increase when mean prey den-
sity (≈) increases and/or when mean predator density
(¥) decreases. The upper plateau is reached more steeply
when predator density is low than when it is high. Figure 3
shows that prey clumping increases the predation rate Traditionally, predator–prey (or parasitoid–host) theory
at low prey densities and decreases it at high prey densities has presumed that populations are homogeneously dis-
as compared with a random distribution of prey. The tributed in space (e.g. Vandermeer & Goldberg 2003).
result for an even prey distribution is almost identical However, in the real world this will hardly ever be true,
to that of a random prey distribution and is therefore not which means that the performance of the predators will
shown. Higher predator density reduces the per capita depend on how the prey is distributed and how the pred-
predation rate, in particular when the prey is aggregated, ators respond to this distribution. Ignoring these spatial
because of the intense mutual interference experienced effects may seriously bias the estimated predation rates
© 2006 The Author. by the predators when crowding in the same patches as at the population level (Ives et al. 1999) and may lead to
Journal compilation
the prey. A change in prey distribution from random to erroneous conclusions concerning the ability of the pred-
© 2006 British
aggregated also changes the functional response from ators to regulate the density of their prey (Hochberg &
Ecological Society,
type II to type III (Fig. 4), indicating that prey patch- Holt 1999) or exaggerate the risk of prey extinction due
Journal of Animal
iness per se may be able to stabilize the predator–prey to a predator-mediated Allee effect (Gascoigne & Lipcius
Ecology, 75,
dynamics as long as the mean prey density does not 2004). Thus, analyses of host–parasitoid models have
948–958
so it seems unlikely that such a weak type III response can
954
prevent outbreaks of spider mites. In an experimental
G. Nachman
set-up, Ryoo (1996) did not find evidence for a type III
response when he varied the distribution of T. urticae
eggs in an arena, but the area of his system (2500 cm2)
was much smaller than that of a cucumber plant. Besides,
it may not be possible to distinguish between a type II
response and a weak type III response when experi-
mental data are analysed statistically (Juliano 1993).
From an evolutionary point of view, a prey species may
Fig. 5. The effect of predator aggregation and density on the switch from a random distribution to a more clumped
functional response predicted by means of eqn 6 and
distribution in order to reduce the predation pressure,
parameter values given in Table 1.
but this strategy will only work as long as the predators
search at random. Once the predators have responded
shown that nonrandom distribution of parasitism (either by evolving searching behaviours that improve their for-
in time or space) can stabilize the interactions between aging success, the prey should be selected to reduce pre-
hosts and parasitoids with discrete generations (Hassell dation pressure by being evenly distributed at low density
et al. 1991). It seems likely that this also applies to predator– and being more patchily distributed at high densities. Prey
prey systems with overlapping generations, but aggregation may lead to the so-called ‘dilution effect’,
because such systems are not as straightforward to analyse where the risk of predation to individual prey declines
mathematically as host–parasitoid systems where the CV2 > with the number of surrounding conspecifics because
1 criterion can be used (Hassell et al. 1991; Taylor 1993), of relatively fewer predators per prey and/or because the
the analyses of predator–prey interactions have often been individual predator becomes satiated (Turchin & Kareiva
limited to whether the functional response is potentially 1989). In fact, many prey species tend to become more
stabilizing or not (Murdoch 1977; Murdoch & Stewart- aggregated as the mean density increases [see Taylor
Oaten 1989). A full analysis of predator–prey stability (1984) for a review]. This also applies to T. urticae (Nach-
based on Kolmogorov’s theorem would require additional man 1981). Hence, prey patchiness and predator aggre-
information about the relationship between prey density gation might be seen as the outcome of a coevolutionary
and the numerical response of the predators (May 1974). arms race (Janzen 1980; Abrams 2000; Lima 2002),
The present model shows that prey patchiness can although other factors than avoidance of predation may
markedly reduce the feeding rate of a predator individual also contribute to prey aggregation, e.g. the quality and
unless it is able to compensate by adopting a nonran- distribution of the prey’s food, females laying eggs in
dom searching behaviour. Such behaviour will lead to a clusters, and offspring tending to stay close to their
positive aggregative response where the majority of natal site (Turchin & Kareiva 1989; Begon et al. 1996).
predators will cluster in the most profitable prey patches.
If the degree of prey aggregation is high, the predators
may actually be able to achieve a higher predation rate
than they would obtain if the prey had been evenly dis-
tributed, but only as long as prey density is low. At high The present model represents progress relative to other
prey densities, a relatively large part of the predator analytical functional response models (e.g. Murdoch &
population will waste time by searching patches with Oaten 1975; Hassell et al. 1991; Ives et al. 1999; Williams
prey densities below average, whereas the remaining part & Martinez 2004) because it integrates prey patchiness,
spends time in patches with high prey density. Owing to predator aggregation and mutual interference between
satiation of the predators staying in the most profitable predators into a single analytical model, which can be
patches, the arithmetic mean of the predation rates aver- solved either numerically (eqn 6) or explicitly (eqn 8),
aged over all patches will thus be lower than if all pred- depending on the validity of the approximations. The
ators had been exposed to the mean prey density ≈. The model is general in the sense: (1) that it allows for var-
fact that prey aggregation benefits the predators at low iation in mean densities of both prey and predators; (2)
densities may also have implications for biological con- it includes the most common (type II) functional response
trol, because it slows down the growth rate of the prey and type at the patch level; (3) it allows for different types of
© 2006 The Author. helps the predators to survive during periods of prey prey distribution (even, random and clumped); (4) it
Journal compilation
scarcity (Murdoch & Briggs 1996). However, when the allows for different types of prey density-dependent aggre-
© 2006 British
model was applied to the P. persimilis–T. urticae system, gative responses (none, linear, hyperbolic, sigmoid); and (5)
Ecological Society,
the type III response was only regulatory up to a prey it allows the predators to exhibit mutual interference. Fur-
Journal of Animal
density of about 0·02 spider mites per cm2 (correspond- thermore, the model’s parameters can be estimated exper-
Ecology, 75,
ing to about 200 –400 individuals on a cucumber plant), imentally and interpreted in a biologically meaningful
948–958
way. Finally, the model can be used to make specific pre- computing time. The model presented here may fulfil
955
dictions that can be compared against empirical data as both purposes: On one hand, it can be used to analyse
Functional
demonstrated for the T. urticae–P. persimilis system. the importance of incorporating spatial heterogeneity
response in a patchy
The model also has limitations: in population models and to analyse the costs and ben-
environment
1. It focuses on only one prey and one predator species. efits of being patchily distributed. On the other hand,
2. It ignores variation among individuals within species. the model may be included as a submodel in a metap-
3. Patches are assumed to be identical except for vari- opulation model of predator–prey dynamics to predict
ation in number of prey and predators. the current predation rate (i.e. the predation rate from
time t to time t + ∆t where ∆t has to be so short that
4. The functional response, the aggregative response and
mutual interference are modelled as empirically based prey distribution can be regarded as constant) of the
(phenomenological) mathematical relationships that do predators occupying a patch (e.g. a plant), which in
not necessarily incorporate the underlying biological itself consists of patches (e.g. leaves). For instance, in
mechanisms in a realistic way (van der Meer & Ens 1997). the T. urticae–P. persimilis system at least three differ-
5. The model is computational complex and contains ent spatial scales can be identified: leaves, plants and an
many parameters that have to be estimated from sepa- entire field (or a greenhouse). Plants within a field can
rate experiments. If predation could be observed directly, be modelled in a metapopulation context (Nachman
an alternative way of estimating the parameters might 2001), but computing time would increase dramatically
be to fit the complete model to the observed predation if the within-plant dynamics should be modelled explicitly.
rates, but the large number of parameters makes it unlikely Instead, the present model may be used as a short-cut
that this will yield values that are biologically meaningful by modelling within-plant processes implicitly.
(see, e.g. Jost & Arditi 2001).
6. The number of patches and their distribution in space
Acknowledgements
are not modelled explicitly. In principle, the model pre-
sumes that any patch in the habitat is equally likely to The author thanks Professor Koos Boomsma, Institute
be found by a predator, allowing the predators to redis- of Biology, University of Copenhagen, for his many
tribute in response to a changing prey distribution. This constructive comments to improve the manuscript.
may be a reasonable assumption for a mobile predator
inhabiting a small system, such as P. persimilis and
References
T. urticae occupying a single cucumber plant, although
Abrams, P.A. (2000) The evolution of predator–prey interac-
the large scatter about the aggregative response indi-
tions: theory and evidence. Annual Review of Ecology and
cates that redistribution is far from being perfect. At a
Systematics, 31, 79 –105.
spatial scale larger than a single plant, the predators Bancroft, J.S. & Margolies, D.C. (1999) An individual-based
are likely to show even less coupling with their prey, model of an acarine tritrophic system: lima bean, Phaseolus
making the prey density-dependent component of the lunatus L., twospotted spider mite, Tetranychus urtieae
(Acari: Tetranychidae), and Phytoseiulus persimilis (Acari:
aggregative response less clear (i.e. c decreases) while
Phytoseiidae). Ecological Modelling, 123, 161–181.
the prey density-independent component becomes more
Beddington, J.R. (1975) Mutual interference between para-
pronounced (i.e. 1/κ increases). sites or predators, and its effect on searching efficiency.
7. Transit time for individuals moving from one patch Journal of Animal Ecology, 44, 331– 340.
to another is not included. If transit time is long, it will Begon, M., Harper, J.L. & Townsend, C.R. (1996) Ecology:
Individuals, Populations, Communities, 3rd edn. Blackwell
delay the redistribution of individuals and thereby
Scientific Publications, Oxford.
contribute to the scatter in the aggregative response.
Casas, J. (1990) Multidimensional host distribution and non-
Transit time may also reduce the overall foraging random parasitism: a case study and a stochastic model.
efficiency (Ryoo 1996). Ecology, 71, 1893 –1903.
Chesson, P.L. & Murdoch, W.W. (1986) Aggregation of risk:
relationships among host-parasitoid models. American
Naturalist, 127, 696 –715.
Gascoigne, J.C. & Lipcius, R.N. (2004) Allee effects driven by
Incorporating spatial effects into population models predation. Journal of Applied Ecology, 41, 801–810.
tends to render such models analytically intractable, Hassell, M.P. (1978) Arthropod Predator–prey Systems.
unless several simplifying assumptions are being made. Monographs in Population Biology, 13. Princeton Univer-
sity Press, Princeton, NJ.
Therefore, alternative approaches, such as simulation
Hassell, M.P. & Rogers, D.J. (1972) Insect parasite responses
of population-based (e.g. Nachman 2001) and individual-
in the development of population models. Journal of Animal
based (e.g. Casas 1990; Bancroft & Margolies 1999) Ecology, 41, 661– 672.
models, are usually preferred in order not to sacrifice Hassell, M.P. & Wilson, H.B. (1997) The dynamics of spatially
© 2006 The Author. realism. However, instead of regarding analytical models distributed host-parasitoid systems. Spatial Ecology (eds
Journal compilation D. Tilman & P. Kareiva), pp. 75 –110. Princeton University
as alternatives to simulation models, they may rather
© 2006 British Press, Princeton, NJ.
be considered as supplements serving two purposes: (1)
Ecological Society, Hassell, M.P., May, R.M., Pacala, S.W. & Chesson, P.L.
they are better to provide theoretical insight into the
Journal of Animal (1991) The persistence of host–parasitoid associations in
factors that affect population dynamics, and (2) they
Ecology, 75, patchy environments. I. A general criterion. American
can serve as submodels in simulation models to save Naturalist, 138, 568 – 583.
948–958
956 Hochberg, M.E. & Holt, R.D. (1999) The uniformity and Nachman, G. (1981) A mathematical model of the functional
density of pest exploitation as guides to success in biologi- relationship between density and spatial distribution of a
G. Nachman
cal control. Theoretical Approaches to Biology Control (eds population. Journal of Animal Ecology, 50, 453 –460.
B.A. Hawkins & H.V. Cornell), pp. 71–88. Cambridge Uni- Nachman, G. (1984) Estimates of mean population density
versity Press, Cambridge. and spatial distribution of Tetranychus urticae and Phyto-
Holling, C.S. (1959) Some characteristics of simple types of pre- seiulus persimilis based upon the proportion of empty
dation and parasitism. Canadian Entomologist, 91, 385 –398. sampling units. Journal of Applied Ecology, 21, 903–913.
Ives, A.R. (1992) Density-dependent and density-independent Nachman, G. (2001) Predator–prey interactions in a nonequi-
parasitoid aggregation in model host-parasitoid systems. librium context: the metapopulation approach to modelling
American Naturalist, 140, 912 – 937. ‘hide-and-seek’ dynamics in a spatially explicit tri-trophic
Ives, A.R., Schooler, S.S., Jagar, V.J., Knuteson, S.E., Grbic, system. Oikos, 94, 72 – 88.
M. & Settle, W.H. (1999) Variability and parasitoid forag- Nachman, G. (2006) The effect of prey patchiness, predator
ing efficiency: a case study of pea aphids and Aphidius ervi. aggregation, and mutual interference on the functional
American Naturalist, 154, 652 – 673. response of Phytoseiulus persimilis feeding on Tetranychus
Ivlev, V.S. (1961) Experimental Ecology of the Feeding of urticae (Acari: Phytoseiidae, Tetranychidae). Experimental
Fishes. Yale University Press, New Haven, CT. and Applied Acarology, 38, 87–111.
Janzen, D.H. (1980) When is it coevolution? Evolution, 34, Royama, T. (1992) Analytical Population Dynamics.
611– 612. Chapman & Hall, London.
Jeschke, J.J., Kopp, M. & Tollrian, R. (2002) Predator func- Ryoo, M.I. (1996) Influence of the spatial distribution pattern
tional responses: discriminating between handling and of prey among patches and spatial coincidence on the func-
digesting prey. Ecological Monographs, 72, 95 –112. tional and numerical response of Phytoseiulus persimilis
Jost, C. & Arditi, R. (2001) From pattern to process: identi- (Acarina, Phytoseiidae). Journal of Applied Entomology,
fying predator–prey models from time series data. Popula- 120, 187 –192.
tion Ecology, 43, 229 – 243. Solomon, M.E. (1949) The natural control of animal popu-
Juliano, S.A. (1993) Nonlinear curve fitting: predation and lations. Journal of Animal Ecology, 18, 1– 35.
functional response curves. Design and Analysis of Ecolog- Taylor, L.R. (1984) Assessing and interpreting the spatial dis-
ical Experiments (eds S.M. Scheiner & J. Gurevitch), pp. tributions of insect populations. Annual Review of Entomol-
159 –182. Chapman & Hall, New York. ogy, 29, 321– 357.
Lima, S.L. (2002) Putting predators back into behavioral Taylor, A.D. (1993) Heterogeneity in host–parasitoid interac-
tions: ‘Aggregation of risk’ and the ‘CV2 > 1 rule’. Trends in
predator–prey interactions. Trends in Ecology and Evolution,
17, 70 –75. Ecology and Evolution, 8, 400– 405.
May, R.M. (1974) Stability and Complexity in Model Ecosys- Turchin, P. & Kareiva, P. (1989) Aggregation in Aphis varians:
tems. Monographs in Population Biology, 2nd edn. Princ- an effective strategy for reducing predation risk. Ecology,
eton University Press, Princeton, NJ. 70, 1008 –1016.
Murdoch, W.W. (1973) The functional response of predators. Van der Meer, J. & Ens, B.J. (1997) Models of interference
Journal of Applied Ecology, 14, 335 – 341. and their consequences for the spatial distribution of ideal
Murdoch, W.W. (1977) Stabilizing effects of spatial heteroge- and free predators. Journal of Animal Ecology, 66, 846–
neity in predator–prey systems. Theoretical Population 858.
Biology, 11, 252 – 273. Vandermeer, J.H. & Goldberg, D.E. (2003) Population
Murdoch, W.W. & Briggs, C.J. (1996) Theory for biological Ecology – First Principles. Princeton University Press,
control: recent developments. Ecology, 77, 2001– 2013. Princeton, NJ.
Murdoch, W.W. & Oaten, A. (1975) Predation and population Williams, R.J. & Martinez, N.D. (2004) Stabilization of
stability. Advances in Ecological Research, 9, 1–125. chaotic and non-permanent food-web dynamics. European
Murdoch, W.W. & Stewart-Oaten, A. (1989) Aggregation by Physical Journal of B, 38, 297 –303.
parasitoids and predators: effects on equilibrium and sta-
bility. American Naturalist, 134, 288 – 310. Received 30 August 2005; accepted 28 February 2006
Appendix 1. The aggregative response
∞
Integration of eqn 3 and constraining the solution so that Σ x=0 p( x ) ¥x = ¥ and ¥x ≥ 0 for x ≥ 0 yield the following
solutions:
1 No response (i.e. c = 0):
¥x = ¥ eqn A1
2 Linear (i.e. c > 0, m = 0, µ = 0):
¥
¥ x = ¥ c( x − ≈ ) c ≤
+ eqn A2
≈
3 Accelerating (i.e. c > 0, m = 1, µ = 0):
© 2006 The Author.
c 2
∞
Journal compilation 2¥
c2
¥ x = ¥ + x − ∑ p( x )x = ¥ + [x − ( σ + ≈ )] c ≤ 2
2 2 2
eqn A3
© 2006 British 2
σ +≈
2 2
x=0
Ecological Society,
Journal of Animal
where δ2 is the sportial variance of the prey population.
Ecology, 75,
4 Convex (i.e. c > 0, m = 0, µ < 0):
948–958
957
c µx
∞
Functional µ¥
∑ p(x )e
µx
c ≤
¥ x = ¥ + e − eqn A4
response in a patchy ∞
µ
∑ p(x )e
x=0 µx
−1
environment
x=0
5 Sigmoid (i.e. c > 0, m = 1, µ < 0):
c µx 1 µx
∞ ∞ 2
µ¥
∑ p(x )xe ∑ p(x )e c ≤
µx µx
¥ x = ¥ + xe − − e − eqn A5
∞
µx
µ µ
n
1 + µ ∑ p( x )xe − ∑ p(x )e
x=0 x=0 µx
x=0
j =1
Appendix 2. The aggregative response for specific prey distributions
Replacing p(x) in eqns A4 and A5 with the individual terms of the NBD (see Appendix 3) yields
−k
c µx µ µ¥
≈
m = 0: ¥ x = ¥ e − 1 + (1 − e ) c ≤
+
−k
µ eqn A6
k µ
≈
1 + (1 − e ) − 1
k
c µx
−k
1 µx
− ( k+1)
xe − ( ≈e µ )1 + ≈ (1 − e µ ) µ
≈
− e − 1 + (1 − e )
m = 1: ¥ x = ¥ +
µ
µ
k k
eqn A7
2
µ¥
c ≤
−k
− ( k+1)
µ µ µ
≈ ≈
− 1 + (1 − e )
1 + µ≈e 1 + (1 − e )
k k
Replacing p(x) with the individual terms of the Poisson distribution (Appendix 3) yields:
µ¥
c µx − ≈ (1 − e )
µ
m = 0: ¥ x = ¥ + e −e c ≤ − ≈ (1 − e µ ) eqn A8
µ −1
e
c µx − ≈ (1 − e )
2
µ¥
1 µx
µ µ
µ − ≈ (1 − e )
m = 1: ¥ x = ¥ + xe − ≈e − (e − e ) c ≤ µ eqn A9
µ
µ µ µ − ≈ (1 − e ) − ≈ (1 − e )
1 + µ ≈e −e
As x is equal to ≈ in all patches, the expected number of predators per patch is given by eqn A1.
Appendix 3
∞ −k ∞ − ( k+1)
ax a ax a a
Derivation of ∑ x=0 p( x )e = (1 + ( ≈ /k )(1 − e )) and ∑ x=0 p( x )xe = ≈e (1 + ( ≈ /k )(1 − e ))
Defining p = k/(x + k) and q = 1 − p = x/(x + k) means that the individual terms of the NBD can be written
as p(x) = [Γ(x + k)/x!Γk]pkqx. This means that
k∞
∞ ∞
Γ( x + k ) k x ax p Γ( x + k )
∑ ∑ x!Γk p q e = p′ ∑
ax k x
= ( p′ ) (q′ )
p( x )e eqn A10
x!Γk
x=0 x=0 x=0
© 2006 The Author.
Journal compilation
where q′ = qea and p′ = 1 − q′. Provided q′ < 1, which requires that a < ln(1 + k/≈), the right-hand sum in eqn A10
© 2006 British
will converge to unity, yielding
Ecological Society,
Journal of Animal −k
k
∞
p a
≈
∑ p(x )e
ax
= = 1 + (1 − e )
Ecology, 75, eqn A11
p′
k
948–958 x=0
As (1 + 1/u)u = e for u → ∞, we may define 1/u = x/k(1 − ea) yielding −k = −≈(1 − ea)u. Hence, if p(x) is distributed
958
according to the Poisson distribution (corresponding to k → ∞) we get
G. Nachman
a
−k − ≈ (1 − e )u
∞
a 1
≈
∑ p(x )e
a
− ≈ (1 − e )
ax
= 1 + (1 − e ) = 1 + =e eqn A12
u
k
x=0
∞ ∞ ∞
ax ax k k x
To solve ∑ x=0 p( x )xe , apply eqn A10 to write ∑ x=0 p( x )xe = ( p/ p ′ ) ∑ x=0 [Γ( x + k )]/( x!Γk )( p ′ ) ( q ′ ) x , but as
∞ ∞ k x
≈ ∑ x=0 p( x )x = kq/ p, ∑ x=0 [Γ( x + k )]/( xΓk )( p ′ ) ( q ′ ) x can be replaced by kq′/p′ (provided q′ < 1), which leads to
=
− ( k+1)
∞
p k kq ′ a a
≈
∑ p(x )xe
ax
= = ≈e 1 + (1 − e ) eqn A13
p′ p′
k
x=0
For p(x) being distributed according to the Poisson distribution, eqn A13 reduces to
− ( k+1)
∞
a a
≈
∑
a
a − ≈ (1 − e )
ax
= ≈e 1 + (1 − e ) = ≈e
p( x )xe eqn A14
k
x=0
Appendix 4
− ( κ +1)
∞ − ε ( y−1)/ A − ε/ A
Derivation of ∑ x=0 p( y|x ) ye = ¥x [1 + ( ¥x /κ )(1 − e if p(y | x) is a NBD with parameters ¥x and κ.
)]
− ε/ A
In line with Appendix 3, we define p = κ /( ¥x + κ ), q = 1 − p = ¥x /( ¥x + κ ), q ′ = ( ¥x e )/( ¥x + κ ) , and p′ = 1 − q′.
∞ − ε ( y−1)/ A
∑ x=0 p( y | x ) ye can therefore be written as
κ∞
∞
p Γ( κ + y )
∑ p( y | x ) ye ∑
− ε ( y − 1)/ A ε/ A κ y
=e ( p′ ) (q′ ) y
y!Γκ
p′
y=0 y=0
κ −κ −1
p κq ′ ε/ A − ε/ A − ε/ A
¥ ¥
ε/ A − ε/ A
=e = e 1 + x (1 − e ) 1 + x (1 − e ) ¥x e eqn A15
κ κ
p′ p′
− ( κ +1)
− ε/ A
¥
= ¥x 1 + x (1 − e )
κ
If the predators are randomly distributed for a given ¥x (i.e. κ → ∞), eqn A15 reduces to
∞
− ε /A
∑ p( y | x ) ye
−¥x (1 − e
− ε ( y − 1)/ A )
= ¥x e eqn A16
y=0
Appendix 5. The Q-function
The Q-function is defined (see also Appendix 3) as
− ( k+b )
∞
a b a
≈
∑
b ax
Qb ( a ) = = ( ≈e ) 1 + (1 − e ) ( a < ln(1 + k / ≈ ); b = 0 or 1)
p( x )x e eqn A17
k
x=0
which gives
−k
a
≈
Q0 ( a ) = 1 + (1 − e ) eqn A18
k
− ( k+1)
a a
≈
Q1( a ) = ≈e 1 + (1 − e ) eqn A19
k
for b = 0 and b = 1, respectively, if p(x) follows the NBD with parameters ≈ and k, and
∞
∑ p(x )x e
a
a b − ≈ (1 − e )
b ax
Qb ( a ) = = ( ≈e ) e ( −∞ < a < ∞; b = 0 or 1) eqn A20
x=0
which becomes
© 2006 The Author.
Journal compilation a
− ≈ (1 − e )
Q0 ( a ) = e eqn A21
© 2006 British
Ecological Society, a
a − ≈ (1 − e )
Q1( a ) = ≈e eqn A22
Journal of Animal
Ecology, 75,
if p(x) follows the Poisson distribution with parameter ≈.
948–958
Journal of Animal
Blackwell Publishing Ltd
Ecology 2006
foraging in a patchy habitat
75, 948–958
GÖSTA NACHMAN
Department of Population Biology, Institute of Biology, University of Copenhagen, Denmark, Universitetsparken 15,
DK 2100 Ø Copenhagen
Summary
1. Functional response models (e.g. Holling’s disc equation) that do not take the spatial
distributions of prey and predators into account are likely to produce biased estimates
of predation rates.
2. To investigate the consequences of ignoring prey distribution and predator aggre-
gation, a general analytical model of a predator population occupying a patchy envi-
ronment with a single species of prey is developed.
3. The model includes the density and the spatial distribution of the prey population,
the aggregative response of the predators and their mutual interference.
4. The model provides explicit solutions to a number of scenarios that can be inde-
pendently combined: the prey has an even, random or clumped distribution, and the
predators show a convex, sigmoid, linear or no aggregative response.
5. The model is parameterized with data from an acarine predator–prey system consist-
ing of Phytoseiulus persimis and Tetranychus urticae inhabiting greenhouse cucumbers.
6. The model fits empirical data quite well and much better than if prey and predators
were assumed to be evenly distributed among patches, or if the predators were distributed
independently of the prey.
7. The analyses show that if the predators do not show an aggregative response it will
always be an advantage to the prey to adopt a patchy distribution. On the other hand,
if the predators are capable of responding to the distribution of prey, then it will be an
advantage to the prey to be evenly distributed when its density is low and switch to a
more patchy distribution when its density increases. The effect of mutual interference is
negligible unless predator density is very high.
8. The model shows that prey patchiness and predator aggregation in combination can
change the functional response at the population level from type II to type III, indicat-
ing that these factors may contribute to stabilization of predator–prey dynamics.
Key-words: aggregative response, mutual interference, spatial distribution, stability
Journal of Animal Ecology (2006) 75, 948–958
doi: 10.1111/j.1365-2656.2006.01114.x
involved. The most important component in such
Introduction
models is the density of prey that determines the func-
The capacity of predators to find, kill and consume tional response, i.e. the ability of a predator individual
prey plays a fundamental role in shaping the trophic to adjust its feeding rate to changes in prey density
interactions of food webs (Begon, Harper & Townsend (Solomon 1949). Phenomenological functional response
1996). The success of an individual predator depends models like the ones by Holling (1959) and Ivlev (1961)
on a combination of prey and predator traits that need (see Jeschke, Kopp & Tollrian 2002 for a review) predict
to be incorporated in a predation model to understand the predation rate as a function of prey density only.
fully the temporal and spatial dynamics of the species This will be true only if the prey is evenly distributed in
space, which limits the general applicability of these
models as most prey populations occur in aggregated
© 2006 The Author. Correspondence: Department of Population Biology,
patterns (Turchin & Kareiva 1989). To improve the pre-
Journal compilation Institute of Biology, University of Copenhagen, Denmark,
cision and realism of functional response models it is
© 2006 British Universitetsparken 15, DK 2100 Ø Copenhagen.
Ecological Society thus necessary to include the spatial distributions of
E-mail: gnachman@bi.ku.dk
allows for a rather general analysis of the ecological
both prey and predators, especially the ability of the
949
and evolutionary consequences of nonhomogeneous
predator to aggregate in patches with abundant prey
Functional
spatial distributions of prey and predators.
(Murdoch & Stewart-Oaten 1989; Ives 1992; Ives et al.
response in a patchy
1999). Furthermore, as many predators tend to con-
environment
centrate their searching efforts to patches where prey is
The empirical background
plentiful, mutual interference between searching pred-
ators also needs to be included as this is likely to deter- Spider mites and predatory mites occurring on all leaves
mine the limit beyond which further aggregation no from six cucumber plants Cucumis sativus L. were
longer pays-off (Beddington 1975; Hassell & Rogers counted, amounting to 92 293 spider mites and 24 801
1972; Van der Meer & Ens 1997). predatory mites distributed over 440 leaves. For each
As pointed out by Ives et al. (1999), a functional plant consisting of n leaves, the instantaneous per capita
predation rate f was computed as
response depends on the spatial scale on which it is
measured. Thus, the average per capita predation rate n
1
∑ ∑ y jt f jt =
f obs =
of predators occupying a complex environment with a
Y j =1 t
nonhomogeneous distribution of prey is likely to differ
′
n
ast x js y jt / Aj
1
∑∑∑
from the per capita predation rate of a single predator
Y
individual living in a small homogeneous area, even
1 + ∑ast hst x js + ηt ( y jt − 1) Aj
j =1 t s
′
s
though the mean prey density is the same. Ives et al. (1999)
eqn 1
coined the former ‘the population functional response’
and the latter ‘the behavioural functional response’. It where Y denotes the total number of predators on a plant
is usually the behavioural response that is studied ex- and fjt is the per capita predation rate of a predator indi-
perimentally, while it is the population response that is of vidual of stage t staying on leaf j with area Aj. Leaf j is
inhabited by xjs prey of stage s (s = eggs, juveniles, adults)
ecological interest, because it affects the population
and yjt predators of stage t (t = eggs, juveniles, adults). ast ′
dynamics of predators and prey.
Recently, Williams & Martinez (2004) suggested a and hst denote the stage-specific parameters of Holling’s
generalized version of Holling’s (1959) disc equation, (1959) disc equation, i.e. the attack rate and the handling
which incorporates an extra parameter (q) that, depending time, respectively, of stage t predators attacking prey of
stage s. Finally, ηt is a constant expressing the time a pred-
on its value, can change the functional response from
being of type II (convex) to type III (sigmoid). The ator of stage t wastes per encounter with conspecifics
model shows that even small deviations from a type II due to mutual interference (Murdoch 1973; Beddington
response in the direction of a type III response can have 1975; Hassell 1978). If stage structure was omitted from
profound effects on food-chain dynamics. However, eqn 1 (by assuming a fixed stage distribution across
leaves), the values of f were rather close to those obtained
the problem is that q has no obvious biological inter-
pretation as it is unclear how it relates to prey distribu- from eqn 1. This indicates that stage structure can be
tion and/or predator behaviour. Thus, it seems likely modelled implicitly by using lumped parameters derived
that q is not merely a predator-specific constant but a by weighting the relative contributions of the different
variable that depends on how the prey is distributed. stages to the parameter values (Nachman 2006).
In order to investigate this in more detail, I present a
functional response model that explicitly includes the
The model
density and spatial distribution of prey, and the aggre-
gative response and mutual interference of predators. Omitting stage distributions from eqn 1, and replacing
n
the observed distributions of X prey ( X = ∑ j=1 x j ) and
Thus, the contribution of each component to the over-
n
Y predators (Y = ∑ j=1 y j) with probability distributions
all functional response can be studied separately and
its potentially stabilizing effect on the predator–prey lead to the generalized functional response model
interaction can be assessed. ∞ ∞
1
∑ p(x ) ∑ p( y | x ) yf (x, y)
f=
In another study (Nachman 2006), empirical data eqn 2
¥ x=0 y=0
were used to quantify the spatial distributions of two-
spotted spider mites Tetranychus urticae Koch and preda- where p(x) denotes the probability that a patch is
inhabited by exactly x(x = 0, 1, 2, ... , ∞) prey and p(y | x)
tory mites Phytoseiulus persimilis Athias-Henriot, and
to estimate the per capita predation rate of P. persimilis the conditional probability that patches with x prey are
occupied by y( y = 0, 1, 2, ... , ∞) predators. f (x, y) is the
taking the actual distributions of the two species on
greenhouse cucumbers into consideration. The results instantaneous predation rate of a predator individual
indicate that prey patchiness may be an advantage staying in a patch occupied by x prey and y predators.
© 2006 The Author. to the prey at high prey densities but an advantage to f(x, y) depends on the size of a patch (which for simplicity
Journal compilation
the predator at low prey densities, primarily because is assumed to be the same for all patches), the number
© 2006 British
the specialist predator does not search at random but of prey (the behavioural functional response), the
Ecological Society,
spends relatively more time in the most profitable patches. number of predators (the aggregative response), and the
Journal of Animal
Applying the functional response model developed in the interaction between predators (mutual interference).
Ecology, 75,
present paper to the Tetranychus–Phytoseiulus system, The challenge is to find mathematical expressions to
948–958
substitute p(x), p( y | x) and f (x, y) in eqn 2. These
950
expressions should be simple enough to provide explicit
G. Nachman
convergent solutions to the infinite summations and at
the same time be reasonably realistic.
Three different types of prey distribution to substitute
p(x) will be considered in the following: (1) a clumped
distribution described by the negative binomial distri-
bution (NBD) with clumping parameter k; (2) a random
distribution described by the Poisson distribution with
Fig. 1. The expected number ( ¥x) of P. persimilis on a leaf
parameter ≈ (the mean number of prey per patch); and
inhabited by x T. urticae for two different combinations of
(3) an even distribution where all patches are occupied
mean prey (≈) and predator ( ¥) density. (a) Full line: ≈ = 400 T.
by the same density of prey. The clumping parameter k urticae/leaf and ¥ = 120 P. persimilis/leaf; (b) broken line:
can be considered either as a constant or as a function ≈ = 200 T. urticae/leaf and ¥ = 60 P. persimilis/leaf. Curves are
of prey density (Nachman 2006). computed by means of eqn A7 with parameter values given in
Table 1, except that c in graph (b) is constrained to 1·298 to
The predator distribution has two components: an
prevent ¥x from becoming negative.
aggregative response determined by the local density of
prey and some amount of spatial variation that cannot
be attributed to prey density (Chesson & Murdoch 1986; response levels off at high prey densities, only cases in
which µ is negative will be considered, encompassing con-
Hassell & Wilson 1997). The prey density-dependent
vex (m = 0) and sigmoid aggregative responses (m = 1).
component expresses the expected value of y at a given
density of prey, i.e. E( y | x) = ¥x. The prey density-
independent component expresses the variance of y for
-
a fixed x and is described by the conditional probability
2
function p(y | x) with variance σ y | x.
Three specific distributions of p(y | x) are considered to
describe the prey-independent component of predator
aggregation: (1) a clumped distribution characterized
2
by σ 2 | x > ¥x ; (2) a random distribution ( σ y | x = ¥x ); and
The expected number of predators inhabiting a patch y
(3) an even distribution ( σ y | x = 0). Applying the
2
with x prey is assumed to increase with x according to
the general expression NBD to represent the clumped distribution yields
2 2
σ y | x = ¥x + ¥x /κ , where κ expresses the tendency of
d ¥x m µx
= cx e eqn 3 the predators to clump independently of the prey. The
dx
random (Poisson) distribution appears as a special case
where m, c and µ are constants (Nachman 2006). of the NBD as κ → ∞.
Five main types of aggregative response (see, e.g.
Van der Meer & Ens 1997) can be produced by eqn 3
depending on its parameter values: (i) c = 0: the preda-
tors do not show any aggregative response; (ii) c > 0, The behavioural functional response of a predator
m = 0, µ = 0: the aggregative response increases linearly individual staying in a patch with x prey is assumed
with prey density; (iii) c > 0, m = 0, µ > 0: the response to be convex (type II). A model that describes such
accelerates with prey density; (iv) c > 0, m = 0, µ < 0, response is given by Ivlev (1961) as
the response increases with decelerating slope and
− ψx/ A
approaches an upper asymptote: (v) c > 0, m = 1, µ < 0: f ( x ) = f m (1 − e ) eqn 4
the response is sigmoid. Type (ii), (iv) and (v) corre-
where fm is the maximal predation rate per individual, ψ
spond to what Gascoigne & Lipcius (2004) classify as
type I, II and III aggregative response, respectively. a positive constant expressing the efficiency of the
General solutions for ¥x are found by integration of eqn predators to find and attack prey and A the patch area.
3 (Appendix 1). The higher the value of c, the more will For simplicity, A is assumed to be the same for all
patches (and equal to the mean patch area A ).
the predators tend to aggregate in patches with abun-
dant prey, whereas ¥x will decrease in patches with few
© 2006 The Author. prey. As ¥x cannot be negative, it sets an upper limit to
Journal compilation
how large c can be (see Fig. 1).
© 2006 British
To obtain explicit solutions to the aggregative response
Ecological Society,
function, it is necessary to specify whether the prey dis- Mutual interference among the y predators in a patch
Journal of Animal
tribution p(x) is clumped, random or even (Appendix 2). with x prey may tend to reduce the per capita predation
Ecology, 75,
Furthermore, as it is most realistic that the aggregative rate. This effect is included in eqn 4 as (cf. Royama 1992)
948–958
¥x in eqn 7 depends on prey distribution and on the
951 − ψx/ A
− ε ( y − 1)/ A
f ( x, y ) = f me (1 − e ) eqn 5
shape of the aggregative response (Appendix 2). If
Functional
where ε is a positive constant expressing the intensity of the prey is either clumped or randomly distributed and
response in a patchy
mutual interference (ε = 0 implies no mutual interference). the predators show a convex aggregative response,
environment
substitution of ¥x in eqn 7 by eqn A6 gives
fm c µx
∞
− ψx/ Á c
∑
I p( x )(1 − e
f= ) ¥ Q0 (µ ) + e
−
µ µ
¥ x = 0
Replacing f(x,y) in eqn 2 with eqn 5 yields
which leads to
∞ ∞
f
∑ p(x )(1 − e )∑ p( y | x ) ye
− ε ( y − 1)/ A
− ψx/ A
f= m eqn 6
ψ
ψ c
¥ ψ
x=0
y=0
f = fm I 1 − Q0 − + Q0 (µ ) Q0 − − Q0 µ −
A µ ¥ A A
According to Appendix 4, eqn 6 reduces to eqn 8a
− ( κ + 1)
∞
¥ where the functions Q0(·) and Q1(·) depend on the den-
fm
∑ p(x)(1 − e − ψx/ A
) ¥x 1 + x (1 − e− /A )
ε
eqn 6a
f=
κ
¥ sity and distribution of the prey (Appendix 5).
x=0
If the predators show a sigmoid aggregative response,
if p(y | x) follows a NBD with parameters ( ¥x,κ), and to
and the prey is either clumped or randomly distributed,
¥x is replaced by eqn A7, yielding
∞
fm − ε /A
∑ ( p(x )(1 − e
−¥x (1−e
− ψx/ A )
f= ) ¥x e ) eqn 6b
¥ x=0
fm
∞
∑
I p( x )(1 − e ψ )
− x/ A
f=
if p(y | x) follows a Poisson distribution with mean ¥x. ¥ x=0
Finally, if the aggregative response has no prey density-
c µx
− ( k + 1)
µ µ
≈
2
independent component (i.e. σ y | x = 0), all patches with ¥ + xe − (≈e ) 1 + (1 − e )
µ
k
x prey are expected to be inhabited by ¥x predators, so
− ε ( ¥ −1)/ A − ε ¥ /A
∑∞=0 p( y | x ) ye− ε( y−1)/A = ¥x e x ∑∞=0 p ( y | x ) ≈ ¥x e x and
1 x
−k
≈
y y
− eµ − 1 + (1 − eµ )
eqn 6 therefore becomes µ
k
∞
c ψ
ψ ψ
fm
∑( p(x )(1 − e
− ε ¥x / A
− ψx/ A
f= ) ¥x e = fm I 1 − Q0 − + Q1(µ )Q0 − − Q1 µ −
) eqn 6c
A µ ¥ A A
¥
x=0
ψ
c ψ
+ 2 Q0 µ − − Q0 (µ )Q0 −
The three special cases of eqn 6 can be solved numeri-
µ¥ A A
cally, because the infinite sums converge as x → ∞.
eqn 8b
However, approximate analytical solutions are derived
by making some simplifying assumptions. The terms If the predators show no aggregative response and prey
− ε /A
(1 + (¥x /κ)(1 − e−ε/A))−(κ+1) in eqn 6a, e−¥x (1 − e ) in eqn 6b, distribution is either clumped or random, ¥x will be
and e−ε¥ /A in eqn 6c will be equal to 1 when ε = 0 and less equal to ¥ for all x. Setting c = 0 in eqn 8b yields
x
than 1 when ε > 0. Hence, these terms represent the
ψ
inhibiting effect of mutual interference on the predation
f = f mI 1 − Q0 − eqn 8c
A
rate. As long as ε/A is small, all three terms will be close
to unity, but decrease with an increase in ¥x, which in
Finally, if the prey is evenly distributed, so that prey
turn increases with ¥. Therefore, as an approximation, the
density is ≈/A in all patches, eqn 7 reduces to
overall inhibiting effect of mutual interference is assumed
to depend on the overall density of predators ( ¥) rather − ψ ≈ /A
f = f mI (1 − e ) eqn 8d
than on ¥x. Equation 6 can therefore be replaced by
∞
fm
I ∑( p( x )(1 − e
− ψx/ A
Analyses
f= ) ¥x ) eqn 7
¥ x=0
where I is a factor accounting for mutual interference
(0 < I ≤ 1), calculated as Iclumped = (1 + (¥/κ)(1 − e−ε/A))−(κ+1),
− ε /A
−¥ (1 − e
and Ieven = e−ε¥/A, depending on whe-
)
I random = e The models require three state variables (≈, ¥, and A)
and seven parameters ( fm, ψ, ε, c, µ, k, and κ) to com-
ther the predators for a given value of x are clumped,
© 2006 The Author. pute f. Furthermore, as k may be density-dependent it
randomly or evenly distributed, respectively. It appears
Journal compilation
that for ε > 0, Iclumped < Irandom < Ieven and that Iclumped could be replaced by a function of ≈, which would add
© 2006 British
declines with decreasing κ. This implies that scatter two extra parameters to the model (Nachman 1984).
Ecological Society,
Table 1 gives the estimated parameter values for the T.
(prey density-independent variation) in the aggregative
Journal of Animal
urticae–P. persimilis system if k is density-independent,
response will increase mutual interference and thereby
Ecology, 75,
otherwise see Nachman 2006).
reduce the overall predation rate.
948–958
tional response has the capacity to be stabilizing. It
952 1. Parameter values of the model. Parameters allowing for density dependence in
Table
G. Nachman Nachman (2006)
k are given in should be emphasized that even if this criterion is met,
it does not necessarily mean that the predator–prey system
Parameter Parameter description Value
will be stable, because other conditions have to be ful-
filled as well (see Discussion). However, applying the above
7·787 day−1
Maximal predation rate
fm
ψ 0·887 cm2 criterion to the model can reveal whether incorpora-
Attack efficiency
ε 0·0402 cm2
Mutual interference coefficient tion of spatial heterogeneity, mutual interference and
} 2·220
c
predator aggregation change an otherwise type II func-
Aggregative response of P. persimilis
µ − 0·0834
tional response to a type III.
Aggregation index of T. urticae 0·091
k
κ Aggregative response of P. persimilis conditioned on ¥x 0·450
Results
The model-predicted predation rates ( fpred) were com-
pared with the empirical values ( fobs), which were obtained
by means of eqn 1. It should be emphasized that observed
and predicted values of f are not completely independ- The correlation coefficient (R) between the observed pre-
dation rates ( fobs) and the model-predicted rates ( fpred)
ent, as data from the same plants were used to calculate
the values of fobs and to estimate the model’s parameters. expresses the qualitative fit of a model to data. R was
However, no attempts were made afterwards to calibrate high and almost the same irrespective of whether k was
the parameters to the observed predation rates so as to assumed to be density-dependent (using the relationship
β
−α ≈
(1 + x/k)−k = e , where α and β are positive constants
improve the agreement between fobs and fpred.
(cf. eqn A3 in Nachman 2006) or density-independent
provided the prey was clumped and the predators showed
an aggregative response. The correlation decreased when
the predators were assumed to lack an aggregative response
The functional response of a predator has a stabilizing or when the prey was assumed to be evenly distributed
effect on the dynamics of the prey if the predation risk (Table 2).
Table 2 also shows how much of the variation in fobs
per prey increases with increasing prey density (Murdoch
& Oaten 1975). Convex functional response curves (type that could be explained by the respective models. The best
quantitative agreement (88·9%) between fobs and fpred
II), like Holling’s disc equation or Ivlev’s predation model,
lead to a monotonous decline in the predation risk for was achieved by eqn 6 assuming density independence
a given predator density, and therefore do not posses the of k, although k was found to vary with prey density
potential of being stabilizing, whereas a sigmoid (type III) (Nachman 2006). However, density dependence of k will
response may be stabilizing up to a certain prey density. only be of importance when prey density is either very low
Murdoch (1977) suggests using ∂f/∂x > f/x as a crite- or very high, but not within the range of observed densi-
ties (0·002–3·758 T. urticae cm−2). Equation 8 explained
rion to identify the density interval for which the func-
Table 2. Comparisons between observed ( fobs) and predicted predation rates ( fpred). The models are compared with respect to the correlation between fobs
and fpred and the percentage variation in fobs explained by means of the model
Predicted predation rates ( fpred)
Prey clumped Prey even
Greenhouse data Predator aggregation No predator aggregation
k density-
dependent k density-independent
Average leaf
P. per./leaf ( ¥) fobs
size (cm2) T. urt./leaf (≈)
Plant Eqn 6* Eqn 6† Eqn 8‡ Eqn 8§ Eqn 8¶
1 166·8 416·3 73·52 3·932 4·676 4·569 5·006 1·871 6·558
2 168·3 17·53 7·096 3·318 2·772 2·508 2·574 0·475 0·680
3 160·6 0·316 4·776 0·096 0·186 0·075 0·078 0·013 0·014
4 234·6 58·33 99·43 1·471 2·217 2·260 2·547 0·780 1·460
5 191·1 0·262 2·857 0·014 0·185 0·068 0·070 0·009 0·009
6 2006 The175·5
© Author. 659·6 187·4 3·362 4·299 3·388 3·832 1·914 6·580
Journal compilation
Correlation between fobs and fpred (R) 0·959 0·949 0·937 0·828 0·771
© 2006 British fobs explained by fpred
% variation in 84·7 88·9 79·6 1·6 0
Ecological Society,
Journal of Animal
*Predictions based on eqn 6 and parameter values in Table 1, but with a density-dependent prey clumping parameter (k). †Predictions based on eqn 6 but
Ecology, 75,
with density-independent k. ‡Predictions based on eqn 8 but with the same assumptions as in the previous column. §Predators are assumed to search
948–958
independently of prey distribution. ¶Prey is assumed to be evenly distributed.
953
Functional
response in a patchy
environment
Fig. 3. The effect of prey distribution and predator density on
Fig. 2. The predicted per capita predation rate ( fpred) for the functional response predicted by means of eqn 6 and
parameter values given in Table 1.
different densities of T. urticae and P. persimilis based on
eqn 6 and parameter values given in Table 1.
79·6%, which indicates that the approximations made
to derive this equation from eqn 6 were acceptable. The last
two columns of Table 2 demonstrate the importance of
taking the aggregative response and distribution of prey
into consideration when predicting predation rates.
Without predator aggregation (i.e. c = 0), the predicted
predation rates become consistently lower than fobs,
which shows the advantage to the predators of being
able to adopt nonrandom search. On the other hand, if
the prey is assumed to be evenly distributed the pre-
dicted predation rates will be lower than fobs at low prey
densities but higher than fobs when prey density is high.
Fig. 4. The instantaneous risk of predation calculated from
Fig. 2, assuming a mean predator density of one individual
per cm2.
exceed c.0·02 prey cm−2. Nonrandom search for a patchily
distributed prey (k = 0·091) leads to higher predation
The parameters obtained from the P. persimilis–T. urticae
system were used as the standard case to which changes rates than random search, in particular at low predator
in the model assumptions could be compared. All densities (Fig. 5).
predictions of f were obtained from eqn 6 assuming
density independence of k.
Discussion
Figure 2 shows that the functional response of P.
persimilis is predicted to increase when mean prey den-
sity (≈) increases and/or when mean predator density
(¥) decreases. The upper plateau is reached more steeply
when predator density is low than when it is high. Figure 3
shows that prey clumping increases the predation rate Traditionally, predator–prey (or parasitoid–host) theory
at low prey densities and decreases it at high prey densities has presumed that populations are homogeneously dis-
as compared with a random distribution of prey. The tributed in space (e.g. Vandermeer & Goldberg 2003).
result for an even prey distribution is almost identical However, in the real world this will hardly ever be true,
to that of a random prey distribution and is therefore not which means that the performance of the predators will
shown. Higher predator density reduces the per capita depend on how the prey is distributed and how the pred-
predation rate, in particular when the prey is aggregated, ators respond to this distribution. Ignoring these spatial
because of the intense mutual interference experienced effects may seriously bias the estimated predation rates
© 2006 The Author. by the predators when crowding in the same patches as at the population level (Ives et al. 1999) and may lead to
Journal compilation
the prey. A change in prey distribution from random to erroneous conclusions concerning the ability of the pred-
© 2006 British
aggregated also changes the functional response from ators to regulate the density of their prey (Hochberg &
Ecological Society,
type II to type III (Fig. 4), indicating that prey patch- Holt 1999) or exaggerate the risk of prey extinction due
Journal of Animal
iness per se may be able to stabilize the predator–prey to a predator-mediated Allee effect (Gascoigne & Lipcius
Ecology, 75,
dynamics as long as the mean prey density does not 2004). Thus, analyses of host–parasitoid models have
948–958
so it seems unlikely that such a weak type III response can
954
prevent outbreaks of spider mites. In an experimental
G. Nachman
set-up, Ryoo (1996) did not find evidence for a type III
response when he varied the distribution of T. urticae
eggs in an arena, but the area of his system (2500 cm2)
was much smaller than that of a cucumber plant. Besides,
it may not be possible to distinguish between a type II
response and a weak type III response when experi-
mental data are analysed statistically (Juliano 1993).
From an evolutionary point of view, a prey species may
Fig. 5. The effect of predator aggregation and density on the switch from a random distribution to a more clumped
functional response predicted by means of eqn 6 and
distribution in order to reduce the predation pressure,
parameter values given in Table 1.
but this strategy will only work as long as the predators
search at random. Once the predators have responded
shown that nonrandom distribution of parasitism (either by evolving searching behaviours that improve their for-
in time or space) can stabilize the interactions between aging success, the prey should be selected to reduce pre-
hosts and parasitoids with discrete generations (Hassell dation pressure by being evenly distributed at low density
et al. 1991). It seems likely that this also applies to predator– and being more patchily distributed at high densities. Prey
prey systems with overlapping generations, but aggregation may lead to the so-called ‘dilution effect’,
because such systems are not as straightforward to analyse where the risk of predation to individual prey declines
mathematically as host–parasitoid systems where the CV2 > with the number of surrounding conspecifics because
1 criterion can be used (Hassell et al. 1991; Taylor 1993), of relatively fewer predators per prey and/or because the
the analyses of predator–prey interactions have often been individual predator becomes satiated (Turchin & Kareiva
limited to whether the functional response is potentially 1989). In fact, many prey species tend to become more
stabilizing or not (Murdoch 1977; Murdoch & Stewart- aggregated as the mean density increases [see Taylor
Oaten 1989). A full analysis of predator–prey stability (1984) for a review]. This also applies to T. urticae (Nach-
based on Kolmogorov’s theorem would require additional man 1981). Hence, prey patchiness and predator aggre-
information about the relationship between prey density gation might be seen as the outcome of a coevolutionary
and the numerical response of the predators (May 1974). arms race (Janzen 1980; Abrams 2000; Lima 2002),
The present model shows that prey patchiness can although other factors than avoidance of predation may
markedly reduce the feeding rate of a predator individual also contribute to prey aggregation, e.g. the quality and
unless it is able to compensate by adopting a nonran- distribution of the prey’s food, females laying eggs in
dom searching behaviour. Such behaviour will lead to a clusters, and offspring tending to stay close to their
positive aggregative response where the majority of natal site (Turchin & Kareiva 1989; Begon et al. 1996).
predators will cluster in the most profitable prey patches.
If the degree of prey aggregation is high, the predators
may actually be able to achieve a higher predation rate
than they would obtain if the prey had been evenly dis-
tributed, but only as long as prey density is low. At high The present model represents progress relative to other
prey densities, a relatively large part of the predator analytical functional response models (e.g. Murdoch &
population will waste time by searching patches with Oaten 1975; Hassell et al. 1991; Ives et al. 1999; Williams
prey densities below average, whereas the remaining part & Martinez 2004) because it integrates prey patchiness,
spends time in patches with high prey density. Owing to predator aggregation and mutual interference between
satiation of the predators staying in the most profitable predators into a single analytical model, which can be
patches, the arithmetic mean of the predation rates aver- solved either numerically (eqn 6) or explicitly (eqn 8),
aged over all patches will thus be lower than if all pred- depending on the validity of the approximations. The
ators had been exposed to the mean prey density ≈. The model is general in the sense: (1) that it allows for var-
fact that prey aggregation benefits the predators at low iation in mean densities of both prey and predators; (2)
densities may also have implications for biological con- it includes the most common (type II) functional response
trol, because it slows down the growth rate of the prey and type at the patch level; (3) it allows for different types of
© 2006 The Author. helps the predators to survive during periods of prey prey distribution (even, random and clumped); (4) it
Journal compilation
scarcity (Murdoch & Briggs 1996). However, when the allows for different types of prey density-dependent aggre-
© 2006 British
model was applied to the P. persimilis–T. urticae system, gative responses (none, linear, hyperbolic, sigmoid); and (5)
Ecological Society,
the type III response was only regulatory up to a prey it allows the predators to exhibit mutual interference. Fur-
Journal of Animal
density of about 0·02 spider mites per cm2 (correspond- thermore, the model’s parameters can be estimated exper-
Ecology, 75,
ing to about 200 –400 individuals on a cucumber plant), imentally and interpreted in a biologically meaningful
948–958
way. Finally, the model can be used to make specific pre- computing time. The model presented here may fulfil
955
dictions that can be compared against empirical data as both purposes: On one hand, it can be used to analyse
Functional
demonstrated for the T. urticae–P. persimilis system. the importance of incorporating spatial heterogeneity
response in a patchy
The model also has limitations: in population models and to analyse the costs and ben-
environment
1. It focuses on only one prey and one predator species. efits of being patchily distributed. On the other hand,
2. It ignores variation among individuals within species. the model may be included as a submodel in a metap-
3. Patches are assumed to be identical except for vari- opulation model of predator–prey dynamics to predict
ation in number of prey and predators. the current predation rate (i.e. the predation rate from
time t to time t + ∆t where ∆t has to be so short that
4. The functional response, the aggregative response and
mutual interference are modelled as empirically based prey distribution can be regarded as constant) of the
(phenomenological) mathematical relationships that do predators occupying a patch (e.g. a plant), which in
not necessarily incorporate the underlying biological itself consists of patches (e.g. leaves). For instance, in
mechanisms in a realistic way (van der Meer & Ens 1997). the T. urticae–P. persimilis system at least three differ-
5. The model is computational complex and contains ent spatial scales can be identified: leaves, plants and an
many parameters that have to be estimated from sepa- entire field (or a greenhouse). Plants within a field can
rate experiments. If predation could be observed directly, be modelled in a metapopulation context (Nachman
an alternative way of estimating the parameters might 2001), but computing time would increase dramatically
be to fit the complete model to the observed predation if the within-plant dynamics should be modelled explicitly.
rates, but the large number of parameters makes it unlikely Instead, the present model may be used as a short-cut
that this will yield values that are biologically meaningful by modelling within-plant processes implicitly.
(see, e.g. Jost & Arditi 2001).
6. The number of patches and their distribution in space
Acknowledgements
are not modelled explicitly. In principle, the model pre-
sumes that any patch in the habitat is equally likely to The author thanks Professor Koos Boomsma, Institute
be found by a predator, allowing the predators to redis- of Biology, University of Copenhagen, for his many
tribute in response to a changing prey distribution. This constructive comments to improve the manuscript.
may be a reasonable assumption for a mobile predator
inhabiting a small system, such as P. persimilis and
References
T. urticae occupying a single cucumber plant, although
Abrams, P.A. (2000) The evolution of predator–prey interac-
the large scatter about the aggregative response indi-
tions: theory and evidence. Annual Review of Ecology and
cates that redistribution is far from being perfect. At a
Systematics, 31, 79 –105.
spatial scale larger than a single plant, the predators Bancroft, J.S. & Margolies, D.C. (1999) An individual-based
are likely to show even less coupling with their prey, model of an acarine tritrophic system: lima bean, Phaseolus
making the prey density-dependent component of the lunatus L., twospotted spider mite, Tetranychus urtieae
(Acari: Tetranychidae), and Phytoseiulus persimilis (Acari:
aggregative response less clear (i.e. c decreases) while
Phytoseiidae). Ecological Modelling, 123, 161–181.
the prey density-independent component becomes more
Beddington, J.R. (1975) Mutual interference between para-
pronounced (i.e. 1/κ increases). sites or predators, and its effect on searching efficiency.
7. Transit time for individuals moving from one patch Journal of Animal Ecology, 44, 331– 340.
to another is not included. If transit time is long, it will Begon, M., Harper, J.L. & Townsend, C.R. (1996) Ecology:
Individuals, Populations, Communities, 3rd edn. Blackwell
delay the redistribution of individuals and thereby
Scientific Publications, Oxford.
contribute to the scatter in the aggregative response.
Casas, J. (1990) Multidimensional host distribution and non-
Transit time may also reduce the overall foraging random parasitism: a case study and a stochastic model.
efficiency (Ryoo 1996). Ecology, 71, 1893 –1903.
Chesson, P.L. & Murdoch, W.W. (1986) Aggregation of risk:
relationships among host-parasitoid models. American
Naturalist, 127, 696 –715.
Gascoigne, J.C. & Lipcius, R.N. (2004) Allee effects driven by
Incorporating spatial effects into population models predation. Journal of Applied Ecology, 41, 801–810.
tends to render such models analytically intractable, Hassell, M.P. (1978) Arthropod Predator–prey Systems.
unless several simplifying assumptions are being made. Monographs in Population Biology, 13. Princeton Univer-
sity Press, Princeton, NJ.
Therefore, alternative approaches, such as simulation
Hassell, M.P. & Rogers, D.J. (1972) Insect parasite responses
of population-based (e.g. Nachman 2001) and individual-
in the development of population models. Journal of Animal
based (e.g. Casas 1990; Bancroft & Margolies 1999) Ecology, 41, 661– 672.
models, are usually preferred in order not to sacrifice Hassell, M.P. & Wilson, H.B. (1997) The dynamics of spatially
© 2006 The Author. realism. However, instead of regarding analytical models distributed host-parasitoid systems. Spatial Ecology (eds
Journal compilation D. Tilman & P. Kareiva), pp. 75 –110. Princeton University
as alternatives to simulation models, they may rather
© 2006 British Press, Princeton, NJ.
be considered as supplements serving two purposes: (1)
Ecological Society, Hassell, M.P., May, R.M., Pacala, S.W. & Chesson, P.L.
they are better to provide theoretical insight into the
Journal of Animal (1991) The persistence of host–parasitoid associations in
factors that affect population dynamics, and (2) they
Ecology, 75, patchy environments. I. A general criterion. American
can serve as submodels in simulation models to save Naturalist, 138, 568 – 583.
948–958
956 Hochberg, M.E. & Holt, R.D. (1999) The uniformity and Nachman, G. (1981) A mathematical model of the functional
density of pest exploitation as guides to success in biologi- relationship between density and spatial distribution of a
G. Nachman
cal control. Theoretical Approaches to Biology Control (eds population. Journal of Animal Ecology, 50, 453 –460.
B.A. Hawkins & H.V. Cornell), pp. 71–88. Cambridge Uni- Nachman, G. (1984) Estimates of mean population density
versity Press, Cambridge. and spatial distribution of Tetranychus urticae and Phyto-
Holling, C.S. (1959) Some characteristics of simple types of pre- seiulus persimilis based upon the proportion of empty
dation and parasitism. Canadian Entomologist, 91, 385 –398. sampling units. Journal of Applied Ecology, 21, 903–913.
Ives, A.R. (1992) Density-dependent and density-independent Nachman, G. (2001) Predator–prey interactions in a nonequi-
parasitoid aggregation in model host-parasitoid systems. librium context: the metapopulation approach to modelling
American Naturalist, 140, 912 – 937. ‘hide-and-seek’ dynamics in a spatially explicit tri-trophic
Ives, A.R., Schooler, S.S., Jagar, V.J., Knuteson, S.E., Grbic, system. Oikos, 94, 72 – 88.
M. & Settle, W.H. (1999) Variability and parasitoid forag- Nachman, G. (2006) The effect of prey patchiness, predator
ing efficiency: a case study of pea aphids and Aphidius ervi. aggregation, and mutual interference on the functional
American Naturalist, 154, 652 – 673. response of Phytoseiulus persimilis feeding on Tetranychus
Ivlev, V.S. (1961) Experimental Ecology of the Feeding of urticae (Acari: Phytoseiidae, Tetranychidae). Experimental
Fishes. Yale University Press, New Haven, CT. and Applied Acarology, 38, 87–111.
Janzen, D.H. (1980) When is it coevolution? Evolution, 34, Royama, T. (1992) Analytical Population Dynamics.
611– 612. Chapman & Hall, London.
Jeschke, J.J., Kopp, M. & Tollrian, R. (2002) Predator func- Ryoo, M.I. (1996) Influence of the spatial distribution pattern
tional responses: discriminating between handling and of prey among patches and spatial coincidence on the func-
digesting prey. Ecological Monographs, 72, 95 –112. tional and numerical response of Phytoseiulus persimilis
Jost, C. & Arditi, R. (2001) From pattern to process: identi- (Acarina, Phytoseiidae). Journal of Applied Entomology,
fying predator–prey models from time series data. Popula- 120, 187 –192.
tion Ecology, 43, 229 – 243. Solomon, M.E. (1949) The natural control of animal popu-
Juliano, S.A. (1993) Nonlinear curve fitting: predation and lations. Journal of Animal Ecology, 18, 1– 35.
functional response curves. Design and Analysis of Ecolog- Taylor, L.R. (1984) Assessing and interpreting the spatial dis-
ical Experiments (eds S.M. Scheiner & J. Gurevitch), pp. tributions of insect populations. Annual Review of Entomol-
159 –182. Chapman & Hall, New York. ogy, 29, 321– 357.
Lima, S.L. (2002) Putting predators back into behavioral Taylor, A.D. (1993) Heterogeneity in host–parasitoid interac-
tions: ‘Aggregation of risk’ and the ‘CV2 > 1 rule’. Trends in
predator–prey interactions. Trends in Ecology and Evolution,
17, 70 –75. Ecology and Evolution, 8, 400– 405.
May, R.M. (1974) Stability and Complexity in Model Ecosys- Turchin, P. & Kareiva, P. (1989) Aggregation in Aphis varians:
tems. Monographs in Population Biology, 2nd edn. Princ- an effective strategy for reducing predation risk. Ecology,
eton University Press, Princeton, NJ. 70, 1008 –1016.
Murdoch, W.W. (1973) The functional response of predators. Van der Meer, J. & Ens, B.J. (1997) Models of interference
Journal of Applied Ecology, 14, 335 – 341. and their consequences for the spatial distribution of ideal
Murdoch, W.W. (1977) Stabilizing effects of spatial heteroge- and free predators. Journal of Animal Ecology, 66, 846–
neity in predator–prey systems. Theoretical Population 858.
Biology, 11, 252 – 273. Vandermeer, J.H. & Goldberg, D.E. (2003) Population
Murdoch, W.W. & Briggs, C.J. (1996) Theory for biological Ecology – First Principles. Princeton University Press,
control: recent developments. Ecology, 77, 2001– 2013. Princeton, NJ.
Murdoch, W.W. & Oaten, A. (1975) Predation and population Williams, R.J. & Martinez, N.D. (2004) Stabilization of
stability. Advances in Ecological Research, 9, 1–125. chaotic and non-permanent food-web dynamics. European
Murdoch, W.W. & Stewart-Oaten, A. (1989) Aggregation by Physical Journal of B, 38, 297 –303.
parasitoids and predators: effects on equilibrium and sta-
bility. American Naturalist, 134, 288 – 310. Received 30 August 2005; accepted 28 February 2006
Appendix 1. The aggregative response
∞
Integration of eqn 3 and constraining the solution so that Σ x=0 p( x ) ¥x = ¥ and ¥x ≥ 0 for x ≥ 0 yield the following
solutions:
1 No response (i.e. c = 0):
¥x = ¥ eqn A1
2 Linear (i.e. c > 0, m = 0, µ = 0):
¥
¥ x = ¥ c( x − ≈ ) c ≤
+ eqn A2
≈
3 Accelerating (i.e. c > 0, m = 1, µ = 0):
© 2006 The Author.
c 2
∞
Journal compilation 2¥
c2
¥ x = ¥ + x − ∑ p( x )x = ¥ + [x − ( σ + ≈ )] c ≤ 2
2 2 2
eqn A3
© 2006 British 2
σ +≈
2 2
x=0
Ecological Society,
Journal of Animal
where δ2 is the sportial variance of the prey population.
Ecology, 75,
4 Convex (i.e. c > 0, m = 0, µ < 0):
948–958
957
c µx
∞
Functional µ¥
∑ p(x )e
µx
c ≤
¥ x = ¥ + e − eqn A4
response in a patchy ∞
µ
∑ p(x )e
x=0 µx
−1
environment
x=0
5 Sigmoid (i.e. c > 0, m = 1, µ < 0):
c µx 1 µx
∞ ∞ 2
µ¥
∑ p(x )xe ∑ p(x )e c ≤
µx µx
¥ x = ¥ + xe − − e − eqn A5
∞
µx
µ µ
n
1 + µ ∑ p( x )xe − ∑ p(x )e
x=0 x=0 µx
x=0
j =1
Appendix 2. The aggregative response for specific prey distributions
Replacing p(x) in eqns A4 and A5 with the individual terms of the NBD (see Appendix 3) yields
−k
c µx µ µ¥
≈
m = 0: ¥ x = ¥ e − 1 + (1 − e ) c ≤
+
−k
µ eqn A6
k µ
≈
1 + (1 − e ) − 1
k
c µx
−k
1 µx
− ( k+1)
xe − ( ≈e µ )1 + ≈ (1 − e µ ) µ
≈
− e − 1 + (1 − e )
m = 1: ¥ x = ¥ +
µ
µ
k k
eqn A7
2
µ¥
c ≤
−k
− ( k+1)
µ µ µ
≈ ≈
− 1 + (1 − e )
1 + µ≈e 1 + (1 − e )
k k
Replacing p(x) with the individual terms of the Poisson distribution (Appendix 3) yields:
µ¥
c µx − ≈ (1 − e )
µ
m = 0: ¥ x = ¥ + e −e c ≤ − ≈ (1 − e µ ) eqn A8
µ −1
e
c µx − ≈ (1 − e )
2
µ¥
1 µx
µ µ
µ − ≈ (1 − e )
m = 1: ¥ x = ¥ + xe − ≈e − (e − e ) c ≤ µ eqn A9
µ
µ µ µ − ≈ (1 − e ) − ≈ (1 − e )
1 + µ ≈e −e
As x is equal to ≈ in all patches, the expected number of predators per patch is given by eqn A1.
Appendix 3
∞ −k ∞ − ( k+1)
ax a ax a a
Derivation of ∑ x=0 p( x )e = (1 + ( ≈ /k )(1 − e )) and ∑ x=0 p( x )xe = ≈e (1 + ( ≈ /k )(1 − e ))
Defining p = k/(x + k) and q = 1 − p = x/(x + k) means that the individual terms of the NBD can be written
as p(x) = [Γ(x + k)/x!Γk]pkqx. This means that
k∞
∞ ∞
Γ( x + k ) k x ax p Γ( x + k )
∑ ∑ x!Γk p q e = p′ ∑
ax k x
= ( p′ ) (q′ )
p( x )e eqn A10
x!Γk
x=0 x=0 x=0
© 2006 The Author.
Journal compilation
where q′ = qea and p′ = 1 − q′. Provided q′ < 1, which requires that a < ln(1 + k/≈), the right-hand sum in eqn A10
© 2006 British
will converge to unity, yielding
Ecological Society,
Journal of Animal −k
k
∞
p a
≈
∑ p(x )e
ax
= = 1 + (1 − e )
Ecology, 75, eqn A11
p′
k
948–958 x=0
As (1 + 1/u)u = e for u → ∞, we may define 1/u = x/k(1 − ea) yielding −k = −≈(1 − ea)u. Hence, if p(x) is distributed
958
according to the Poisson distribution (corresponding to k → ∞) we get
G. Nachman
a
−k − ≈ (1 − e )u
∞
a 1
≈
∑ p(x )e
a
− ≈ (1 − e )
ax
= 1 + (1 − e ) = 1 + =e eqn A12
u
k
x=0
∞ ∞ ∞
ax ax k k x
To solve ∑ x=0 p( x )xe , apply eqn A10 to write ∑ x=0 p( x )xe = ( p/ p ′ ) ∑ x=0 [Γ( x + k )]/( x!Γk )( p ′ ) ( q ′ ) x , but as
∞ ∞ k x
≈ ∑ x=0 p( x )x = kq/ p, ∑ x=0 [Γ( x + k )]/( xΓk )( p ′ ) ( q ′ ) x can be replaced by kq′/p′ (provided q′ < 1), which leads to
=
− ( k+1)
∞
p k kq ′ a a
≈
∑ p(x )xe
ax
= = ≈e 1 + (1 − e ) eqn A13
p′ p′
k
x=0
For p(x) being distributed according to the Poisson distribution, eqn A13 reduces to
− ( k+1)
∞
a a
≈
∑
a
a − ≈ (1 − e )
ax
= ≈e 1 + (1 − e ) = ≈e
p( x )xe eqn A14
k
x=0
Appendix 4
− ( κ +1)
∞ − ε ( y−1)/ A − ε/ A
Derivation of ∑ x=0 p( y|x ) ye = ¥x [1 + ( ¥x /κ )(1 − e if p(y | x) is a NBD with parameters ¥x and κ.
)]
− ε/ A
In line with Appendix 3, we define p = κ /( ¥x + κ ), q = 1 − p = ¥x /( ¥x + κ ), q ′ = ( ¥x e )/( ¥x + κ ) , and p′ = 1 − q′.
∞ − ε ( y−1)/ A
∑ x=0 p( y | x ) ye can therefore be written as
κ∞
∞
p Γ( κ + y )
∑ p( y | x ) ye ∑
− ε ( y − 1)/ A ε/ A κ y
=e ( p′ ) (q′ ) y
y!Γκ
p′
y=0 y=0
κ −κ −1
p κq ′ ε/ A − ε/ A − ε/ A
¥ ¥
ε/ A − ε/ A
=e = e 1 + x (1 − e ) 1 + x (1 − e ) ¥x e eqn A15
κ κ
p′ p′
− ( κ +1)
− ε/ A
¥
= ¥x 1 + x (1 − e )
κ
If the predators are randomly distributed for a given ¥x (i.e. κ → ∞), eqn A15 reduces to
∞
− ε /A
∑ p( y | x ) ye
−¥x (1 − e
− ε ( y − 1)/ A )
= ¥x e eqn A16
y=0
Appendix 5. The Q-function
The Q-function is defined (see also Appendix 3) as
− ( k+b )
∞
a b a
≈
∑
b ax
Qb ( a ) = = ( ≈e ) 1 + (1 − e ) ( a < ln(1 + k / ≈ ); b = 0 or 1)
p( x )x e eqn A17
k
x=0
which gives
−k
a
≈
Q0 ( a ) = 1 + (1 − e ) eqn A18
k
− ( k+1)
a a
≈
Q1( a ) = ≈e 1 + (1 − e ) eqn A19
k
for b = 0 and b = 1, respectively, if p(x) follows the NBD with parameters ≈ and k, and
∞
∑ p(x )x e
a
a b − ≈ (1 − e )
b ax
Qb ( a ) = = ( ≈e ) e ( −∞ < a < ∞; b = 0 or 1) eqn A20
x=0
which becomes
© 2006 The Author.
Journal compilation a
− ≈ (1 − e )
Q0 ( a ) = e eqn A21
© 2006 British
Ecological Society, a
a − ≈ (1 − e )
Q1( a ) = ≈e eqn A22
Journal of Animal
Ecology, 75,
if p(x) follows the Poisson distribution with parameter ≈.
948–958