The Rosenzweig–MacArthur Graphical Criterion for a Predator-Prey Model with Variable Mortality Rate

Abstract

We consider a general modified Gause type model of predation, for which the predator mortality rate can depend on the densities of both species, prey and predator. We give a graphical criterion for the stability of positive hyperbolic equilibria, which is an extension of the well-known Rosenzweig–MacArthur graphical criterion for the case of a constant predator mortality rate. We examine the occurrence of a Poincaré–Andronov–Hopf bifurcation and give an expression for the first Lyapunov coefficient. Our model generalizes several models appearing in the literature. The relevance of our results, i.e. the use of the graphical criterion and the expression for the first Lyapunov coefficient, is tested on these models. The global behavior of the system is illustrated by numerical simulations which confirm the local properties of the models near the equilibria.

Appendices

A Proof of Theorem 4

First, by hypotheses (24) and (25) the necessary conditions (20) for a PAH bifurcation are satisfied. Secondly, we prove the transversality condition

$$\begin{aligned} \dfrac{d}{d x^*} {\text {tr}} \mathcal J(x^*,h(x^*))\vert _{x^*={\tilde{x}}}\ne 0. \end{aligned}$$

From the expression of the trace given in the proof of Theorem 2, we have

$$\begin{aligned} {\text {tr}} {\mathcal {J}}(x^*, h(x^*))=h(x^*) \left( H(x^*)-G(x^*)\right) . \end{aligned}$$

Therefore,

$$\begin{aligned} \dfrac{d}{d x^*} {\text {tr}} {\mathcal {J}}(x^*,h(x^*))= h'(x^*) \left( H(x^*)-G(x^*)\right) +h(x^*) \left( H'(x^*)-G'(x^*)\right) . \end{aligned}$$

By hypothesis (25), we obtain

$$\begin{aligned} \dfrac{d}{d x^*} {\text {tr}} {\mathcal {J}}(x^*, y^*)\vert _{x^*={\tilde{x}}}&=h({\tilde{x}})\left[ H'({\tilde{x}})-G'({\tilde{x}})\right] . \end{aligned}$$

From the assumption (26), the transversality condition is verified. Thirdly, in order to examine the non-degeneracy condition and to compute the first Lyapunov coefficient $\rho$, we introduce the following change of variables

$$\begin{aligned} x=N,\quad y=h'({\tilde{x}})N+ \frac{\omega }{p({\tilde{x}})}P, \end{aligned}$$

(46)

where $\omega =\sqrt{\det {\mathcal {J}}({\tilde{x}},h({\tilde{x}}))}$. The model (2) becomes

