Does mutual interference always stabilize predator–prey dynamics? A comparison of models

Roger Arditi; Jean-Marc Callois; Yuri Tyutyunov; Christian Jost

doi:10.1016/j.crvi.2004.06.007

Biological modelling / Biomodélisation

Does mutual interference always stabilize predator–prey dynamics? A comparison of models

Presented by: Pierre Auger

Roger Arditi ¹ ; Jean-Marc Callois ¹; Yuri Tyutyunov ¹; Christian Jost ¹

¹ « Écologie des populations et communautés », Institut national agronomique Paris-Grignon, 16, rue Claude-Bernard, 75231 Paris cedex 05, France

Comptes Rendus. Biologies, Volume 327 (2004) no. 11, pp. 1037-1057.

Abstracts

English
French

Based on a qualitative analysis of ODE systems, the dynamic properties of alternative predator–prey models with predator-dependent functional response have been compared in order to study the role that predator interference plays in the stabilisation of trophic systems. The models considered for interference have different mathematical expressions and different conceptual foundations. Despite these differences, they give essentially the same qualitative results: when interference is low, increasing it has a positive effect on asymptotic stability and thus on the resilience of the biological system. When it is high, it is the contrary (with logistic prey growth, increasing the interference parameter ensures stability but leads to very small predator densities). Possible consequences on the evolution of the interference level in real ecosystems are discussed.

