## 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:

$$\mathrm{d}N/\mathrm{d}t=F(N)-g(N\text{,}P)P\mathrm{d}P/\mathrm{d}t=eg(N\text{,}P)P-qP$$ | (1) |

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\text{,}P)$. The well-known models of Lotka–Volterra

$$(\mathrm{LV})\phantom{\rule{2em}{0ex}}g(N\text{,}P)=aN$$ |

$$(\mathrm{H})\phantom{\rule{2em}{0ex}}g(N\text{,}P)=\frac{aN}{1+ahN}$$ |

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:

$$(\mathrm{HV})\phantom{\rule{2em}{0ex}}a=\alpha {P}^{-m}$$ | (2) |

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):

$$(\mathrm{HVLV})\phantom{\rule{2em}{0ex}}g(N\text{,}P)=\alpha 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]:

$$(\mathrm{HVH})\phantom{\rule{2em}{0ex}}g(N\text{,}P)=\frac{\alpha N{P}^{-m}}{1+\alpha hN{P}^{-m}}$$ | (3) |

$$(\mathrm{RD})\phantom{\rule{2em}{0ex}}g(N\text{,}P)=\frac{\alpha N/P}{1+\alpha hN/P}$$ | (4) |

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\text{,}P)=\frac{aN}{1+ahN+aw(P-{P}_{0})}$$ | (5) |

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:

$$(\mathrm{BDA})\phantom{\rule{2em}{0ex}}g(N\text{,}P)=\frac{aN}{1+ahN+awP}$$ | (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\text{,}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:

$$\mathrm{d}P/\mathrm{d}t=e\frac{\alpha N{P}^{-m}}{1+\alpha hN{P}^{-m}}P-qP=0$$ |

$${P}^{m}=\frac{\alpha (e-hq)}{q}N$$ |

$$e-hq>0$$ | (7) |

Fig. 1a–e shows the shape of the predator isocline $({\mathrm{Iso}}_{\mathrm{P}})$ in the phase plane $(N\text{,}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).

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

$$\mathrm{d}P/\mathrm{d}t=e\frac{aN}{1+ahN+awP}P-qP=0$$ |

$$P=\frac{e-hq}{wq}N-\frac{1}{aw}$$ |

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}_{\mathrm{C}}=q/(a(e-hq))$, 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\mathbf{J}<0\text{;}\phantom{\rule{2em}{0ex}}\mathrm{det}\mathbf{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\text{,}0)$ and $(K\text{,}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\u2a7dm<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:

$$\mathrm{d}N/\mathrm{d}t=rN-\frac{\alpha N{P}^{-m}}{1+\alpha hN{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-hq>0$, there exists a unique non-trivial equilibrium, except when $m=1$. It is given by:

$${N}^{*}=\frac{q}{er}{\left(\frac{er}{\alpha (e-hq)}\right)}^{1/(1-m)}$$ |

$${P}^{*}={\left(\frac{er}{\alpha (e-hq)}\right)}^{1/(1-m)}$$ |

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

$$0<\frac{rh}{e-hq}<m<1$$ | (8) |

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 $er>\alpha (e-hq)$, the equilibrium densities rise with m, when $0<m<1$, and decrease when $m\u2a7e1$. If $er<\alpha (e-hq)$, it is the opposite.

#### 3.1.2 BDA model

For the BDA model, the prey equation is:

$$\mathrm{d}N/\mathrm{d}t=rN-\frac{aN}{1+ahN+awP}P$$ |

$$w<\frac{e-hq}{er}$$ | (9) |

If the equilibrium exists, it is given by

$${N}^{*}=\frac{q}{a(e-hq-erw)}$$ |

$${P}^{*}=\frac{er}{a(e-hq-erw)}$$ |

$$\frac{h}{e}<w<\frac{e-hq}{er}$$ | (10) |

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-hq)/er$, 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{\mathrm{d}N}{\mathrm{d}t}=rN(1-N/K)-\frac{\alpha N{P}^{-m}}{1+\alpha hN{P}^{-m}}P$$ |

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

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-hq)}<K$$ | (11) |