$$\begin{aligned} \left\{ \begin{array}{lll} {\dot{N}}=-\omega P+F(N, P),\\ {\dot{P}}=\omega N+G(N, P), \end{array} \right. \end{aligned}$$

(47)

with

$$\begin{aligned} F(N, P)=p(N)\left[ h(N)-h'({\tilde{x}})N \right] -\omega P\left[ \dfrac{p(N)}{p({\tilde{x}})}-1\right] , \end{aligned}$$

and

$$\begin{aligned} \begin{aligned} G(N, P)&=\dfrac{1}{\omega }\Bigl [q(N)-D\bigl (N, h'({\tilde{x}})N+ \dfrac{\omega }{p({\tilde{x}})}P \bigr )-m\Bigr ][p({\tilde{x}})h'({\tilde{x}})N+\omega P]\\&\quad -\dfrac{h'({\tilde{x}})p(N)}{\omega }\Bigl [p({\tilde{x}})h(N)-p({\tilde{x}})h'({\tilde{x}})N-\omega P \Bigr ]-\omega N. \end{aligned} \end{aligned}$$

The parameter $\rho$ is given by

$$\begin{aligned} \rho =\dfrac{1}{16}A_1+\dfrac{1}{16\omega }\left[ A_2-A_3-A_4\right] , \end{aligned}$$

(48)

where (see [14], page 169 or [10], page 152)

$$\begin{aligned} A_1=F_{NNN}+F_{NPP}+G_{NNP}+G_{PPP},\quad A_2=F_{NP}(F_{NN}+F_{PP}), \end{aligned}$$

$$\begin{aligned} A_3=G_{NP}(G_{NN}+G_{PP}),\quad A_4=F_{NN}G_{NN}-F_{PP}G_{PP}. \end{aligned}$$

where $F_{NN}$ denotes $\frac{\partial ^2F}{\partial N \partial N}({\tilde{N}},{\tilde{P}})$, $F_{NP}$ denotes $\frac{\partial ^2F}{\partial N \partial P}({\tilde{N}},{\tilde{P}})$, and similarly for all other partial derivatives. Since F is linear in P, we have $F_{PP}=F_{NPP}=0$. Notice that from the change of variables (46), we have ${\tilde{N}}={\tilde{x}}$ and ${\tilde{P}}=\frac{p({\tilde{x}})}{\omega }\left( {\tilde{y}}-h'\left( {\tilde{x}}\right) {\tilde{x}}\right)$.

After calculations, and using the notations (21), we find that $A_i$, $i=1,\cdots ,4$ are given by the following expressions

$$\begin{aligned} \begin{aligned} A_1&=\dfrac{1}{p_0^2}\Bigl [ p_0^2p_2h_1+3p_0^2p_1h_2+p_0^3h_3 +p_0^2q_2 -p_0^2d_{11} -4p_0^2h_1d_{12} -3\left( \omega ^2+p_0^2h_1^2\right) d_{22}\\&\qquad -p_0^2h_0d_{112}-\left( \omega ^2+p_0^2h_1^2\right) h_0d_{222}-2p_0^2h_0h_1d_{122}\Bigr ],\\ A_2&=-\omega p_1h_2,\\ A_3&=\dfrac{1}{p_0\omega } \Bigl [ p_1h_1+q_1-d_1-2h_1d_2-h_0d_{12}-h_0h_1d_{22} \Bigr ] \Bigl [ p_0^2q_2h_0-p_0^3h_1h_2+2p_0^2q_1h_1\\&\qquad -2p_0^2h_1d_1-2\left( \omega ^2+p_0^2h_1^2\right) d_2 -p_0^2h_0d_{11}-2p_0^2h_0h_1d_{12} -\left( \omega ^2+p_0^2h_1^2\right) h_0d_{22} \Bigr ],\\ A_4&=\dfrac{p_0^2h_2}{\omega } \Bigl [ q_2h_0-p_0h_1h_2+2q_1h_1-2h_1d_1-2h_1^2d_2 -h_0d_{11}-2h_0h_1d_{12}-h_0h_1^2d_{22} \Bigr ]. \end{aligned} \end{aligned}$$

Replacing these expressions in (48) we obtain the following formula for $\rho$

$$\begin{aligned} \rho =\frac{1}{16p_0^2\omega ^2}\left( a_0+a_2\omega ^2+a_4\omega ^4\right) , \end{aligned}$$

where $a_0$, $a_2$ and $a_4$ do not depend on $\omega$. Now, using MAPLE [23], we can replace $\omega ^2$ by its expression (22) and $h_1$ by $h_1=d_2h_0/p_0$, which follows from the condition $H({\tilde{x}})=G({\tilde{x}})$. We obtain the expression of $\rho$ given by (23).

Under the assumption that $\rho \ne 0$ we conclude that (2) undergoes a non degenerate PAH bifurcation at ${\tilde{E}}({\tilde{x}}, {\tilde{y}})$, see Theorem 3.4.2 in [10].

B Some Biological Explanations

It is not useless to give briefly some biological explanations for the models we used to illustrate our results (see for example Table 1 in [38]). In RMA and Gause models, the predator density is regulated only by food limitation ($d(x,y)=m$). In the Hsu model, the mortality rate is a decreasing function of the prey density. In [13], Hsu indicates for example the case $d(x)=(ex+f)/(rx+s)$ with $fe>se$ which can be written in the form $m+\alpha /(\delta +x)$. In this reference, there is no justification for this rate. However, it can be interpreted naturally by the fact that the more prey is available, the more the predator mortality rate is small. Maybe a better justification is to consider d(x)y as an emigration term, in the sense that predators get out of the ecosystem when the density of the prey is small. The Hsu model contains also a model due to Minter et al. [24] where a detailed biological justification leads to a mortality rate of the form $d(x)=\alpha /(\delta +x)$ for which d(x) decreases to zero when x goes to infinity. Bazykin [2, 3, 21] introduced the regulation by intraspecific mechanisms, that is a competition among predators for resources other than prey. To do this, he subtracted a quantity $\alpha y^2$ from the predator equation ($d(x,y)=m+\alpha y$). Even if not concerned by our study, it is also worth mentioning the recent article [16] which deals with Bazykin model with ratio-dependent functional response of Arditi-Ginzburg. If the predator density-dependence parameter $\alpha$ is made inversely proportional to resource availability, the dynamics are described by the variable-territory model of Turchin-Batzli [37, 38] ($d(x,y)=m+\alpha y/x$, where $\alpha$ is called the prey/predator ratio at equilibrium). The term $d(x,y)y=\alpha y^2/x$ in the predator equation represents the self-limitation of the predator. Nevertheless, a self-limitation should be biologically limited. It is not the case for Turchin-Batzli model where $\alpha y^{2}/x$ could become quite large in case of predators can subsist on few prey (i.e. $\alpha y^2/x\rightarrow \infty$ as $x\rightarrow 0$). This leads to the well-defined modified VT model for which $d(x,y)=m+\frac{\alpha y }{\delta +x}$ where $\delta >0$ is fixed [15]. Not only the self limitation of the predator becomes less than $\alpha y^{2}/d$ but also the singularity in $x=0$ is avoided. Cavani and Farkas [6, 8] proposed the modified RMA model for which the mortality of the predator in the absence of the prey is a growing and bounded function of the predator quantity. More exactly $d(x,y)=m+\alpha y/(1+y)$ where, contrary to what we assumed for technical reasons, $\alpha =\delta -m$ depends on m. Here, $m>0$ is the mortality at low density and is less than $\delta$ which is the limiting, maximal mortality. Except for Gause, Hsu and Bazykin models, the theoretical study of these models is not, to our knowledge, always exhaustive. Authors have sometimes directed their investigation to a particular aspect, especially the appearance of periodic orbits. Indeed, all the mentioned models can generically exhibit limit cycles that can appear or disappear by PAH bifurcation.