À l'aide d'une analyse qualitative de systèmes d'EDO, ont été comparées les propriétés dynamiques de modèles proie-prédateur différents, incluant une réponse fonctionnelle prédateur-dépendante, dans le but d'étudier le rôle que joue l'interférence mutuelle dans la stabilisation des systèmes trophiques. Les modèles d'interférence considérés ont des expressions mathématiques et des fondements conceptuels différents. Malgré ces divergences, ils donnent pour l'essentiel les mêmes résultats qualitatifs : lorsque le degré d'interférence est bas, son augmentation exerce un effet favorable sur la stabilité asymptotique et donc sur la résilience du système biologique. Lorsqu'il est élevé, c'est le contraire (avec une croissance logistique des proies, augmenter le paramètre d'interférence maintient la stabilité, mais conduit à de très faibles densités des prédateurs). Sont ensuite discutées les conséquences possibles sur l'évolution du degré d'interférence dans les écosystèmes réels.

Metadata

Received: 2004-03-19
Accepted: 2004-06-03
Published online: 2004-11-01

PMID

DOI: 10.1016/j.crvi.2004.06.007

Keywords: trophic interactions, population dynamics, functional response, predator-dependence, ratio-dependence, enrichment
Mot clés : interactions trophiques, dynamique des populations, réponse fonctionnelle, prédateur-dépendance, ratio-dépendance, enrichissement

Author's affiliations:

Roger Arditi ¹; Jean-Marc Callois ¹; Yuri Tyutyunov ¹; Christian Jost ¹

¹ « Écologie des populations et communautés », Institut national agronomique Paris-Grignon, 16, rue Claude-Bernard, 75231 Paris cedex 05, France

@article{CRBIOL_2004__327_11_1037_0,
     author = {Roger Arditi and Jean-Marc Callois and Yuri Tyutyunov and Christian Jost},
     title = {Does mutual interference always stabilize predator{\textendash}prey dynamics? {A} comparison of models},
     journal = {Comptes Rendus. Biologies},
     pages = {1037--1057},
     publisher = {Elsevier},
     volume = {327},
     number = {11},
     year = {2004},
     doi = {10.1016/j.crvi.2004.06.007},
     language = {en},
}

TY  - JOUR
AU  - Roger Arditi
AU  - Jean-Marc Callois
AU  - Yuri Tyutyunov
AU  - Christian Jost
TI  - Does mutual interference always stabilize predator–prey dynamics? A comparison of models
JO  - Comptes Rendus. Biologies
PY  - 2004
SP  - 1037
EP  - 1057
VL  - 327
IS  - 11
PB  - Elsevier
DO  - 10.1016/j.crvi.2004.06.007
LA  - en
ID  - CRBIOL_2004__327_11_1037_0
ER  -

%0 Journal Article
%A Roger Arditi
%A Jean-Marc Callois
%A Yuri Tyutyunov
%A Christian Jost
%T Does mutual interference always stabilize predator–prey dynamics? A comparison of models
%J Comptes Rendus. Biologies
%D 2004
%P 1037-1057
%V 327
%N 11
%I Elsevier
%R 10.1016/j.crvi.2004.06.007
%G en
%F CRBIOL_2004__327_11_1037_0

Roger Arditi; Jean-Marc Callois; Yuri Tyutyunov; Christian Jost. Does mutual interference always stabilize predator–prey dynamics? A comparison of models. Comptes Rendus. Biologies, Volume 327 (2004) no. 11, pp. 1037-1057. doi : 10.1016/j.crvi.2004.06.007. https://comptes-rendus.academie-sciences.fr/biologies/articles/10.1016/j.crvi.2004.06.007/

Version originale du texte intégral

1 Introduction

In the founding models of predator–prey theory in the 1920s (the Lotka–Volterra model in continuous time or its Nicholson–Bailey equivalent in discrete time), it was assumed that the instantaneous rate of encounter of prey by a single predator individual was independent of the predator density. More formally, this rate, known as the ‘searching efficiency’ a is defined as the proportion of prey encountered per predator per unit of searching time. Although this assumption was occasionally questioned (e.g., by Volterra himself [1]), it was only much later that Hassell and Varley [2] provided empirical evidence for an adverse influence of predator density on the searching efficiency. The evidence rested on data on insect predators and parasitoids.

Hassell and Varley suggested first that this effect was due to direct behavioural interference between searching individuals but it was later shown that other mechanisms (e.g., predator aggregation) could lead to a similar effect (e.g., [3]). For this reason, we use in this paper the word ‘interference’ in a generic sense, designating any mechanism by which predator density depresses individual predatory performance. In this context, ‘density-dependence’ can be used as a synonym to ‘interference’.

In their paper, Hassell and Varley [2] quantified the intensity of interference by a parameter m, the negative slope of the log–log regression of searching efficiency against predator density. However, their estimates of searching efficiency ignored the saturation effect that superabundant prey can have on predators [4]. Arditi and Akçakaya [5] proved that this leads to underestimation of the interference constant m and that, using new estimates, interference was even more frequent and more intense. The descriptive approach initiated by Hassell and Varley has been very useful for empiricists, allowing us to assess interference with the single parameter m.

Using a mechanistic approach to describe predator behaviour, Beddington [6] developed an explicit mathematical model to describe the effect of individual interference on the consumption rate. Extending the assumptions of Holling's ‘disc model’ [4], the decline of the predator efficiency is due to the time wasted by predator encounters. Quite independently, DeAngelis et al. [7] developed almost the same model with a perspective in predator–prey population dynamics. It is worth mentioning that DeAngelis et al. followed a phenomenological approach on the scale of populations, contrary to Beddington, who reasoned in terms of individual behaviour.

A number of studies have investigated the effect of mutual interference on population dynamics and population stability. Hassell and Rogers [8] and Hassell and May [9] studied the effect of the Hassell–Varley model on a discrete-time model of the Nicholson–Bailey type. Free et al. [3] performed a similar analysis with the Beddington model. Rogers and Hassell [10] proposed another interaction model of the same kind as the Beddington model. DeAngelis et al. [7] studied the dynamic properties of a continuous-time autonomous model incorporating their interference model. This model was taken up recently by Hwang [11,12] to establish that locally stable equilibria are also globally stable, and that periodic orbits, if existing, are unique. Further mathematical results on this model are provided by Fan and Kuang [13] who studied the non-autonomous case and by Cantrell and Cosner [14,15] who studied the persistence or extinction of species in spatially explicit reaction–diffusion models. From these and other theoretical studies, and from some empirical evidence, a consensus has emerged to consider that interference has a stabilizing influence on population dynamics (e.g., [16, p. 383]), although Hassell and May [9] pointed out that there was an upper limit on the interference constant beyond which the dynamics became unstable.

In the present paper, we intend to perform a unified and comparative analysis of the effect of the major interference models on population dynamics. In particular, we will fully analyse the dynamics of a model that incorporates the Hassell–Varley interference model together with Holling's saturation in a continuous time setting.

We shall focus on the study of model equilibria and their asymptotic stability. This is motivated by the usual assumption that mathematical stability of a trivial equilibrium indicates population extinction, while a nontrivial equilibrium is indicative of the averages of real population abundances. Weak, aperiodic fluctuations or fluctuations correlated with environmental changes suggest the existence of a stable nontrivial equilibrium. Thus, the existence of a stable nontrivial equilibrium is required, in practice, for the persistence of the populations. However, even a stable equilibrium may sometimes lead to strong short-term fluctuations [17]. Periodic fluctuations with no environmental correlates are suggestive of limit cycles around an unstable equilibrium. But limit cycles may in practice lead to extinction due to demographic and environmental stochasticity, even though some of the most popular records of predator–prey interactions concern periodic oscillations of vertebrate populations that are not due to seasonal changes. A large limit cycle that periodically brings either population close to zero implies a high probability of its extinction, because of stochastic processes or simply because of external impacts that are not taken into account by the model.

We shall examine the effect of two main parameters: first, the intensity of predator interference, as it is generally suggested that it should be a stabilizing factor [6,7,18] second, the prey productivity, as it is generally suggested that it should be a destabilizing factor [19].

2 The models and their analysis

We consider deterministic aggregated predator–prey models written in the following ‘canonical’ form:

d N / d t = F (N) - g (N, P) P d P / d t = e g (N, P) P - q P

(1)

where N and P are the densities (or biomasses) of prey and of predators, respectively. The production of prey in the absence of predators is described by the function

F (N)

, whereas

g (N, P)

is the functional response (number of preys eaten per predator per unit time). Natural mortality of prey is considered to be negligible compared to mortality due to predation. The constant e is the trophic efficiency and predators are assumed to die with a constant death rate q. The function F will be taken either as the Malthusian growth

F (N) = r N

or as the logistic model

F (N) = r N (1 - N / K)

.

Many questions in predator–prey theory, including the question of interference between predators, revolve around the expression that is used for the functional response $g (N, P)$ . The well-known models of Lotka–Volterra

(LV) g (N, P) = a N

and Holling [4]

(H) g (N, P) = \frac{a N}{1 + a h N}

do not include interference since they do not actually depend on predator density. They are of the form

g = g (N)

(termed ‘prey-dependent’ by Arditi and Ginzburg [20]).

Using empirical estimates of the searching efficiency a, Hassell and Varley [2] observed that this quantity declined with predator density almost linearly on log–log scales. Therefore, they suggested the model:

(HV) a = α P^{- m}

(2)

This expression is dimensionally unbalanced [21] and it should rather be written so that m be the exponent of a dimensionless quantity like, for example,

α {(P / P_{0})}^{- m}

, where

P_{0}

is the density of a ‘population’ consisting of a single predator in the whole habitat. There is no loss of generality to choose the units for measuring habitat surface so that

P_{0} = 1

and the original expression can be used safely.

Hassell and Varley ignored saturation of the functional response and used a Lotka–Volterra-type response that was modified by replacing the attack rate a by Eq. (2):

(HVLV) g (N, P) = α N P^{- m}

If, however, saturation occurs (e.g., in the form of a handling time h), the Hassell–Varley model must be introduced into Holling's type II model. This is the model that was actually considered in the work of Arditi and Akçakaya [5]:

(HVH) g (N, P) = \frac{α N P^{- m}}{1 + α h N P^{- m}}

(3)

This work showed that this correction leads, in general, to higher estimates of the interference constant m. In most cases considered in [5], the value of m was found to be between 0.5 and 1 and, quite often, not significantly different from 1. For the special case when

m = 1

, the functional response becomes ratio-dependent [20]:

(RD) g (N, P) = \frac{α N / P}{1 + α h N / P}

(4)

In this case, interference has a perfectly compensatory effect, while it becomes overcompensating if

m > 1

. If

m < 0

, the model has generally been interpreted as describing cooperation instead of interference. This case will not be examined here.

In contrast with Hassell and Varley's empirical model, Beddington [6] proposed a simple mechanism for interference. Extending Holling's assumption that predators spent some time h on encounters with prey, he assumed that they also wasted a time w on encounters with other predators. This yields the expression:

g (N, P) = \frac{a N}{1 + a h N + a w (P - P_{0})}

(5)

where

P_{0}

is, as above, the ‘population density’ of a single predator. This expression accounts for the fact that a given predator cannot encounter itself. Note that for

w = {(a P_{0})}^{−1}

, this function becomes ratio-dependent.

Independently, DeAngelis et al. [7] proposed the same model with purely phenomenological arguments, and omitting the constant $P_{0}$ . Under reasonable conditions, this omission should be acceptable and we will consider here only the simpler form:

(BDA) g (N, P) = \frac{a N}{1 + a h N + a w P}

(6)

DeAngelis et al. [7] performed a stability analysis of the continuous-time dynamic model (1) with the functional response (6) and logistic growth. In this paper, our primary aim will be to undertake a similar analysis with the functional responses based on the Hassell–Varley model. It will be interesting to compare the dynamic properties generated by the two approaches, as there are no clear direct empirical arguments favouring one or the other.

Note that the BDA model tends to the ratio-dependent model if the interference parameter w is large, provided that predator density is not negligible: the constant 1 in the denominator of Eq. (6) can then be neglected. Thus, the ratio-dependent model can be considered as the limit of either the HVH model or the BDA model when the interference effect becomes very strong. Note that the value $m = 1$ in the HVH model corresponds to a large interference, while its BDA equivalent is a very large value of w. Since the HVLV and RD models are special cases of the HVH model, the analysis will focus on the HVH and BDA models. The results for those two special cases will also be discussed.

In the special cases of zero interference parameters (i.e., when $w = 0$ or $m = 0$ ), both BDA and HVH models become the prey-dependent Holling model (H), or even the LV model if $h = 0$ .

2.1 The isoclines

Continuous-time models, as those studied here, have the advantage over discrete-time models that they can be studied with the isocline method, thus providing a clear graphic perception of the dynamic properties. The zero isoclines define domains in the phase space $(N, P)$ in which the population rate of change does not change sign. They give a straightforward way of finding the equilibrium (intersection of the isoclines). They also give indications on the stability of the equilibria and on the global behaviour of the dynamics.

In this section, we will establish some general properties about the shape of the predator isocline and its consequences on the dynamics. Additionally, a sufficient condition for stability will be stated. The prey isocline depends on the production function $F (N)$ and will be studied later. The expression of the predator isocline holds for all types of prey growth. In particular, it does not depend on the potential productivity of the prey.

To find the predator isocline for the HVH model, we just have to set:

d P / d t = e \frac{α N P^{- m}}{1 + α h N P^{- m}} P - q P = 0

hence, besides the trivial isocline

P = 0

, there is a non-zero one:

P^{m} = \frac{α (e - h q)}{q} N

This function is positive if and only if

e - h q

is positive. If this were negative, the predator would always decrease and could never persist. Biologically, it means that the handling time and the death rate (i.e., the losses of time and energy) must not be too large compared with the ecological efficiency. The quantity

e - h q

can be considered the ‘net efficiency’ of the metabolism of the predator. In the rest of the paper, we shall always assume that it is positive,

e - h q > 0

(7)

Fig. 1a–e shows the shape of the predator isocline $({Iso}_{P})$ in the phase plane $(N, P)$ :

• if $m = 0$ , the predator dependence disappears from the expression for the functional response and the isocline is a vertical line, as in all prey-dependent models (Fig. 1a);
• if $0 < m < 1$ , the isocline is a convex curve with a horizontal slope at its origin (Fig. 1b);
• if $m = 1$ (ratio-dependent model), the isocline is a straight line that passes through the origin (Fig. 1c–d);
• if $m > 1$ , the isocline is a concave curve with vertical slope at the origin (Fig. 1e).

Fig. 1
Phase plane portraits with Malthusian prey growth. ${Iso}_{N}$ and ${Iso}_{P}$ are prey and predator isoclines respectively, the arrows indicate the vector field determined by the model on the (N,P) phase plane: (a) the prey-dependent case (m=0), (b) the low-interference case m<1, (c) the ratio-dependent case with explosion, (d) the ratio-dependent case with extinction, (e) the high-interference case (m>1) and (f) the BDA model.

For the BDA model, the predator isocline is always a straight line (Fig. 1f), since:

d P / d t = e \frac{a N}{1 + a h N + a w P} P - q P = 0

implies:

P = \frac{e - h q}{w q} N - \frac{1}{a w}

As in the HVH model, condition (7) is also necessary in the BDA model for the existence of a non-trivial equilibrium.

In the HVH model (except for $m = 0$ ), at any prey density, there are always predator densities low enough that the predator can grow. In the BDA model there is a threshold $N_{C} = q / (a (e - h q))$ , below which the predator cannot grow. This expression is independent of the interference constant w. This contrasts with the HVH model, in which predator growth at low prey density becomes easier when the interference parameter increases. Note that increasing the interference constant has a negative effect on the domain in which the predator increases, by decreasing the slope of the predator isocline.

3 Remarks on stability

The local stability of an equilibrium can be studied by analysing the Jacobian J of the system, i.e., the matrix of partial derivatives of (1), calculated at equilibrium. In practice, a stable non-trivial equilibrium (with sufficiently large basin of attraction) implies the persistence of both populations. The necessary and sufficient conditions for local stability (Routh–Hurwith criteria) in two dimensions are:

Tr J < 0; \det J > 0

A general expression for the matrix J, valid for any prey growth function, can be derived (Appendix A). An interesting and classical result that can be deduced from this expression is that, except when $m > 1$ in the HVH model, a sufficient condition for stability is that the isoclines intersect in a portion of the prey isocline where it has a declining slope. It must be pointed out that this condition is not necessary in general. Only in a prey-dependent model is this condition both necessary and sufficient. Thus, when the functional response depends on P, the equilibrium may be stable even if the prey isocline has a positive slope.

Biologically, the fact that the prey isocline is decreasing means that with more prey, a smaller predator density is sufficient to control the prey. This is due to density dependence in the prey, and generates a negative feedback that ensures stability. The lack of density dependence in Malthusian growth is the reason why the prey isocline is always increasing, except for $m > 1$ (see Fig. 1a–d).

The stability of the trivial equilibria $(0, 0)$ and $(K, 0)$ is also important. It characterises the ability or inability of the species to colonise the system. Generally, the isocline patterns allow us to study stability graphically. In our case, the predator isocline shows that in both the BDA and the HVH models (for $0 ⩽ m < 1$ ), the predator can never become established when the prey density is too small. In the BDA model, this is the threshold effect, and in the HVH model with $0 < m < 1$ , the slope at the origin is horizontal, so the region where the predator increases becomes infinitely small when the prey density tends toward zero.

These general properties will facilitate the study of the different cases below, and will help to determine the specificity of each.

3.1 Malthusian growth

In this section, we examine the case in which the prey growth rate is density-independent (i.e., constant). This Malthusian growth is important to study, because it permits to identify the cases in which stability is solely due to control by the predator. This is particularly important in a biological control perspective, since the prey (the pest) must be maintained at very low densities, at which the density dependence modelled by the logistic model does not operate.

3.1.1 HVH model

The prey equation in the HVH model reads:

d N / d t = r N - \frac{α N P^{- m}}{1 + α h N P^{- m}} P

Fig. 1a–e shows the shape of the isoclines for various values of m. The details of the analysis are given in Appendix B.

If $e - h q > 0$ , there exists a unique non-trivial equilibrium, except when $m = 1$ . It is given by:

N^{*} = \frac{q}{e r} {(\frac{e r}{α (e - h q)})}^{1 / (1 - m)}

P^{*} = {(\frac{e r}{α (e - h q)})}^{1 / (1 - m)}

When

m = 1

(the ratio dependent case), the only equilibrium is

(0, 0)

(Fig. 1c and d), and according to the values of the parameters and of the initial conditions, there is either explosion of both densities (the predator does not control the prey and ‘follows’ the prey to explosion) or convergence to

(0, 0)

[22,23] (the predator consumes all the prey and then dies of starvation, which is the classical laboratory dynamics [24,25]).

A necessary and sufficient condition for stability (see Appendix B for details) is:

0 < \frac{r h}{e - h q} < m < 1

(8)

It is a classical result that the Holling functional response (i.e.,

m = 0

) with a Malthusian growth for the prey leads to an unstable equilibrium. This is because the prey isocline is an increasing curve and the predator isocline is vertical (Fig. 1a). If

m > 0

, interference can stabilise the equilibrium, provided that r is not too large (Figs. 1b and 3b). Moreover, there is always an upper limit for h, under which we have stability: reducing the handling time contributes to stability by facilitating prey growth control. In particular, when

h = 0

, which corresponds to the HVLV model, the equilibrium always exists and is always stable.

Fig. 3
Existence and stability of a positive nontrivial equilibrium in relation to growth rate r and interference parameter w (BDA model, (a) and (c)) or m (HVH model, (b) and (d)). Figures (a) and (b) are for Malthusian prey growth, (c) and (d) for logistic prey growth. The model parameters are K=70; a=α=0.5; h=0.1; e=0.2; q=0.9.

It should also be emphasised that for $m > 1$ (Fig. 3b), the equilibrium is never stable. According to the initial conditions, there is either extinction or explosion, while when $0 < m < 1$ , interference has a stabilising influence, the interference parameter becomes destabilising when $m > 1$ . This result had already been found by Hassell and May [9] for discrete time models.

It is also noteworthy that the effect of m on the values of equilibrium densities depends on the parameters. In particular, if $e r > α (e - h q)$ , the equilibrium densities rise with m, when $0 < m < 1$ , and decrease when $m ⩾ 1$ . If $e r < α (e - h q)$ , it is the opposite.

3.1.2 BDA model

For the BDA model, the prey equation is:

d N / d t = r N - \frac{a N}{1 + a h N + a w P} P

The isoclines are two lines (Fig. 1f) that intersect at positive densities if and only if

w < \frac{e - h q}{e r}

(9)

(note that this is equivalent to the slope of the predator isocline being superior to the slope of the prey isocline). If this equation is not fulfilled, then there is no non-trivial equilibrium and both densities explode. Thus, contrary to the HVH model, there is an upper limit for the interference parameter beyond which there is no non-trivial equilibrium. In common with the HVH model, a too intense interference between predators leads to explosion (unstable equilibrium, the predators are not able to control the prey), as does a too-high prey growth rate.

If the equilibrium exists, it is given by

N^{*} = \frac{q}{a (e - h q - e r w)}

P^{*} = \frac{e r}{a (e - h q - e r w)}

Raising the parameter w increases these values. Unlike the HVH model, this tendency does not depend on the other parameters. The stability criterion is:

\frac{h}{e} < w < \frac{e - h q}{e r}

(10)

Thus, interference w and efficiency e are stabilising factors, whereas the handling time h has a destabilising influence. Note that the upper limit on the interference parameter beyond which the equilibrium is unstable coincides with the upper limit for existence of the equilibrium.

In sum, the same qualitative conclusion about the effects of the parameters on stability emerges from both models: the interference parameter is stabilising up to a certain limit, beyond which there is instability or even disappearance of the positive equilibrium. Parameters r, h and q are destabilising, while e is stabilising. Moreover, a (or α) has no effect on stability in either model. A difference on the effect of the interference parameter is that, in the BDA model, increasing the interference parameter rises the equilibrium densities when it is under the threshold value $w = (e - h q) / e r$ , while in the HVH model, the effect depends on the parameter values.

3.2 Logistic growth

In the case of logistic growth, the system cannot explode because the prey cannot grow beyond its carrying capacity K. If an equilibrium is stable in the logistic case and not in the Malthusian case, its stability is basically due to the limitation of the prey by its nutrients, and neither by the control of the prey by the predator, nor by feedback on the predator due to interference.

3.2.1 HVH model

The prey growth equation is now:

\frac{d N}{d t} = r N (1 - N / K) - \frac{α N P^{- m}}{1 + α h N P^{- m}} P

Fig. 2a–f shows the different shapes of the isoclines. The details of the analysis are given in Appendix C.

Fig. 2
Phase plane portraits with logistic prey growth: (a) is the prey-dependent case, (b) and (c) are low-interference cases (m<1) without or with a hump in the prey isocline, (d)–(f) are high-interference cases with either two or no non-trivial equilibria, (g) and (h) the BDA model, in the latter case the ‘limited-predation’ case.

The case $m = 0$ is the classical model by Rosenzweig and MacArthur [26] (for its analysis see, e.g., [27]). In summary, there is a non-trivial (positive) equilibrium if:

\frac{q}{a (e - h q)} < K

(11)

Otherwise, the predator cannot persist and

(K, 0)

is a global attractor. This equilibrium is stable if and only if:

\frac{q}{a (e - h q)} < K < \frac{(e + h q)}{a h (e - h q)}

(12)

which is equivalent to say that the equilibrium is on the decreasing branch of the prey isocline (see Fig. 2a). If the equilibrium is stable, it is a global attractor. Otherwise, the global attractor of the system is a stable limit cycle that emerges from the equilibrium as a result of a Poincaré–Andronov–Hopf bifurcation at the critical value

K = (e + h q) / a h (e - h q)

.

For $0 < m < 1$ , there is one non-trivial equilibrium (see Fig. 2b–c). When K goes to infinity, it tends to a limit identical to the Malthusian growth one. If condition (8) is fulfilled, then the equilibrium is stable. If not, there is a limit value $K_{0}$ such that for $K < K_{0}$ it is stable and for $K > K_{0}$ it is not. Again, if the equilibrium is stable, it is a global attractor. Otherwise, there is a stable limit cycle that becomes a global attractor. Note that for $h = 0$ , which corresponds to the original Hassell–Varley model, we always have stability (see Appendix C for details).

For $m = 1$ , besides equilibria $(0, 0)$ and $(K, 0)$ , there is a non-trivial equilibrium if

α (e - h q) < r e

The equilibrium densities are proportional to K and are given by:

N^{*} = K (1 - \frac{α (e - h q)}{r e})

P^{*} = N^{*} \frac{α (e - h q)}{q}

The non-trivial equilibrium is stable if and only if:

r > \frac{α (e - h q)}{e} (1 + q (\frac{h}{e} - \frac{1}{α}))

This condition is independent of K. Moreover,

r > α

is a sufficient condition of stability (in the Malthusian case,

r > α

always implies explosion). In the ratio-dependent case, the prey is controlled by its carrying capacity only. It is somewhat puzzling that in this case r is a stabilising factor, which seems contradictory with the case

0 < m < 1

with a Malthusian growth. This is to be discussed in the next section.

If the equilibrium is stable, stability is not always global [22,23]. According to initial conditions, there is either convergence to $(0, 0)$ (when the initial prey density is too small compared to the initial predator density), or to the non-trivial equilibrium (see Fig. 5.3 in [22]). Note that if there is a stable limit cycle in the ratio-dependent model, it is not always a global attractor either [22,23]. Thus, depending on the model parameters and on the initial values, the instability of the non-trivial equilibrium leads either to species extinction or to sustained oscillations.

For $m > 1$ , there is a critical value $K_{1} (m)$ such that for $K < K_{1} (m)$ there is no non-trivial equilibrium (see Fig. 2e) and for $K > K_{1} (m)$ there are two non-trivial equilibria (Fig. 2d and f). The greater prey density at equilibrium asymptotically tends to K when $K \to + \infty$ , and the smaller decreases to the Malthusian value (see Appendix C for details). The lower equilibrium is always unstable and the higher always stable.

When $K < K_{1} (m)$ , $(0, 0)$ is a global attractor and the whole system becomes extinct. When $K > K_{1} (m)$ , according to initial conditions, the dynamics converge either to the higher equilibrium, or to $(0, 0)$ (in particular when prey density is too low, see Fig. 2d and f).

Finally, some words should be said about deterministic extinction in the HVH model. The phase portraits show that both populations can go extinct for certain initial conditions when $m ⩾ 1$ , but never for $m < 1$ . The phenomenon of deterministic extinction was first studied in the context of the ratio-dependent model [22,28,29] but criticised for not being ‘generic’ (i.e., using the intermediate original BDA model (5), deterministic extinction only occurs for the limiting case $w = {(a P_{0})}^{−1})$ . Our analysis of the HVH model (which is another intermediate model) shows that this phenomenon is ‘half-generic’ in the sense that it happens for all $m ⩾ 1$ . However, the important observation is that for values of $m < 1$ but close to 1, the dynamics can come very close to the axes, and extinction may occur due to demographic or environmental stochasticity.

3.2.2 BDA model

Now the prey growth equation is:

\frac{d N}{d t} = r N (1 - N / K) - \frac{a N}{1 + a h N + a w P} P

The isoclines of this model are drawn in Fig. 2g and h. Again, the details of the analysis can be found in Appendix C. A necessary and sufficient condition to have a non-trivial equilibrium is:

\frac{q}{a (e - h q)} < K

This condition is identical to (11) in the HVH model for

m = 0

. If it is not fulfilled, the dynamics converge to

(K, 0)

.

The non-trivial equilibrium is given by:

N^{*} = \frac{K}{2 e r} (e r - \frac{e - h q}{w} + \sqrt{{(e r - \frac{e - h q}{w})}^{2} + \frac{4 e r q}{K w a}})

P^{*} = \frac{1}{w a q} (a (e - h q) N^{*} - q)

If

w < (e - h q) / (e r)

(which is also condition (9) of the equilibrium positivity in the Malthusian case), then

N^{*}

tends to a limit when

K \to + \infty

:

N_{\infty}^{*} = \frac{q}{a (w e r - e + h q)}

If

w > (e - h q) / e r

then, for high K values:

N^{*} \approx (1 - \frac{e - h q}{e r w}) K

(13)

Thus, when the condition for existence of an equilibrium in the Malthusian case is fulfilled, it corresponds to the limiting case when

K \to + \infty

. The boundary value

w = (e - h q) / (e r)

is the equivalent of

m = 1

in the HVH model: beyond it, the prey is no longer predator-controlled but resource-controlled.

Note that when $w \to + \infty$ , $P^{*} \to 0$ (because N is bounded with respect to w): as in the HVH model, when interference is too large, the predator density at equilibrium is so low that extinction is almost certain.

The non-trivial equilibrium is stable if and only if:

q (\frac{h}{e} - w) - (1 + q (\frac{h}{e} - w)) N^{*} / K < 0

In particular, if

w > h / e

, it is always stable, because

1 - N^{*} / K > 0

. If the equilibrium is stable, it is a global attractor. Otherwise, there is a stable limit cycle.

The same type of conclusions about both equilibrium and stability arise in both models. Increasing interference is at first stabilising, then it prevents the predator from controlling the prey and, eventually, produces the extinction of the predator, even if, mathematically, we have a stable equilibrium. There is a difference concerning the existence of a non-trivial equilibrium. In the HVH model, when $m < 1$ , we always have an equilibrium, no matter what the values of the parameters are (provided that $e - h q > 0$ ). In the BDA model, on the contrary, there is a critical value of K under which the predator cannot persist because of insufficient prey supply. However, in practice, when m is small, the predator density at equilibrium quickly tends to zero when K decreases, and we may consider that we should have extinction because the conditions of validity of the continuous deterministic models are not fulfilled anymore. Therefore, this difference between the models is not qualitatively significant. Another feature of the BDA model is that it yields conditions of existence of the non-trivial equilibrium quite similar to the HVH model when $m = 0$ .

4 Effects of K, r, m, w on stability

Let us consider a predator with fixed values of a (α in the HVH model), h, q and e. These can be considered as intrinsic parameters that are assumed not to depend on environmental factors. We will focus now on the effects of the crucial parameters K, r, m, w on stability. Parameters K and r relate to the growth of the prey, and can be controlled experimentally (for instance by modifying the input of nutrients) or influenced by the environment (for instance by temperature). Parameters m and w have the same phenomenological function, they represent the degree of interference in the two models HVH and BDA, respectively. They can be considered to be intrinsic to the predator, and their effect must be studied because interference is the topic of this paper.

The Paradox of Enrichment [19] is based on the mathematical finding that increasing K has a negative effect on stability in many models, which contradicts natural observations. On the other hand, interference between predators is thought to be a stabilising factor, because it reduces predator reproduction as predators become more numerous and thus acts as a negative feedback. The effect of the maximal prey growth rate r on stability is often thought to be destabilising, because under high r value, the prey is able to escape the control of predation. However, in the simple logistic equation, increasing r makes only the convergence to the carrying capacity K faster. Thus, increasing r should have a stabilising influence in the neighbourhood of the carrying capacity K.

Fig. 3 shows the shape of the domains of stability in the plane $(r, m)$ for the HVH model (Fig. 3b and d), and in the plane $(r, w)$ for the BDA model (Fig. 3a and c).

In the HVH model, condition (8) provides a straightforward way of getting the domains of stability in the Malthusian case (Fig. 3b). In the logistic case, the limits of the domains are intractable analytically. Fig. 3d shows how these domains, obtained numerically, look like (see also Appendix D). There is a stable domain for small r, and another one for large r. The interpretation is that a large maximum growth rate ensures a fast return to equilibrium as in the simple logistic equation (for large r, the equilibrium is close to K: the prey is controlled by its carrying capacity and not by the predator). For small r, the dominating phenomenon is the stabilising effect of the interference that controls the predator. The prey is controlled by the predator, and a too large value of r prevents this control. When $m > 1$ , there is a threshold value of r below which there is no equilibrium, and above which there are two equilibria, one unstable, and one stable.

In the BDA model, there is no stability domain for small r values. Increasing parameter r has a stabilising influence. The domains are drawn in Fig. 3c (see also Appendix D).

An important remark must be made about the positive effect of r on stability. Classically, we use the logistic form for the prey growth rate assuming that the time delay T of the feedback of density on reproduction is negligible. However, this is only true when the product rT is small. If r is large, this amplifies the effect of the time delay, and instability around K might appear [30] even for very small values of T. This stabilising effect of r can thus be seen as a mathematical artefact.

The effect of K on stability is essentially the same in both models. Either the equilibrium is stable for all K, or there is a limit on K beyond which it becomes unstable (see Fig. 4a and b). The results for the Malthusian case may be deduced just by setting $K = + \infty$ . There is no discontinuity at infinity for K. Fig. 4 shows the shape of the stability portraits in the plane $(K, m)$ for the HVH model, and in the plane $(K, w)$ for the BDA model.

Fig. 4
Existence and stability of a positive nontrivial equilibrium in relation to carrying capacity K and interference parameter w (BDA model with parameter values a=2.0; h=0.5; e=0.3; q=0.3; r=1.0, (a)) and m (HVH model with parameter values α=0.5; h=0.1; e=0.2; q=0.9; r=1.0, (b)). The dashed line in (b) is the asymptotic limit of the curved line separating the stable and unstable regions.

In the HVH model, when $m < 1$ , there is always stability for small K, and depending on the values of the parameters, either there is an upper limit on K for all m or there is a region $m_{C} < m < 1$ in which there is always stability (with $m_{C} = r h / (e - h q)$ , see Fig. 4b). When $m > 1$ , there is a threshold on K below which there is no equilibrium, and above which there are two equilibria, one unstable, and one stable.

In the BDA model, we have a lower limit on K for existence of a non-trivial positive equilibrium. For $w > h / e$ , we always have stability. For smaller w, there is an upper limit on K for stability (similar to condition (12) when $w = 0$ ). One can define a curve $w (K)$ of the lower limit for stability for a given K. It is an increasing curve with a limit when $K \to + \infty$ which is less than $h / e$ (Fig. 4a).

In conclusion, both models lead to qualitatively identical results. The most important difference between the models is the effect of the interference parameter on stability. In the BDA model, there is always a value of the interference parameter above which the Paradox of Enrichment disappears. In the HVH model on the contrary, the existence of the discontinuity at $m = 1$ implies that a too strong interference leads to instability of the dynamics, at least for small K, whereas, for values $m < 1$ , interference is stabilising, as in the BDA model. However, condition (13) implies that, with a too great interference parameter w, the densities at equilibrium are proportional to the carrying capacity K. Even if this has no effect on stability, it shows that even in the BDA model, a large interference makes the predator unable to control the prey.

5 Discussion

The two general models for interference studied here have very different mathematical expressions and very different conceptual foundations. Despite these differences, it is quite satisfactory that they give essentially the same qualitative results concerning the effects of the parameters on the stability of the system. The interference parameter has a peculiar characteristic: when it is low, increasing it has a positive effect on the stability and thus on the persistence of the system. When it is large, it is the contrary (with logistic growth, increasing the interference parameter ensures stability, but leads to very small predator densities).

Another important point is that the Paradox of Enrichment [19] still holds in both types of model for small values of the interaction coefficients. The observation that diversity decreases in highly productive sites in both experimental manipulations [31] and in situ recordings (e.g., [32,33]) has led to many theories aiming at explaining it [34]. Rosenzweig [34] argues that the effect of the Paradox of Enrichment for predation (and of its equivalent for competition [35]) acts on a short time scale before natural selection occurs, and could provide an alternative to the more popular interpretation, based on Tilman's competition theory [31].

Is it relevant to take interference into account in the models? Which values of interference are to be expected in nature? Since the very beginning of interference theory, it was pointed out that the values found in the laboratory could not be extrapolated to field populations [36,37] because densities used in the laboratory were far higher than in the field. However, some field studies reported values of the parameter m higher than 0.5 and even higher than 1 (examples are given in [18,38,39]). Free et al. [3] suggested that this was not due to direct interference, but was rather a consequence of the tendency of the predators to aggregate in patches with high prey density. As mentioned in the Introduction, we disregard in this paper the precise mechanisms leading to interference, but we can conclude that there are field situations where we need a high interference parameter to explain the dynamics.

The question of which of the models studied here should be used in practice must be addressed. As mentioned in the introduction, the HVH model was proposed after empirical evidence of such a form of the functional response. It is widely used in the study of interfering populations, from insects to vertebrates [40]. In a direct comparison, it won over the BDA model describing soil fauna [38]. However, up to now, this form was rarely incorporated into dynamic models at the level of populations. No easy empirical test can be made for discriminating between the two classes of models studied here, but varying the nutrient input, for instance, gives a way of identifying the predator isocline and hence, indirectly, the functional response and the interference parameter. Indeed, the values of the densities at equilibrium, plotted in the phase space, must all be on the predator isocline. This method can be applied when it is difficult to measure predator consumption directly. An advantage of the HVH model is that the predator isocline is sufficient to determine the interference parameter, while in the BDA model, one needs to identify other parameters to measure the interference parameter.

Other considerations than experimental ones can also help determine which model is the most proper to use. An argument that is raised against the Hassell–Varley and ratio-dependent models is that they do not rest on mechanistic arguments. However, it is interesting to note that in the modelling of bioreactors (that describe processes that are close to predator–prey interaction: the prey is the substrate, the predator the microorganism), no mechanical approach has ever been used to find a form for the functional response [41,42]. On the contrary, the approach has always been entirely empirical, even though the processes should be much simpler than in a predator–prey system. In fact, the HVH model is satisfactory for several theoretical reasons. An advantage of using a functional response of the form $g (N P^{- m})$ is that, with a small number of assumptions about the shape of the function g, it is possible to obtain interesting qualitative predictions. Moreover, the parameter m can be easily interpreted and has been directly measured in various experiments (examples in [5,18,40,43]). Some authors pointed out that the Hassell–Varley model did not accurately represent the results of some experiments: they seldom yield straight lines in the $\log N$ – $\log P$ plane (e.g., [3]). However, the BDA model does not do any better here. Another advantage of the HVH model is that there is a fixed boundary $m = 1$ that corresponds to the value of the interference parameter above which the predator can never control the prey. In the BDA model, on the contrary, this value depends on the other parameters (condition (10)). Moreover, the value $m = 1$ corresponds to the ratio-dependent model, which can naturally be interpreted as the boundary over which interference leads to overcompensation (increasing predator density leads to a loss of per capita consumption more than proportional).

In both HVH and BDA models, the ratio-dependent model is obtained with a high interference. It seems logical that a ratio-dependent functional response would arise when interference between predators is strong (i.e., the quantity of prey eaten per predator only depends on the relative abundance of prey and predators). Indeed, a ratio-dependent consumption means that there is prey sharing among predator individuals. This phenomenon is directly related with interference (see preceding paragraph), which prevents each predator to behave as if there were no other predators, which is the case in prey-dependent models.

We will conclude with a few evolutionary speculations suggested by the effects of the parameters on stability. Regarding the predator parameters other than interference, the interest of both the individual and of the population consists in having high searching efficiency, high conversion efficiency, low handling time, and low mortality. All these parameters have monotonous effects on predator population fitness. This is definitely not the case for the interference parameter. Both models studied here predict that both too high and too low interference should lead to instability or extinction of the population. However, the interest of individuals is always to increase its interference with the other predators, e.g., being more aggressive to obtain more food, or aggregating in patches with high prey density [3]. The results obtained here imply that this tendency would eventually lead the predator population to collapse. Although some studies giving estimates of the parameter m using the Hassell–Varley formalism report values of m larger than 1, our analysis suggests that these situations should be rare in the field: populations with an interference parameter larger than 1 are at risk of becoming unstable if the carrying capacity of the prey is too low (e.g., in a ‘bad year’). Moreover, predator equilibrium densities rapidly decrease when m rises. Given that for $m < 1$ , condition (8) shows that the interaction between predators is a stabilising factor, one could thus argue that evolution has selected species of predators that have values of m close to, but smaller than 1. This could be an explanation of why Arditi and Akcakaya [5] found that, in many sets of experimental data, m was close to 1.

Finally, we have seen that in both models, below large values of the interaction parameter (m or w), the prey density is proportional (or nearly proportional) to K. This means that, as individual selection tends to favour high values of this parameter, control by the lower trophic level ought to be common. However, it can be physically impossible to achieve high interaction, as it is the case for many herbivorous zooplankton species, which are known to have a significant effect on phytoplankton density. Parasites are also known to have a significant effect on the control of some populations, and they are likely to have a small interference parameter. Thus, predators that can achieve a high (direct or indirect) interference should not control their prey: this implies bottom-up regulation. Predators that cannot interfere much should lead to top-down regulation.

Appendix A Derivation of the Jacobians and of the local stability criterion

A.1 HVH model

For the HVH model, the dynamical equations are of the form:

d N / d t = f_{N} (N, P) = F (N) - g (N P^{- m}) P

d P / d t = f_{P} (N, P) = e g (N P^{- m}) P - q P

The Jacobian of the system is given in Fig. 5. We will use that at the non-trivial equilibrium, we have $g (N P^{- m}) = q / e$ because $d P / d t = 0$ .

Moreover, using the type II (saturated) form of g,

g (θ) = \frac{α θ}{1 + α h θ} (with g^{'} (θ) = \frac{α}{{(1 + α h θ)}^{2}})

we get at equilibrium

N P^{- m} = q / α (e - h q)

. Since we are only interested in positive equilibria, this implies:

e - h q > 0

. Together with

g^{'} (N P^{- m}) = \frac{α}{{(1 + h α N P^{- m})}^{2}} = \frac{1}{α {(N P^{- m})}^{2}} {(g (N P^{- m}))}^{2} = \frac{(α {(e - h q))}^{2}}{α e^{2}}

the Jacobian at the non-trivial equilibrium (if it exists) is:

J = [\begin{matrix} F^{'} (N) - P^{1 - m} \frac{(α {(e - h q))}^{2}}{α e^{2}} & - \frac{q}{e} [1 - m + m \frac{h q}{e}] \\ P^{1 - m} \frac{(α {(e - h q))}^{2}}{α e} & - m \frac{q}{e} (e - h q) \end{matrix}]

(A.1)

For m> 1 this Jacobian has the sign structure:

[\begin{matrix} * & < 0 \\ > 0 & < 0 \end{matrix}]

From the implicit function theorem, we know that the slope of the prey isocline (expressed with P as a function of N) is given by $- (\partial f_{N} / \partial N) / (\partial f_{N} / \partial P)$ and since $(\partial f_{N} / \partial P) < 0$ , a negative slope ensures that $(\partial f_{N} / \partial N) < 0$ , which in turn results in a sign structure with $Tr J < 0$ and $\det J > 0$ . Thus, a negative slope of the prey isocline guarantees stability of the non-trivial equilibrium.

Note that this result can be generalised for any system where g has the form $g (N P^{- m})$ provided that $0 < m < 1$ , and that g be an increasing concave function with $g (0) = 0$ and $0 < g^{'} (0) < \infty$ . To see this, all we have to do is to find the sign of:

\partial f_{N} / \partial P = m N P^{- m} g^{'} (N P^{- m}) - g (N P^{- m})

The derivative of the function $x \mapsto m x g^{'} (x) - g (x)$ is $m x g^{″} (x) - (1 - m) g^{'} (x)$ . If we assume that g is concave (i.e. $g^{″} < 0$ ) and increasing (i.e. $g^{'} > 0$ ) and $0 < m < 1$ , then this latter expression is negative. The function $x \mapsto m x g^{'} (x) - g (x)$ thus takes negative values for $x > 0$ , because it is a decreasing function that takes the value 0 at $x = 0$ . In conclusion, if the prey isocline decreases at the equilibrium, then $(\partial f_{N} / \partial N) < 0$ and thus the stability criterion is met.

A.2 BDA model

Now the Jacobian is (see Fig. 6).

Using that at equilibrium

\frac{a N}{1 + a w P + a h N} = \frac{q}{e} and

\frac{a P}{1 + a w P + a h N} = F (N) / N

we get:

J = [\begin{matrix} F^{'} (N) - \frac{e - h q}{e} \frac{F (N)}{N} & \frac{q}{e} (w \frac{F (N)}{N} - 1) \\ (e - h q) \frac{F (N)}{N} & - q w \frac{F (N)}{N} \end{matrix}]

(A.2)

Note that we have $\partial f_{N} / \partial P < 0$ because $w F (N) / N = a w P / (1 + a h N + a w P) < 1$ at equilibrium. The same stability condition as in the HVH model can thus be stated.

Appendix B Derivation of the isoclines and equilibrium properties for Malthusian growth

Let $F (N) = r N$ .

B.1 HVH model

The prey isocline is given by

N = \frac{1}{α h r} (α P - r P^{m})

and the predator isocline by

N = \frac{q}{α (e - h q)} P^{m}

Since we are only interested in systems that have a non-trivial equilibrium, the predator isocline must pass through the positive quadrant of the phase plane. We will thus further assume here that $e - h q > 0$ . Three cases have to be distinguished with respect to the value of m (presented graphically in the $(N, P)$ phase plane, see Fig. 1):

– if $0 ⩽ m < 1$ , the prey isocline implicitly defines a monotonous concave function of $N \in [0, + \infty [$ that increases from ${(r / α)}^{1 / (1 - m)}$ to +∞, while the predator isocline is either a vertical line $(m = 0)$ or a convex increasing function through the origin (Fig. 1a and b);
– if $m = 1$ , both isoclines are straight lines through the origin. $e > h q$ and $α > r$ are required for them to pass through the positive quadrant; two cases may be distinguished according to the slopes of the two isoclines, $e r > α (e - h q)$ and $e r < α (e - h q)$ ; in both cases, there is either extinction or explosion, depending on the initial conditions (Fig. 1c and d);
– if $m > 1$ , both isoclines are essentially non-linear (see Fig. 1e); note that, due to the fact that the prey isocline crosses the P-axis in the origin and in $\sqrt[1 - m]{r / α}$ , there is only a bounded domain inside which the prey decreases.

The shape of the isoclines suggests that there is one and only one non-trivial equilibrium provided that $m \neq 1$ . Analytically, we obtain:

N^{*} = \frac{q}{e r} {(\frac{e r}{α (e - h q)})}^{1 / (1 - m)}

P^{*} = {(\frac{e r}{α (e - h q)})}^{1 / (1 - m)}

Note that the ratio $N^{*} / P^{*}$ is equal to the value $q / (e r)$ and thus does not depend either on α or h.

Fig. 1e suggests that the non-trivial equilibrium is never stable if $m > 1$ . We can see graphically that depending on the initial conditions in the neighbourhood of the equilibrium, a trajectory can either explode to the infinity or converge to $(0, 0)$ , as in the case $m = 1$ .

At the equilibrium, using (A.1), we get:

J = [\begin{matrix} r \frac{h q}{e} & - \frac{q}{e} ((1 - m) + m \frac{h q}{e}) \\ r (e - h q) & - m \frac{q}{e} (e - h q) \end{matrix}]

Thus:

(B.1)

A necessary condition for local stability is

\det J > 0

, i.e.

m < 1

. The condition for the trace

(Tr J < 0)

is equivalent to:

m > \frac{r h}{e - h q} \Leftrightarrow h < \frac{m e}{r + m q}

(B.2)

B.2 BDA model

The non-trivial equilibrium is

N^{*} = \frac{q}{a (e - h q - e r w)}, P^{*} = \frac{e r}{q} N^{*}

and it exists if and only if

w < (e - h q) / e r

. Note that this implies

1 - r w > 0

. The prey isocline, given by

P = \frac{r (1 + a h N)}{a (1 - w r)}

is a straight line with positive slope and is shown in Fig. 1f.

Using (A.2) with $F (N) = r N$ yields the expression of the Jacobian

J = [\begin{matrix} r \frac{h q}{e} & \frac{q}{e} (w r - 1) \\ r (e - h q) & - q w r \end{matrix}]

and the criterion of local stability resolves to

\frac{h}{e} < w < \frac{e - h q}{e r}

Appendix C Derivation of the isoclines and equilibrium properties for logistic growth

Let $F (N) = r N (1 - N / K)$ . We will first discuss the general form of the isoclines and then discuss existence and stability of the nontrivial equilibrium.

C.1 HVH model

Since the results for $m = 0$ are classical, we shall only deal with $m > 0$ .

Since the predator isocline remains the same as in the case of Malthusian growth, we only have to discuss the form of the prey isocline, which is now given by:

\frac{r α h}{K} N^{2} + r (\frac{1}{K} P^{m} - α h) N + α P - r P^{m} = 0

This expression can be rewritten as:

N = φ_{\pm} (P) = \frac{K}{2 r α h} [r (α h - \frac{1}{K} P^{m}) \pm \sqrt{r^{2} {(α h - \frac{1}{K} P^{m})}^{2} - \frac{4 r α h}{K} (α P - r P^{m})}]

(C.1)

Let first $0 < m < 1$ . There are two possibilities:

– ${(α h K)}^{1 - m} < {(r / α)}^{m}$ ; $φ_{-}$ is always negative and the isocline is given by $φ_{+}$ only, which is positive between $P = 0$ and $P = {(r / α)}^{1 / (1 - m)}$ where it decreases from $N = K$ to $N = 0$ (see Fig. 2b);
– ${(α h K)}^{1 - m} > {(r / α)}^{m}$ ; a graphical argument shows that the expression under the square root becomes zero at some value $\bar{P}$ where $φ_{-}$ and $φ_{+}$ are positive and equal, and that corresponds to a maximum of the isocline as expressed as P as a function of N. Fig. 2c shows that case.

Let now $m = 1$ . This case, better known as a ratio-dependent model, has been studied extensively (see main text for the references), so we will only summarise here the results. The predator isocline is always a straight line through the origin. For the prey isocline, two cases arise:

– if $α > r$ then the isocline is a parabola-like curve that passes through the trivial equilibria $(0, 0)$ and $(K, 0)$ . In that case, a non-trivial equilibrium exists if $α (e - h q) < r e$ ;
– if $α < r$ then it is a decreasing curve defined on $[K \frac{r - α}{r}, K]$ that passes through $(K, 0)$ providing the existence of the non-trivial equilibrium.

Finally, let $m > 1$ . Eq. (C.1) still holds, but now we can see that for large P, positive values of N can be defined. To study the shape of the curve, it is more illustrative to make qualitative remarks than to try to study its derivative.

We will refer to the isocline as a set of two functions $N = φ_{-} (P)$ and $N = φ_{+} (P)$ . The following remarks can be made:

– $φ_{-} (0) = 0$ and $φ_{+} (0) = K$ ;
– $φ_{-}$ increases for small P; this is because the term in P under the square root dominates the terms in $P^{m}$ ;
– $φ_{-}$ only crosses the axe $N = 0$ in two points: $P = 0$ and $P = {(α / r)}^{1 / (m - 1)}$ ; beyond this domain, $φ_{-} < 0$ if it exists;
– $φ_{+} < K$ (seen in Eq. (C.1)); it tends to K when P tends to infinity.

No trivial criterion can determine if the expression under the square root of (C.1) becomes negative or not. However, it can be said that it can never happen for $P > {(α / r)}^{1 / (m - 1)}$ , whatever the parameters are. It is easily seen from Eq. (C.1) – a graphical argument is shown in Fig. 3f – that for K sufficiently large it never happens. On the contrary, when ${(α h K)}^{1 / m} < {(α / r)}^{1 / (m - 1)}$ , i.e., for sufficiently small K values, the expression becomes negative in a bounded interval of $] 0, + \infty]$ containing ${(α h K)}^{1 / m}$ . Again, at the endpoints of this interval, $φ_{-} = φ_{+}$ .

Now one can have an idea of the possible shapes of the isoclines in the $(N, P)$ phase plane. Fig. 3e shows the case of low K values when there is no equilibrium, and Fig. 3d the case of high K when there are two equilibria, and expression (C.1) is defined for all P. The intermediate case (two equilibria and a domain where expression (C.1) is not defined) is represented in Fig. 3f.

The predator and the prey equation resolve to

P^{m} = \frac{α e - h α q}{q} N and P = \frac{r}{q} N (1 - \frac{N}{K})

respectively. Combining them gives the equation:

N^{* 2} - K N^{*} + K \frac{q}{e r} {(\frac{α e - α h q}{q} N^{*})}^{1 / m} = 0

(C.2)

This equation is of the form:

ψ (N, K, γ) = N^{2} - K N + K γ β^{1 / m} N^{1 / m} = 0

(C.3)

To study the stability, Eq. (A.1) implies:

J = [\begin{matrix} r (1 - 2 \frac{N^{*}}{K}) - P^{* 1 - m} \frac{{(α e - α h q)}^{2}}{α e^{2}} & - \frac{q}{e} (1 - m + m \frac{h q}{e}) \\ e P^{* 1 - m} \frac{{(α e - α h q)}^{2}}{α e^{2}} & - \frac{q}{e} m (e - h q) \end{matrix}]

Using again the equations above, we get the following relation:

r (1 - 2 N^{*} / K) = 2 r / N^{*} (N^{*} - N^{* 2} / K) - r = 2 \frac{q}{e} {(\frac{α e - α h q}{q})}^{1 / m} N^{* (1 - m) / m} - r

Replacing

P^{* 1 - m}

by

({(α e - α h q) / q)}^{(1 - m) / m} N^{(1 - m) / m}

, we eventually get

J = [\begin{matrix} \frac{q}{e^{2}} {(\frac{α (e - h q)}{q})}^{\frac{1}{m}} (e + h q) N^{* \frac{1 - m}{m}} - r & - \frac{q}{e} (1 - m + m \frac{h q}{e}) \\ \frac{q}{e} {(\frac{α (e - h q)}{q})}^{\frac{1}{m}} (e - h q) N^{* \frac{1 - m}{m}} & - \frac{q}{e} m (e - h q) \end{matrix}]

(C.4)

Case $0 < m < 1$ . Eq. (C.3) can be rewritten as:

\frac{{(N / K)}^{1 - m}}{{(1 - N / K)}^{m}} = \frac{1}{γ^{m} β K^{1 - m}}

Under this form, it is clear there is one and only one equilibrium for each K, and that $N^{*} / K$ is a function of K that decreases from 1 to 0 when K varies from 0 to infinity. That comes from the properties of the function:

θ \in] 0, 1 [\mapsto \frac{θ^{1 - m}}{{(1 - θ)}^{m}}

that increases from 0 to infinity.

Moreover, applying the Implicit Function Theorem to Eq. (C.3) gives:

\frac{d N}{d K} = - \frac{\partial ψ / \partial K}{\partial ψ / \partial N} = \frac{m N^{2}}{K (m N + (1 - m) (K - N))} > 0

Thus,

N^{*}

as a function of K increases. Moreover it tends to a limit that is identical to the value given by Eq. (C.3) when

1 / K = 0

, i.e.,

N^{*} = {(γ β^{1 / m})}^{m / (1 - m)}

.

For $K \to 0$ we have $N^{*} \to 0$ . Now, when $N^{*} \to 0$ in (C.4), we have:

J \to [\begin{matrix} - r & - \frac{q}{e} (1 - m + m \frac{h q}{e}) \\ 0 & - \frac{q}{e} m (e - h q) \end{matrix}]

which corresponds obviously to a stable equilibrium (

Tr J < 0

and

\det J > 0

). We then have TrJ increasing, and

\det J = \frac{q}{e^{2}} (e - h q) ({(\frac{α (e - h q)}{q})}^{1 / m} (1 - 2 m) N^{* (1 - m) / m} q + e m r)

increases

(m < 0.5)

or decreases

(m > 0.5)

monotonically when K (thus

N^{*}

) increases, the other parameters being held constant. This is because they both are linear functions of

N^{* (1 - m) / m}

. The limit values when

K^{*} \to + \infty

are a positive one for

\det J

(which is therefore positive for any K), and

\frac{q}{e} (r h - m (e - h q))

for TrJ. They correspond to inequality (B.1) found in Appendix B. Thus, criterion (B.2) found for stability in the Malthusian case is a sufficient condition that this system does not destabilise with enrichment (increasing K).

Note that, in particular, when ${(α h K)}^{1 - m} < {(r / α)}^{m}$ , i.e., the prey isocline is decreasing as a function of N, then the non-trivial equilibrium is always stable.

Case $m = 1$ (ratio-dependent). The non-trivial equilibrium is given by:

N^{*} = K (1 - \frac{α (e - h q)}{r e})

P^{*} = N^{*} \frac{α (e - h q)}{q}

therefore, it exists if and only if

α (e - h q) < r e

. In particular, when

r > α

, it always exists.

The Jacobian is, according to (A.1):

J = [\begin{matrix} r (1 - 2 N^{*} / K) - \frac{(α {(e - h q))}^{2}}{α e^{2}} & - \frac{h q^{2}}{e^{2}} \\ \frac{(α {(e - h q))}^{2}}{α e} & - \frac{q}{e} (e - h q) \end{matrix}]

Whence

J = [\begin{matrix} \frac{α (e^{2} - h^{2} q^{2})}{e^{2}} - r & - h \frac{q^{2}}{e^{2}} \\ \frac{α {(e - h q)}^{2}}{e} & - \frac{q}{e} (e - h q) \end{matrix}]

Calculating TrJ and $\det J$ , gives:

Tr J = - r + \frac{α (e^{2} - h^{2} q^{2})}{e^{2}} - \frac{q}{e} (e - h q) = - r + \frac{α (e - h q)}{e} (1 + q (\frac{h}{e} - \frac{1}{α}))

\det J = \frac{q}{e} (e - h q) | \begin{matrix} \frac{α (e^{2} - h^{2} q^{2})}{e^{2}} - r & - h \frac{q^{2}}{e^{2}} \\ \frac{α (e - h q)}{e} & - 1 \end{matrix} | = (r - \frac{α}{e} (e - h q)) q \frac{(e - h q)}{e}

Since the criterion for the positive equilibrium to exist is $r e > a (e - h q)$ , $\det J > 0$ , and the criterion for stability reads:

r > \frac{a (e - h q)}{e} (1 + q (\frac{h}{e} - \frac{1}{a}))

In particular, one can see that $e > h a$ implies stability, because this inequality follows then simply from the existence condition for the equilibrium.

Case $m > 1$ . It is better to rewrite Eq. (C.3) as:

{(N / K)}^{m - 1} {(1 - N / K)}^{m} = \frac{γ^{m} β}{K^{m - 1}}

(C.5)

Now, the function $θ \in] 0, 1 [\to θ^{m - 1} {(1 - θ)}^{m}$ increases from 0 to a maximum in $] 0, 1 [$ , and then decreases to 0. Thus for small K, there is no non-trivial equilibrium, and beyond a critical value that depends on the other parameter there are two equilibria. For $K \to + \infty$ , the larger equilibrium tends to infinity, while the smaller one tends to ${(γ^{m} β)}^{1 / (m - 1)}$ , which is simply the equilibrium in the case of Malthusian growth.

Fig. 2e, f, and d show the evolution of the isoclines when K varies ( $K = 6, 8$ and 10, respectively). One can see on these figures that, in the cases of large K values, i.e., when two non-trivial equilibria exist, the smaller equilibrium is unstable, and the higher one is stable.

C.2 BDA model

Let us first discuss the shape of the prey isocline. The non-trivial isocline for the prey is given by:

P = \frac{r (1 - N / K) (1 + a h N)}{a (1 - r w) + a r w N / K}

As for the case $m = 1$ in the last model, two cases arise:

– $r w < 1$ : we get a parabola like curve that reaches its maximum between 0 and K and that becomes negative for $N > K$ .
– $r w > 1$ : we get a decreasing curve that is positive in $[K \frac{r w - 1}{r w}, K]$ .

The shape of the isoclines suggests that a necessary and sufficient condition to have a non-trivial equilibrium is that the predator isocline (the same as in the case of Malthusian growth) passes through the N-axis between 0 and K, that is:

\frac{q}{a (e - h q)} < K

(C.6)

We have:

r N^{*} (1 - N^{*} / K) = \frac{q}{e} P^{*}

(C.7)

and

P^{*} = \frac{1}{w a q} (a (e - h q) N^{*} - q)

Substituting $P^{*}$ into Eq. (C.7), we get:

\frac{e r}{q K} N^{* 2} + (\frac{e - h q}{w q} - \frac{e r}{q}) N^{*} - \frac{1}{w a} = 0

(C.8)

which has two solutions, but only one is positive under condition (C.6):

N^{*} = \frac{K}{2 e r} (e r - \frac{e - h q}{w} + \sqrt{{(e r - \frac{e - h q}{w})}^{2} + \frac{4 e r q}{K w a}})

(C.9)

There are two possibilities:

– $e r < (e - h q) / w$ , then $N^{*}$ has a limit when $K \to + \infty$ ; de l'Hopital's rule yields:
$N_{\infty} = \frac{q}{w a \sqrt{{(e r - \frac{e - h q}{w})}^{2}}}$
– $e r > (e - h q) / w$ (which is always the case when $r w > 1$ ); then $N^{*}$ diverges and we have:
$N^{*} \approx (1 - \frac{e - h q}{e r w}) K$

The Jacobian at equilibrium is, according to Eq. (A.2):

J = [\begin{matrix} r (- N / K + \frac{h q}{e} (1 - N / K)) & \frac{q}{e} (w r (1 - N / K) - 1) \\ r (e - h q) (1 - N / K) & - q w r (1 - N / K) \end{matrix}]

Thus, the stability conditions are:

(C.10)

Thus,

Tr J < 0

is a necessary and sufficient criterion for stability, because substituting

N^{*}

(Eq. (C.9)) into the expression of

\det J

shows that the latter is always positive. Moreover, taking into account that

N^{*} / K < 1

, one can see that a sufficient condition for stability is

h < e w

.

Appendix D Comments to the stability portraits

D.1 HVH model

The effect of K has already been studied when the relationship of the logistic case with the Malthusian case was discussed.

To get the shape of the domains of stability in the logistic case, only qualitative remarks can be done. Eq. (C.4) shows that, for large r, the non-trivial equilibrium is stable. This is because for large K, the prey densities are bounded, and thus the Jacobian is equivalent to:

[\begin{matrix} γ - r & β \\ σ & δ \end{matrix}]

for certain constants

γ, β, σ, δ

. On the opposite, when r is small, there is a region qualitatively similar to the one of Malthusian case for small values of m.

When $m > 1$ , Eq. (C.5) shows that, as for the case of parameter K, there is a boundary under which there is no non-trivial equilibrium, and over which there are two non-trivial equilibria, one unstable, and the other stable.

D.2 BDA model

Two cases must be distinguished when studying the effect of K on stability. First, when $e r < (e - h q) / w$ , then $N^{*} / K$ decreases from 1 to 0 when K increases from $q / (a (e - h q))$ to +∞. Thus, TrJ varies monotonously from −r to $r q (h / e - w)$ . Now we have two cases:

– if $h < w e$ then we always have stability (we already knew that);
– if $h > w e$ , then there is an upper limit on K for stability.

In the second case, when $e r > (e - h q) / w$ , TrJ varies monotonously from −r to $r [q (\frac{h}{e} - w) \frac{e - h q}{e r w} - \frac{e r w - e + h q}{e r w}]$ . According to the parameter values, there is either always stability (and $h < w e$ is a sufficient condition) or an upper limit on K for stability.

The effects of r are somewhat simpler to understand. For large r, condition (C.10) is always fulfilled (because (C.9) implies $N^{*} \to K$ when r tends to infinity).

When r tends to zero, de l'Hôpital's rule applied to (C.9) shows that N has the limit $q / (a (e - h q))$ , which is independent of w. Thus, when K is close to $q / (a (e - h q))$ , condition (C.10) is always fulfilled and when K is large, it is only fulfilled when $w > h / e$ .

References

[1] V. Volterra Leçons sur la théorie mathématique de la lutte pour la vie, Gauthier-Villars, Paris, 1931

[2] M.P. Hassell; G.C. Varley New inductive population model for insect parasites and its bearing on biological control, Nature, Volume 223 (1969), pp. 1133-1137

[3] C.A. Free; J.R. Beddington; J.H. Lawton On the inadequacy of simple models of mutual interference for parasitism and predation, J. Anim. Ecol., Volume 36 (1977), pp. 375-389

[4] C.S. Holling The components of predation as revealed by a study of small mammal predation of the European pine sawfly, Can. Entomol., Volume 91 (1959), pp. 293-320

[5] R. Arditi; H.R. Akçakaya Underestimation of mutual interference of predators, Oecologia, Volume 83 (1990), pp. 358-361

[6] J.R. Beddington Mutual interference between parasites or predators and its effect on searching efficiency, J. Anim. Ecol., Volume 45 (1975), pp. 331-340

[7] D.L. DeAngelis; R.A. Goldstein; R.V. O'Neill A model for trophic interaction, Ecology, Volume 56 (1975), pp. 881-892

[8] M.P. Hassell; D.J. Rogers Insect parasite responses in the development of population models, J. Anim. Ecol., Volume 41 (1972), pp. 661-676

[9] M.P. Hassell; R.M. May Stability in insect host-parasite models, J. Anim. Ecol., Volume 42 (1973), pp. 693-726

[10] D. Rogers; M.P. Hassell General models for insect parasite and predator searching behaviour: interference, J. Anim. Ecol., Volume 43 (1974), pp. 239-253

[11] T.W. Hwang Global analysis of the predator–prey system with Beddington–DeAngelis functional response, J. Math. Anal. Appl., Volume 281 (2003), pp. 395-401

[12] T.W. Hwang Uniqueness of limit cycles of the predator–prey system with Beddington–DeAngelis functional response, J. Math. Anal. Appl., Volume 290 (2004), pp. 113-122

[13] M. Fan; Y. Kuang Dynamics of a nonautonomous predator–prey system with the Beddington–DeAngelis functional response, J. Math. Anal. Appl., Volume 295 (2004), pp. 15-39

[14] R.S. Cantrell; C. Cosner On the dynamics of predator–prey models with the Beddington–DeAngelis functional response, J. Math. Anal. Appl., Volume 257 (2001), pp. 206-222

[15] R.S. Cantrell; C. Cosner Effects of domain size on the persistence of populations in a diffusive food chain model with DeAngelis–Beddington functional response, Nat. Resour. Model., Volume 14 (2001), pp. 335-367

[16] M. Begon; M. Mortimer; D.J. Thompson Population Ecology: A Unified Study of Animals and Plants, Blackwell Science, Oxford, UK, 1996

[17] M.G. Neubert; H. Caswell Alternatives to resilience for measuring the responses of ecological systems to perturbations, Ecology, Volume 78 (1997), pp. 653-665

[18] M.P. Hassell The Dynamics of Arthropod Predator–Prey Systems, Princeton University Press, Princeton, NJ, USA, 1978

[19] M.L. Rosenzweig Paradox of enrichment: destabilization of exploitation ecosystem in ecological time, Science, Volume 171 (1971), pp. 385-387

[20] R. Arditi; L.R. Ginzburg Coupling in predator–prey dynamics: ratio-dependence, J. Theor. Biol., Volume 139 (1989), pp. 311-326

[21] J.R. Beddington; C.A. Free; J.H. Lawton Characteristics of successful natural enemies in models of biological control of insect pest, Nature (Lond.), Volume 273 (1978), pp. 513-519

[22] C. Jost; O. Arino; R. Arditi About deterministic extinction in ratio-dependent predator–prey models, Bull. Math. Biol., Volume 61 (1999), pp. 19-32

[23] F. Berezovskaya; G. Karev; R. Arditi Parametric analysis of the ratio-dependent predator–prey model, J. Math. Biol., Volume 43 (2001), pp. 221-246

[24] G.F. Gause Experimental demonstrations of Volterra's periodic oscillations in the numbers of animals, J. Exp. Biol., Volume 12 (1935), pp. 44-48

[25] L.S. Luckinbill Coexistence in laboratory populations of Paramecium aurelia and its predator Didinium nasutum, Ecology, Volume 54 (1973), pp. 1320-1327

[26] M.L. Rosenzweig; R.H. MacArthur Graphical representation and stability conditions of predator–prey interactions, Am. Natural., Volume 97 (1963), pp. 217-223

[27] M.R. Myerscough; M.J. Darwen; W.L. Hogarth Stability, persistence and structural stability in a classical predator–prey model, Ecol. Modelling, Volume 89 (1996), pp. 31-42

[28] Y. Kuang; E. Beretta Global qualitative analysis of a ratio-dependent predator–prey system, J. Math. Biol., Volume 36 (1998), pp. 389-406

[29] Y. Kuang Rich dynamics of Gause-type ratio-dependent predator–prey system, Fields Inst. Commun., Volume 21 (1999), pp. 325-337

[30] R.M. May Stability and Complexity in Model Ecosystems, Princeton University Press, Princeton, NJ, USA, 1973

[31] D. Tilman Secondary succession and the pattern of plant dominance along experimental nutrient gradients, Ecol. Monogr., Volume 57 (1987), pp. 189-214

[32] R.L. Haedrich; G.T. Rowe; P.T. Poloni The megabenthic fauna in the deep-sea south of New England, USA, Mar. Biol., Volume 57 (1980), pp. 165-179

[33] J.F. Owen On productivity as predictor of rodent and carnivore diversity, Ecology, Volume 69 (1988), pp. 1161-1165

[34] M.L. Rosenzweig Species Diversity in Space and Time, Cambridge University Press, Cambridge, 1993

[35] J.F. Riebesell Paradox of enrichment in competitive systems, Ecology, Volume 55 (1974), pp. 183-187

[36] K.J. Griffiths; C.S. Holling A competition submodel for parasites and predators, Can. Entomol., Volume 101 (1969), pp. 785-818

[37] M. Munster-Swendsen; G. Nachman Asincrony in insect host-parasite interaction and its effect on stability, studied by a population model, J. Anim. Ecol., Volume 47 (1978), pp. 159-171

[38] S. Ponsard; R. Arditi; C. Jost Assessing top-down and bottom-up control in a litter-based soil macroinvertebrate food chain, Oikos, Volume 89 (2000), pp. 524-540

[39] C. Jost; R. Arditi From pattern to process: identifying predator–prey interactions, Popul. Ecol., Volume 43 (2001), pp. 229-243

[40] W.J. Sutherland From Individual Behaviour to Population Ecology, Oxford University Press, Oxford, UK, 1996

[41] J. Monod Recherches sur la croissance des cultures bactériennes, Hermann, Paris, 1942

[42] C. Jost Predator–prey theory: Hidden twins in ecology and microbiology, Oikos, Volume 90 (2000), pp. 202-208

[43] M.P. Hassel The Spatial and Temporal Dynamics of Host-Parasitoid Interactions, Oxford University Press, Oxford, UK, 2000

Comments - Policy