$$\frac{q}{a(e-hq)}<K<\frac{(e+hq)}{ah(e-hq)}$$ | (12) |

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\text{,}0)$ and $(K\text{,}0)$, there is a non-trivial equilibrium if

$$\alpha (e-hq)<re$$ |

$${N}^{*}=K(1-\frac{\alpha (e-hq)}{re})$$ |

$${P}^{*}={N}^{*}\frac{\alpha (e-hq)}{q}$$ |

$$r>\frac{\alpha (e-hq)}{e}(1+q(\frac{h}{e}-\frac{1}{\alpha}))$$ |

If the equilibrium is stable, stability is not always global [22,23]. According to initial conditions, there is either convergence to $(0\text{,}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\text{,}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\text{,}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\u2a7e1$, 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={\left(a{P}_{0}\right)}^{\mathrm{-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\u2a7e1$. 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{\mathrm{d}N}{\mathrm{d}t}=rN(1-N/K)-\frac{aN}{1+ahN+awP}P$$ |

$$\frac{q}{a(e-hq)}<K$$ |

The non-trivial equilibrium is given by:

$${N}^{*}=\frac{K}{2er}(er-\frac{e-hq}{w}+\sqrt{{(er-\frac{e-hq}{w})}^{2}+\frac{4erq}{Kwa}}\phantom{\rule{0.2em}{0ex}})$$ |

$${P}^{*}=\frac{1}{waq}(a(e-hq){N}^{*}-q)$$ |

$${N}_{\infty}^{*}=\frac{q}{a(wer-e+hq)}$$ |

$${N}^{*}\approx (1-\frac{e-hq}{erw})K$$ | (13) |

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$$ |

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-hq>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\text{,}m)$ for the HVH model (Fig. 3b and d), and in the plane $(r\text{,}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\text{,}m)$ for the HVH model, and in the plane $(K\text{,}w)$ for the BDA model.

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}_{\mathrm{C}}<m<1$ in which there is always stability (with ${m}_{\mathrm{C}}=rh/(e-hq)$, 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 $\mathrm{log}N$–$\mathrm{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:

$$\mathrm{d}N/\mathrm{d}t={f}_{N}(N\text{,}P)=F(N)-g(N{P}^{-m})P$$ |

$$\mathrm{d}P/\mathrm{d}t={f}_{P}(N\text{,}P)=eg(N{P}^{-m})P-qP$$ |

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 $\mathrm{d}P/\mathrm{d}t=0$.

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

$$g(\theta )=\frac{\alpha \theta}{1+\alpha h\theta}\phantom{\rule{1em}{0ex}}(\text{with}\phantom{\rule{0.25em}{0ex}}{g}^{\prime}(\theta )=\frac{\alpha}{{(1+\alpha h\theta )}^{2}})$$ |

$${g}^{\prime}(N{P}^{-m})=\frac{\alpha}{{(1+h\alpha N{P}^{-m})}^{2}}=\frac{1}{\alpha {\left(N{P}^{-m}\right)}^{2}}{(g(N{P}^{-m}))}^{2}=\frac{(\alpha {(e-hq))}^{2}}{\alpha {e}^{2}}$$ |

$$\mathbf{J}=\left[\begin{array}{cc}\hfill {F}^{\prime}(N)-{P}^{1-m}\frac{(\alpha {(e-hq))}^{2}}{\alpha {e}^{2}}\hfill & \hfill -\frac{q}{e}[1-m+m\frac{hq}{e}]\hfill \\ \hfill {P}^{1-m}\frac{(\alpha {(e-hq))}^{2}}{\alpha e}\hfill & \hfill -m\frac{q}{e}(e-hq)\hfill \end{array}\right]$$ | (A.1) |

$$\left[\begin{array}{cc}\hfill *\hfill & \hfill <0\hfill \\ \hfill >0\hfill & \hfill <0\hfill \end{array}\right]$$ |

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\mathbf{J}<0$ and $\mathrm{det}\mathbf{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}^{\prime}(0)<\infty $. To see this, all we have to do is to find the sign of:

$$\partial {f}_{N}/\partial P=mN{P}^{-m}{g}^{\prime}(N{P}^{-m})-g(N{P}^{-m})$$ |

The derivative of the function $x\mapsto mx{g}^{\prime}(x)-g(x)$ is $mx{g}^{\u2033}(x)-(1-m){g}^{\prime}(x)$. If we assume that g is concave (i.e. ${g}^{\u2033}<0$) and increasing (i.e. ${g}^{\prime}>0$) and $0<m<1$, then this latter expression is negative. The function $x\mapsto mx{g}^{\prime}(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{aN}{1+awP+ahN}=\frac{q}{e}\phantom{\rule{1em}{0ex}}\text{and}$$ |

$$\frac{aP}{1+awP+ahN}=F(N)/N$$ |

$$\mathbf{J}=\left[\begin{array}{cc}\hfill {F}^{\prime}(N)-\frac{e-hq}{e}\frac{F(N)}{N}\hfill & \hfill \frac{q}{e}(w\frac{F(N)}{N}-1)\hfill \\ \hfill (e-hq)\frac{F(N)}{N}\hfill & \hfill -qw\frac{F(N)}{N}\hfill \end{array}\right]$$ | (A.2) |

Note that we have $\partial {f}_{N}/\partial P<0$ because $wF(N)/N=awP/(1+ahN+awP)<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)=rN$.

### B.1 HVH model

The prey isocline is given by

$$N=\frac{1}{\alpha hr}(\alpha P-r{P}^{m})$$ |

$$N=\frac{q}{\alpha (e-hq)}{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-hq>0$. Three cases have to be distinguished with respect to the value of m (presented graphically in the $(N\text{,}P)$ phase plane, see Fig. 1):

- – if $0\u2a7dm<1$, the prey isocline implicitly defines a monotonous concave function of $N\in [0\text{,}+\infty [$ that increases from ${(r/\alpha )}^{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>hq$ and $\alpha >r$ are required for them to pass through the positive quadrant; two cases may be distinguished according to the slopes of the two isoclines, $er>\alpha (e-hq)$ and $er<\alpha (e-hq)$; 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/\alpha}$, 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\ne 1$. Analytically, we obtain:

$${N}^{*}=\frac{q}{er}{\left(\frac{er}{\alpha (e-hq)}\right)}^{1/(1-m)}$$ |

$${P}^{*}={\left(\frac{er}{\alpha (e-hq)}\right)}^{1/(1-m)}$$ |

Note that the ratio ${N}^{*}/{P}^{*}$ is equal to the value $q/(er)$ 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\text{,}0)$, as in the case $m=1$.

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

$$\mathbf{J}=\left[\begin{array}{cc}\hfill r\frac{hq}{e}\hfill & \hfill -\frac{q}{e}((1-m)+m\frac{hq}{e})\hfill \\ \hfill r(e-hq)\hfill & \hfill -m\frac{q}{e}(e-hq)\hfill \end{array}\right]$$ |

(B.1) |

$$m>\frac{rh}{e-hq}\iff h<\frac{me}{r+mq}$$ | (B.2) |

### B.2 BDA model

The non-trivial equilibrium is

$${N}^{*}=\frac{q}{a(e-hq-erw)}\text{,}\phantom{\rule{2em}{0ex}}{P}^{*}=\frac{er}{q}{N}^{*}$$ |

$$P=\frac{r(1+ahN)}{a(1-wr)}$$ |

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

$$\mathbf{J}=\left[\begin{array}{cc}\hfill r\frac{hq}{e}\hfill & \hfill \frac{q}{e}(wr-1)\hfill \\ \hfill r(e-hq)\hfill & \hfill -qwr\hfill \end{array}\right]$$ |

$$\frac{h}{e}<w<\frac{e-hq}{er}$$ |

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

Let $F(N)=rN(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\alpha h}{K}{N}^{2}+r(\frac{1}{K}{P}^{m}-\alpha h)N+\alpha P-r{P}^{m}=0$$ |

$$N={\phi}_{\pm}(P)=\frac{K}{2r\alpha h}[r(\alpha h-\frac{1}{K}{P}^{m})\pm \sqrt{{r}^{2}{(\alpha h-\frac{1}{K}{P}^{m})}^{2}-\frac{4r\alpha h}{K}(\alpha P-r{P}^{m})}\phantom{\rule{0.2em}{0ex}}]$$ | (C.1) |

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

- – ${\left(\alpha hK\right)}^{1-m}<{(r/\alpha )}^{m}$; ${\phi}_{-}$ is always negative and the isocline is given by ${\phi}_{+}$ only, which is positive between $P=0$ and $P={(r/\alpha )}^{1/(1-m)}$ where it decreases from $N=K$ to $N=0$ (see Fig. 2b);
- – ${\left(\alpha hK\right)}^{1-m}>{(r/\alpha )}^{m}$; a graphical argument shows that the expression under the square root becomes zero at some value $\overline{P}$ where ${\phi}_{-}$ and ${\phi}_{+}$ 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 $\alpha >r$ then the isocline is a parabola-like curve that passes through the trivial equilibria $(0\text{,}0)$ and $(K\text{,}0)$. In that case, a non-trivial equilibrium exists if $\alpha (e-hq)<re$;
- – if $\alpha <r$ then it is a decreasing curve defined on $[K\frac{r-\alpha}{r}\text{,}K]$ that passes through $(K\text{,}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={\phi}_{-}(P)$ and $N={\phi}_{+}(P)$. The following remarks can be made:

- – ${\phi}_{-}(0)=0$ and ${\phi}_{+}(0)=K$;
- – ${\phi}_{-}$ increases for small P; this is because the term in P under the square root dominates the terms in ${P}^{m}$;
- – ${\phi}_{-}$ only crosses the axe $N=0$ in two points: $P=0$ and $P={(\alpha /r)}^{1/(m-1)}$; beyond this domain, ${\phi}_{-}<0$ if it exists;
- – ${\phi}_{+}<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>{(\alpha /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 ${\left(\alpha hK\right)}^{1/m}<{(\alpha /r)}^{1/(m-1)}$, i.e., for sufficiently small K values, the expression becomes negative in a bounded interval of $]0\text{,}+\infty ]$ containing ${\left(\alpha hK\right)}^{1/m}$. Again, at the endpoints of this interval, ${\phi}_{-}={\phi}_{+}$.

Now one can have an idea of the possible shapes of the isoclines in the $(N\text{,}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{\alpha e-h\alpha q}{q}N\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}P=\frac{r}{q}N(1-\frac{N}{K})$$ |

$${N}^{*2}-K{N}^{*}+K\frac{q}{er}{\left(\frac{\alpha e-\alpha hq}{q}{N}^{*}\right)}^{1/m}=0$$ | (C.2) |

This equation is of the form:

$$\psi (N\text{,}K\text{,}\gamma )={N}^{2}-KN+K\gamma {\beta}^{1/m}{N}^{1/m}=0$$ | (C.3) |

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

$$\mathbf{J}=\left[\begin{array}{cc}\hfill r(1-2\frac{{N}^{*}}{K})-{P}^{*1-m}\frac{{(\alpha e-\alpha hq)}^{2}}{\alpha {e}^{2}}\hfill & \hfill -\frac{q}{e}(1-m+m\frac{hq}{e})\hfill \\ \hfill e{P}^{*1-m}\frac{{(\alpha e-\alpha hq)}^{2}}{\alpha {e}^{2}}\hfill & \hfill -\frac{q}{e}m(e-hq)\hfill \end{array}\right]$$ |

$$r(1-2{N}^{*}/K)=2r/{N}^{*}({N}^{*}-{N}^{*2}/K)-r=2\frac{q}{e}{\left(\frac{\alpha e-\alpha hq}{q}\right)}^{1/m}{N}^{*(1-m)/m}-r$$ |

$$\mathbf{J}=\left[\begin{array}{cc}\hfill \frac{q}{{e}^{2}}{\left(\frac{\alpha (e-hq)}{q}\right)}^{\frac{1}{m}}(e+hq){N}^{*\frac{1-m}{m}}-r\hfill & \hfill -\frac{q}{e}(1-m+m\frac{hq}{e})\hfill \\ \hfill \frac{q}{e}{\left(\frac{\alpha (e-hq)}{q}\right)}^{\frac{1}{m}}(e-hq){N}^{*\frac{1-m}{m}}\hfill & \hfill -\frac{q}{e}m(e-hq)\hfill \end{array}\right]$$ | (C.4) |

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

$$\frac{{(N/K)}^{1-m}}{{(1-N/K)}^{m}}=\frac{1}{{\gamma}^{m}\beta {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:

$$\theta \in \phantom{\rule{0.2em}{0ex}}]0\text{,}1[\phantom{\rule{0.2em}{0ex}}\mapsto \frac{{\theta}^{1-m}}{{(1-\theta )}^{m}}$$ |

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

$$\frac{\mathrm{d}N}{\mathrm{d}K}=-\frac{\partial \psi /\partial K}{\partial \psi /\partial N}=\frac{m{N}^{2}}{K(mN+(1-m)(K-N))}>0$$ |

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

$$\mathbf{J}\to \left[\begin{array}{cc}\hfill -r\hfill & \hfill -\frac{q}{e}(1-m+m\frac{hq}{e})\hfill \\ \hfill 0\hfill & \hfill -\frac{q}{e}m(e-hq)\hfill \end{array}\right]$$ |

**J**increasing, and

$$\mathrm{det}\mathbf{J}=\frac{q}{{e}^{2}}(e-hq)({\left(\frac{\alpha (e-hq)}{q}\right)}^{1/m}(1-2m){N}^{*(1-m)/m}q+emr)$$ |

**J**. 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 ${\left(\alpha hK\right)}^{1-m}<{(r/\alpha )}^{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{\alpha (e-hq)}{re})$$ |

$${P}^{*}={N}^{*}\frac{\alpha (e-hq)}{q}$$ |

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

$$\mathbf{J}=\left[\begin{array}{cc}\hfill r(1-2{N}^{*}/K)-\frac{(\alpha {(e-hq))}^{2}}{\alpha {e}^{2}}\hfill & \hfill -\frac{h{q}^{2}}{{e}^{2}}\hfill \\ \hfill \frac{(\alpha {(e-hq))}^{2}}{\alpha e}\hfill & \hfill -\frac{q}{e}(e-hq)\hfill \end{array}\right]$$ |

$$\mathbf{J}=\left[\begin{array}{cc}\hfill \frac{\alpha ({e}^{2}-{h}^{2}{q}^{2})}{{e}^{2}}-r\hfill & \hfill -h\frac{{q}^{2}}{{e}^{2}}\hfill \\ \hfill \frac{\alpha {(e-hq)}^{2}}{e}\hfill & \hfill -\frac{q}{e}(e-hq)\hfill \end{array}\right]$$ |

Calculating Tr**J** and $\mathrm{det}\mathbf{J}$, gives:

$$Tr\mathbf{J}=-r+\frac{\alpha ({e}^{2}-{h}^{2}{q}^{2})}{{e}^{2}}-\frac{q}{e}(e-hq)=-r+\frac{\alpha (e-hq)}{e}(1+q(\frac{h}{e}-\frac{1}{\alpha}))$$ |

$$\mathrm{det}\mathbf{J}=\frac{q}{e}(e-hq)\left|\begin{array}{cc}\hfill \frac{\alpha ({e}^{2}-{h}^{2}{q}^{2})}{{e}^{2}}-r\hfill & \hfill -h\frac{{q}^{2}}{{e}^{2}}\hfill \\ \hfill \frac{\alpha (e-hq)}{e}\hfill & \hfill -1\hfill \end{array}\right|=(r-\frac{\alpha}{e}(e-hq))q\frac{(e-hq)}{e}$$ |

Since the criterion for the positive equilibrium to exist is $re>a(e-hq)$, $\mathrm{det}\mathbf{J}>0$, and the criterion for stability reads:

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

In particular, one can see that $e>ha$ 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{{\gamma}^{m}\beta}{{K}^{m-1}}$$ | (C.5) |

Now, the function $\theta \in \phantom{\rule{0.2em}{0ex}}]0\text{,}1[\phantom{\rule{0.2em}{0ex}}\to {\theta}^{m-1}{(1-\theta )}^{m}$ increases from 0 to a maximum in $]0\text{,}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 ${\left({\gamma}^{m}\beta \right)}^{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\text{,}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+ahN)}{a(1-rw)+arwN/K}$$ |

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

- – $rw<1$: we get a parabola like curve that reaches its maximum between 0 and K and that becomes negative for $N>K$.
- – $rw>1$: we get a decreasing curve that is positive in $[K\frac{rw-1}{rw}\text{,}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-hq)}<K$$ | (C.6) |

$$r{N}^{*}(1-{N}^{*}/K)=\frac{q}{e}{P}^{*}$$ | (C.7) |

$${P}^{*}=\frac{1}{waq}(a(e-hq){N}^{*}-q)$$ |

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

$$\frac{er}{qK}{N}^{*2}+(\frac{e-hq}{wq}-\frac{er}{q}){N}^{*}-\frac{1}{wa}=0$$ | (C.8) |

$${N}^{*}=\frac{K}{2er}(er-\frac{e-hq}{w}+\sqrt{{(er-\frac{e-hq}{w})}^{2}+\frac{4erq}{Kwa}}\phantom{\rule{0.2em}{0ex}})$$ | (C.9) |

There are two possibilities:

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

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

$$\mathbf{J}=\left[\begin{array}{cc}\hfill r(-N/K+\frac{hq}{e}(1-N/K))\hfill & \hfill \frac{q}{e}(wr(1-N/K)-1)\hfill \\ \hfill r(e-hq)(1-N/K)\hfill & \hfill -qwr(1-N/K)\hfill \end{array}\right]$$ |

Thus, the stability conditions are:

(C.10) |

## 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:

$$\left[\begin{array}{cc}\hfill \gamma -r\hfill & \hfill \beta \hfill \\ \hfill \sigma \hfill & \hfill \delta \hfill \end{array}\right]$$ |

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 $er<(e-hq)/w$, then ${N}^{*}/K$ decreases from 1 to 0 when K increases from $q/(a(e-hq))$ to +∞. Thus, Tr**J** varies monotonously from −r to $rq(h/e-w)$. Now we have two cases:

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

In the second case, when $er>(e-hq)/w$, Tr**J** varies monotonously from −r to $r[q(\frac{h}{e}-w)\frac{e-hq}{erw}-\frac{erw-e+hq}{erw}]$. According to the parameter values, there is either always stability (and $h<we$ 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-hq))$, which is independent of w. Thus, when K is close to $q/(a(e-hq))$, condition (C.10) is always fulfilled and when K is large, it is only fulfilled when $w>h/e$.