A modified May-Holling-Tanner predator-prey model with multiple Allee effects on the prey and an alternative food source for the predator

TL;DR

This paper analyzes a predator-prey model incorporating multiple Allee effects, an alternative food source, and bifurcation phenomena, revealing complex dynamics and stability conditions through analytical and numerical methods.

Contribution

It introduces a modified predator-prey model with multiple Allee effects and an alternative food source, exploring its bifurcation structure and stability properties.

Findings

Existence of at most two equilibrium points in the first quadrant.

Identification of various bifurcations including saddle-node, Hopf, Bogadonov-Takens, and homoclinic bifurcations.

The basin of attraction for the stable equilibrium increases when depensation effects are reduced.

Abstract

We study a predator-prey model with Holling type I functional response, an alternative food source for the predator, and multiple Allee effects on the prey. We show that the model has at most two equilibrium points in the first quadrant, one is always a saddle point while the other can be a repeller or an attractor. Moreover, there is always a stable equilibrium point that corresponds to the persistence of the predator population and the extinction of the prey population. Additionally, we show that when the parameters are varied the model displays a wide range of different bifurcations, such as saddle-node bifurcations, Hopf bifurcations, Bogadonov-Takens bifurcations and homoclinic bifurcations. We use numerical simulations to illustrate the impact changing the predation rate, or the non-fertile prey population, and the proportion of alternative food source have on the basins of attraction of the stable equilibrium point in the first quadrant (when it exists). In particular, we also show that the basin of attraction of the stable positive equilibrium point in the first quadrant is bigger when we reduce the depensation in the model.

Equations

$$\frac{dx}{dt}$$

$$\frac{dy}{dt}$$

$$\tilde{K}(x)=nx\text{is replaced by}\overline{K}(x)=nx+c,$$

$$L_o(x)=rx\left(1-\frac{x}{K}\right)\text{is replaced by}{L}_m(x)=rx\left(1-\frac{x}{K}\right)\left(x-m\right).$$

$$L_b(x)=rx\left(1-\frac{x}{K}\right)\left(\frac{1}{x+b}\right)\left(x-m\right).$$

$$x_1=-b+\sqrt{(b+K)(b+m)},$$

$$x_2=\frac{1}{2}(K+m),$$

$$\frac{dx}{dt}$$

$$\frac{dy}{dt}$$

$$\frac{du}{d\tau }$$

$$\frac{dv}{d\tau }$$

$$\begin{array}{rl}& p(u):=(u-M)(1-u)=Q(u+C)(u+B)=:Qd(u),\text{and}{v}^{\ast }=u^{\ast }+C.\end{array}$$

$$\begin{array}{r}{u}_{1,2}=\frac{1}{2(1+Q)}\left(1+M-Q(B+C)\pm \sqrt{\Delta }\right)\text{with}\\ \Delta =(1+M-Q(B+C){)}^2-4(M+BCQ)(1+Q),\end{array}$$

$$\begin{array}{r}J(u,v)=\left(\begin{array}{cc}{J}_{11}(u,v)+J_{12}(u,v)& -uQh(u)\\ Sv(B+C+2u-v)& S(u+B)(C+u-2v)\end{array}\right),\end{array}$$

$$J(0,0)=\left(\begin{array}{cc}-CM& 0\\ 0& BCS\end{array}\right),$$

$${\psi }_{(0,0)}^1={\left(\begin{array}{cc}1& 0\end{array}\right)}^T\text{and}{\psi }_{(0,0)}^2={\left(\begin{array}{cc}0& 1\end{array}\right)}^T.$$

$$J(1,0)=\left(\begin{array}{cc}(M-1)(C+1)& -Q(B+1)(C+1)\\ 0& S(B+1)(C+1)\end{array}\right),$$

$${\psi }_{(1,0)}^1={\left(\begin{array}{cc}-Q(B+1)/S(B+1)+1-M& 1\end{array}\right)}^T\text{and}{\psi }_{(1,0)}^2={\left(\begin{array}{cc}1& 0\end{array}\right)}^T.$$

$$J(M,0)=\left(\begin{array}{cc}-M(M-1)(C+M)& -MQ(B+M)(C+M)\\ 0& S(B+M)(C+M)\end{array}\right),$$

$${\psi }_{(M,0)}^1={\left(\begin{array}{cc}MQ(B+M)/(M(1-M)-S(B+M))& 0\end{array}\right)}^T\text{and}{\psi }_{(M,0)}^2={\left(\begin{array}{cc}1& 0\end{array}\right)}^T.$$

$$J(0,C)=\left(\begin{array}{cc}-C(M+BQC)& 0\\ BCS& -BCS\end{array}\right),$$

$$J(u_i,u_i+C)=\left(\begin{array}{cc}{J}_{12}(u_i,u_i+C)& -Q{u}_i(u_i+B)(u_i+C)\\ S(u_i+B)(u_i+C)& -S(u_i+B)(u_i+C)\end{array}\right),$$

$$\det (J(u_i,u_i+C))=$$

$$-1+Q(B+C)),$$

$$\mathrm{tr}(J(u_i,u_i+C))=$$

$$\begin{array}{r}f(u_i)=\frac{(u_i(M-2u_i-Q{u}_i+1)-S(B+u_i))}{u_{i}Q}.\end{array}$$

$${\lambda }_{P_1}^{1,2}=-\frac{(u_1+C)\left(BS+u_1(-1-M+S+CQ+u_1(2+Q))\right)\pm \sqrt{p_1(u_1)}}{2}$$

$${\psi }_{P_1}^{1,2}=\left(\begin{array}{c}\frac{BS+u_1(1+M+S-CQ-u_1(2+Q))\pm \sqrt{p_2(u_1)}}{2S(u_1+B)}\\ 1\end{array}\right)$$

$$p_1(u_1)=$$

$$+u_1(2+Q)))+(BS+u_1(-1-M+S+CQ+(2+Q)u){)}^2,$$

$$p_2(u_1)=$$

$$+((1+M-CQ{)}^2-2(-1-M+CQ+B(2+Q)+2Q{u}_1)S$$

$$+S^2)u_1^2-2(2+Q)(1+M+S-CQ)u_1^3+(2+Q{)}^2{u}_1^4).$$

$$-M+2u_1(1+Q)-1+Q(B+C)$$

$$-M+2u_2(1+Q)-1+Q(B+C)=\sqrt{\Delta }>0.$$

Full text

††journal: arXiv

1 Introduction

The goal of analysing the dynamics of complex ecological systems is to better describe the different interactions between species, to understand their longterm behaviour, and to predict how they will respond to management interventions [27, 37]. Current predator-prey dynamics studies often use nonlinear mathematical models to describe the species’ interactions and answer these questions. These models aim to be representative of real natural phenomena and they should capture the essentials of the dynamics. However, new theoretical, empirical, and observational research in ecology is revealing species’ interactions to be much more complicated than previous models admit [23, 24, 43, 47]. Moreover, it is becoming increasingly apparent that our understanding of ecosystem dynamics will depend [50], to some extent, on the particular nature of these interaction processes, such as the functional response or predation rate [12, 45, 47].

The standard approach for using models to understand ecological systems is to construct a model from first principles, and then compare species’ abundance timeseries to the predictions from those models. However, this approach becomes more difficult when we add additional nuance to standard models, making them more complex, more nonlinear, and more difficult to parameterise. For instance, Graham and Lambin [35] showed that field-vole (Microtus agrestis) survival can be affected by reducing the weasel predation. They also demonstrated that weasel proportion was suppressed in summer and autumn, while the voles (Microtus agrestis) population always declined to low density. However, they argued that the underlying model was too hard to study due to the large number of parameters. Some ecologists have attempted to resolve this issue by applying qualitative approaches, which make few assumptions about the models’ functional forms or parameters [10, 32, 42]. However, we can also approach the problem by trying to understand the topology of the associated dynamical system, rather than specific trajectories [41]. Such a topological approach may offer general and global insights into the behaviour of the system without requiring accurate parameter estimates.

The phenomena described above can be observed in predator-prey theory. The original Lotka-Volterra predator prey models [34] were straightforward, with simple functional forms for species’ growth and interactions. Empirical observations required successive changes to these assumptions, leading inter alia to the May–Holling–Tanner model, which is itself a special case of the Leslie–Gower predator-prey model [28, 44]. The May–Holling–Tanner model is described by an autonomous two-dimensional system of ordinary differential equations, where the equations for the growth of the predator and prey are logistic-type functions, where the predator carrying capacity is a prey dependent [1, 19, 47]. The functional response describing the predation is Holling Type I, which, for instance, models filter feeders where searching for food can occur at the same time that the species processes the food [33]. A Holling Type I response function corresponds to a linear increasing function in the prey $H(x)=qx$ [26]. This type of functional response is also called Lotka-Volterra type. In particular, the model is given by

[TABLE]

Here, $x(t)$ and $y(t)$ represent the proportion of the prey respectively predator population at time $t$; $r$ is the intrinsic growth rate for the prey; $s$ is the intrinsic growth rate for the predator; $q$ is the per capita predation rate; $K$ is the prey carrying capacity; $n$ is a measure of the quality of the prey as food for the predator; and $\widetilde{K}(x)=nx$ is the prey dependent carrying capacity of the predator.

However, even model (1) does not take into account that some predators act as generalists [3, 18, 25]. For instance, weasels (Mustela nivalis) in the boreal forest region in Fennoscandia can switch to an alternative food source, although its population growth may still be limited by the fact that its preferred food, voles (Microtus agrestis), are not available abundantly [24, 30, 47]. This characteristic can be modelled by modifying the prey dependent carrying capacity of the predator [9]. That is, in (1)

[TABLE]

where we assumed that the alternative food source is constant, which in turns means that the predator proportion is small in compared to the alternative food source. Model (1) with (2) was studied in [5, 9]. It was shown that, in comparison to the original model (1), there is an extra equilibrium point on the $y$-axis corresponding to the extinction of the prey but not the predator. Moreover, the nodependentn-negative parameter $c$ desingularises the origin of system (1).

Another effect that is not incorporated in (1) is the Allee effect [2]. The Allee effect corresponds to a density-dependent phenomenon in which fitness growth initially increases as population density increases [11, 16, 31, 46]. This effect is usually modelled by adding a factor $(x-m)$ to the logistic function where $m$ is the minimum viable population [2, 48, 53, 54] and $0<m<K$. With the Allee effect included in (1)

[TABLE]

For $0<m<K$, the per-capita grow rate of the the prey population with the Allee effect included is negative, but increasing, for $x\in[0,m)$, and this is referred to as the strong Allee effect. When $m\leq 0$, the per-capita growth rate is positive but increases at low prey population densities and this is referred to as the weak Allee effect [11, 15]. Additionally, the Allee effect can also refer to a decrease in per capita fertility rate at low population densities or a phenomenon in which fitness, or population growth, increases as population density increases [2, 16, 31, 46]. For instance, Ostfeld and Canhan [39] found that the stabilisation of vole (Microtus agrestis) populations in southeastern New York depends on the variation in reproductive rate and recruitment of the population. This effect is referred to as the multiple Allee effect [11], sometimes also called the double Allee effect [4, 22]. To incorporate this multiple Allee effect in (3) $L_{m}(x)$ is replaced by

[TABLE]

Here, $b$ is the non-fertile prey population and $0<m<K$ [2, 48, 53, 54]. The per-capita growth rate for the logistic growth function, strong and weak Allee effect; and the multiple Allee effect are shown in Figure 1. We observe that the multiple Allee effect reduces the region of depensation, that is, the region where the per-capita growth rate is positive and growing, when compared to the strong Allee effect. This effect can be generated by the reduction of the probability of fertilisation at lower population density [33]. This reduction commonly occurs in plants such as Diplotaxis erucoides, Banksia goodii and Clarkia concinna [33]. IIn particular, the depensation region for the multiple Allee effect is given by $(m,x_{1})$ with

[TABLE]

and for the strong Allee effect by $(m,x_{2})$ with

[TABLE]

and $x_{1}\leq x_{2}$ for all values of $b$, see Figure 1.

When the alternative food (2) and the multiple Allee effect (4) are included in the modified May–Holling–Tanner model (1) it becomes

[TABLE]

The aim of this manuscript is to study the dynamics of (7) and, in particular, understanding the change in dynamics the multiple Allee effect and the alternative food source causes. Additionally, models (1) and (7) without alternative food sources revealed that there exists a subset of the system parameters where the predator and prey population goes extinct [36]. However, these models assumed different dynamics at low abundance, and the absence of an alternative prey. We find that the alternative food source desingularises the origin and it prevents the extinction of the predator populations. Moreover, we study the basins of attraction of the stable positive equilibrium point(s) by modifying the predation rate $q$ and/or the alternative food source $c$. Moreover, we will show that the addition of the alternative food source and the multiple Allee effects will lead to complex dynamics, and different types of bifurcations such as Hopf bifurcations, homoclinic bifurcations, saddle-node bifurcations and Bogadonov-Takens bifurcations. This manuscript also extends the properties of the May–Holling–Tanner model with multiple Allee effects studied in [36] that is (7) with $c=0$ by showing the impact of the inclusion of alternative food sources for predators. In addition, it complements the results of the May–Holling–Tanner model considering only alternative food for the predator studied in [7, 20] and the model considering only a single Allee effect on the prey and no alternative food for the predator studied in [49, 52]. Model (1) with functional response Holling type II, i.e. $H(x)=qx/(x+a)$, was studied in [28] where the authors showed that there is a region a parameter space where the unique positive equilibrium point is globally asymptotically stable. This model was also studied in [44] where the authors proved the existence of two limit cycles and the species can thus coexist and oscillate.

The basic properties of the model are briefly described in Section 2. In Section 3 we prove the stability of the equilibrium points and give the conditions for the different types of bifurcations. In addition, we discuss the impact changing the predation rate or the alternative food source has on the basins of attraction of the positive equilibrium point in system (7). We further discuss the results and give the ecological implications in Section 4.

2 The Model

Following [5, 7, 8, 13], we introduce dimensionless variables $(u,v,\tau)$ by the function $\varphi:\breve{\Omega}\times\mathbb{R}\rightarrow\Omega\times\mathbb{R}$, where $\varphi(u,v,\tau)=(x,y,t)=(Ku,nKv,\tau(u+c/(nK))(u+b/K)/r)$, $\Omega=\{(x,y)\in\mathbb{R}^{2},x\geq 0,y\geq 0\}$ and $\breve{\Omega}=\{(u,v)\in\mathbb{R}^{2},u\geq 0,v\geq 0\}$. Additionally, we set $B:=b/K$, $C:=c/(nK)$, $M:=m/K\in(0,1)$, $S:=s/r$ and $Q:=qnK/r$, such that $(M,B,C,S,Q)\in\Pi=(0,1)\times\mathbb{R}^{4}_{+}$. This way, we convert (7) to a topologically equivalent nondimensionalised model given by

[TABLE]

The mapping $\varphi$ is a diffeomorphism which preserve the orientation of time since $\det\varphi(u,v,\tau)=nK^{2}u(u+b/K)/r>0$ [14]. Therefore, system (8) is topologically equivalent to system (7) in $\Omega$. Furthermore, system (8) is of Kolmogorov type since $du/d\tau=uR(u,v)$ and $dv/d\tau=vW(u,v)$, with $R(u,v)=(u+C)(u-M)(1-u)-Q(u+C)(u+B)v$ and $W(u,v)=S(u+B)(u-v+C)$. The $u$-nullcline of system (8) in $\breve{\Omega}$ is $v=(u-M)(1-u)/Q(u+B)$, while the $v$-nullcline in $\breve{\Omega}$ is $v=u+C$. Hence, the equilibrium points in $\breve{\Omega}$ for the system (8) are $(0,0)$, $(M,0)$, $(0,C)$ $(1,0)$ and $(u^{*},v^{*})$, where $u^{*}$ is determined by the roots of the following equation

[TABLE]

We observe that $\lim\limits_{u\rightarrow\pm\infty}p(u)=-\infty$ and $\lim\limits_{u\rightarrow\pm\infty}d(u)=\infty$. Hence, $p(u)$ can intersect $d(u)$ in the first quadrant in two points; one point or not at all, see Figure 2.

The solutions of the equation (9) are given by

[TABLE]

such that $M<u_{1}\leq u_{3}\leq u_{2}<1$, where $u_{3}=(1+M-Q(B+C))/(2(1+Q))$. That is, if (9) has two real-valued solutions then these solutions are in the interval $(M,1)$.

Varying the parameters $Q$ and $C$ modifies the value of $\Delta$ and hence the number of equilibrium points in the first quadrant. Specifically:

1. (a)

System (8) has no positive equilibrium points if $\Delta<0$; 2. (b)

System (8) has two positive equilibrium points $P_{1,2}=(u_{1,2},u_{1,2}+C)$ if $\Delta>0$; and 3. (c)

System (8) has one positive equilibrium point $P_{3}=((u_{3},u_{3}+C))$ (order two) if $\Delta=0$.

3 Main Results

In this section, we discuss the stability of the equilibrium points and their bifurcations.

Theorem 3.1**.**

The region $\Phi=\{(u,v),\ 0\leq u\leq 1,\ 0\leq v\leq 1+C\}$ is an invariant region and attracts all trajectories starting in the first quadrant.

Proof.

We follow the proof of [6] where a Holling–Tanner model with strong Allee effect is studied. The main difference between the system studied in [6] and system (8) is that the equilibrium points are located in $\Phi=\{(u,v),\ 0\leq u\leq 1,\ 0\leq v\leq 1+C\}$. However, the invariant region $\Gamma$ is the same invariant region showed in [6] and the system is also a Kolmogorov type. Therefore, trajectories enter into $\Gamma$ and remain in $\Gamma$, see Figure 3. Moreover, trajectories inside $\Lambda=\{(u,v),\ u>1,\ 0<v<u+C\}$ enter into $\Phi$ or the region $\Theta=\{(u,v),\ u>1,\ v\geq u+C\}$ since $du/d\tau<0$ and $dv/d\tau>0$, see $\Lambda$ and $\Theta$ in Figure 3. The $u$-component of trajectories in $\Theta$ are non-increasing as time increases and then these trajectories enter into $\Gamma\backslash\Phi$. As a result, all trajectories starting outside $\Gamma$ enter into $\Gamma$ and end up in $\Phi$ since if $u<1+C$, then $dv/d\tau<0$. ∎

3.1 Nature of equilibrium points

The Jacobian matrix $J(u,v)$ of system (8) is

[TABLE]

where $J_{11}(u,v)=((1-u)(u-M)-Q(u+B)v)(2u+C)$, $J_{12}(u,v)=(M-2u+1-Qv)(u+C)u$ and $d(u)$ is defined in (9).

Lemma 3.1**.**

The equilibrium points $(0,0)$ and $(1,0)$ are saddle points.

Proof.

The Jacobian matrix evaluated at $(0,0)$ gives

[TABLE]

with eigenvalues $\lambda_{(0,0)}^{1}=-CM<0$ and $\lambda_{(0,0)}^{2}=BCS>0$ and eigenvectors

[TABLE]

Similarly, the Jacobian matrix evaluated at $(1,0)$ gives

[TABLE]

with eigenvalues $\lambda_{(1,0)}^{1}=S(C+1)(B+1)>0$ and, since $0<M<1$, $\lambda_{(1,0)}^{2}=(M-1)(C+1)<0$. The associated eigenvectors are

[TABLE]

Thus, it follows that $(0,0)$ and $(1,0)$ are a saddle points in system (8). ∎

Lemma 3.2**.**

The equilibrium point $(M,0)$ is a repeller.

Proof.

The Jacobian matrix evaluated at $(M,0)$ gives

[TABLE]

with eigenvalues $\lambda_{(M,0)}^{1}=M(1-M)(C+M)>0$ and $\lambda_{(M,0)}^{2}=M(1-M)(C+M)>0$ and eigenvectors

[TABLE]

It follows that $(M,0)$ is a hyperbolic repeller in system (8). ∎

Lemma 3.3**.**

If $\Delta\geq 0$ (10), then the equilibrium point $(0,C)$ is a local attractor. Moreover, if $\Delta<0$ (10), then $(0,C)$ is a global attractor (for positive initial conditions).

Proof.

The Jacobian matrix evaluated in the point $(0,C)$ is

[TABLE]

with eigenvalues $\lambda_{(0,C)}^{1}=-C(BCQ+M)<0$ and $\lambda_{(0,C)}^{2}=-BCS<0$ and eigenvectors $\psi_{(0,C)}^{1}=\begin{pmatrix}-(M+B(CQ-S))/BS&1\end{pmatrix}^{T}$ and $\psi_{(0,C)}^{2}=\begin{pmatrix}0&1\end{pmatrix}^{T}$. It follows that $(0,C)$ is local attractor in system (8). Moreover, if $\Delta<0$ (10), then $(0,C)$ is the only stable equilibrium point in $\Phi$. Hence, by the Poincaré–Bendixson Theorem $(0,C)$ is the unique $\omega$-limit for all trajectories starting in the first quadrant, since by Theorem 3.1 all positive solutions are bounded and eventually end up in $\Gamma$, see Figure 4. ∎

Next, we consider system parameters values such that system (8) has two equilibrium points in the first quadrant, that is, we assume $\Delta>0$ (10). These equilibrium points lie on the line $v=u+C$ such that $Qh(u)=g(u)$ (9) and $J_{11}=0$ (11). Hence, the Jacobian matrix (11) at these equilibrium points simplifies to

[TABLE]

with $J_{12}(u_{i},u_{i}+C)=(M-2u_{i}+1-Q(u_{i}+C))(u_{i}+C)u_{i}$, $i=1,2$ and $u_{i}$ given in (10). The determinant and the trace of the Jacobian matrix (12) are:

[TABLE]

where

[TABLE]

Thus, the sign of the determinant depends on the sign of $-M+2u_{i}(1+Q)-1+Q(B+C)$ and the sign of the trace depends on the sign of $f(u_{i})-C$. Moreover, the eigenvalues of the Jacobian matrix of system (8) evaluate at $P_{1}=(u_{1},u_{1}+C)$ are

[TABLE]

and eigenvectors

[TABLE]

with

[TABLE]

Note that the first element of $\psi_{P_{1}}^{1}>0$ since $1+M-CQ>u_{1}(2+Q)$. Similarly, it turns out that $\psi_{P_{1}}^{2}>0$. This gives the following results.

Lemma 3.4**.**

If $\Delta>0$ (10), then the equilibrium point $P_{1}$ is a saddle point.

Proof.

Evaluating $-M+2u(1+Q)-1+Q(B+C)$ at $u_{1}$ gives:

[TABLE]

Hence $\det(J(P_{1}))<0$ and $P_{1}$ is thus a saddle point, see Figure 5. ∎

Lemma 3.5**.**

If $\Delta>0$ (10), then the equilibrium point $P_{2}$ is:

1. (i)

a repeller if $0<C<C_{H}=f(u_{2})$; and 2. (ii)

an attractor if $C>C_{H}$,

with $f$ defined in (13).

Proof.

Evaluating $-M+2u(1+Q)-1+Q(B+C)$ at $u_{2}$ gives:

[TABLE]

Hence $\det(J(P_{2}))>0$. Evaluating $f(u)-C$ at $u=u_{2}$ gives

[TABLE]

Therefore, the sign of the trace, and thus the behaviour of $P_{2}$, depends on the parity of $f(u_{2})-C_{H}$, see Figure 5. ∎

If $\Delta>0$ and $C>C_{H}$, then system (8) has two stable equilibrium points $(0,C)$ and $P_{2}$. Furthermore, if $C=C_{H}$, then ${\rm{tr}}(J(P_{2}))=0$ and $P_{2}$ undergoes a Hopf bifurcation [14].

Finally, if $\Delta=0$ (10) then the equilibrium points $P_{1}$ and $P_{2}$ collapse and system (8) has a unique equilibrium point in the first quadrant.

Lemma 3.6**.**

If $\Delta=0$ (10), then the equilibrium point $P_{3}$ is:

1. (i)

a saddle-node attractor if $C>C_{SN}=f(u_{3})$; and 2. (ii)

a saddle-node repeller if $C<C_{SN}$,

with $f$ defined in (13).

Proof.

Evaluating $-M+2u(1+Q)-1+Q(B+C)$ at $u=u_{3}$ gives:

[TABLE]

Hence $\det(J(P_{3}))=0$. Evaluating $f(u)-C$ at $u=u_{3}$ gives

[TABLE]

Therefore, the sign of the trace, and thus the behaviour of $P_{3}$, depends on the parity of $f(u_{3})-C$, see Figure 6. ∎

3.2 Bifurcation analysis

In this section we present some of the possible bifurcation scenarios when $\Delta=0$ (10) in system (8).

Theorem 3.2**.**

If $\Delta=0$ (10), then by changing $Q$ system (8) experiences a saddle-node bifurcation at the equilibrium point $P_{3}$.

Proof.

In order to prove the saddle-node bifurcation at $P_{3}=(u_{3},u_{3}+C)$ with $u_{3}=(1+M-Q(B+C))/(2(Q+1))$ we follow the Sotomayor’s Theorem [40]. First, if we consider $\Delta=0$ then system (8) has one positive equilibrium point $P_{3}=(u_{3},u_{3}+C)$. Moreover, in Lemma 3.6 we showed that if $\Delta=0$, then $\det(J(P_{3}))=0$. So, $\lambda=0$ is an eigenvalue of the Jacobian matrix $J(P_{3})$ with eigenvector $U=\begin{pmatrix}1&1\end{pmatrix}^{T}$. Furthermore, we denote $W$ as the eigenvector corresponding to the eigenvalue $\lambda=0$ of the Jacobian matrix $J(P_{3})^{T}$

[TABLE]

The vector form of system (8) is given by

[TABLE]

then differentiating $F$ with respect to the bifurcation parameter $Q$ at $P_{3}$ gives

[TABLE]

with $-u_{3}(u_{3}+B)(u_{3}+C)^{2}=\dfrac{1}{16(1+Q)^{4}}(1+M-BQ-CQ)(1+M+BQ-CQ+2B)(-2C-M+BQ-CQ-1)^{2}.$

Therefore,

[TABLE]

since we assumed $\Delta=0$ and $u_{3}>0$.

Next, we analyse the expression $W\cdot[D^{2}F(P_{3};Q)(U,U)]$. Therefore, we first compute the Hessian matrix at the equilibrium point $P_{3}$

[TABLE]

Hence, since $M\in(0,1)$, we get

[TABLE]

again since we assumed $\Delta=0$ and $u_{3}>0$.

Thus, the conditions of Sotomayor’s Theorem [40] are satisfied. Hence, system (8) experiences a saddle-node bifurcation at the equilibrium point $P_{3}$. ∎

If $\Delta=0$ (10) and $C=C_{SN}=f(u_{3})$, then the equilibrium points collapse and system (8) has one positive equilibrium point $P_{3}$. This equilibrium point is a cusp point given that a non-degeneracy condition is met. To show that, we first translate the equilibrium point $P_{3}=(u_{3},u_{3}+C)$ to the origin by setting $X=u-u_{3}$ and $Y=v-u_{3}-C$ and expand system (8) in a power series around the origin. System (8) can now be written as

[TABLE]

Making the affine transformation

[TABLE]

system (14) becomes

[TABLE]

By Lemma $3.1$ presented in [51] we obtain an equivalent system of (15) as follows

[TABLE]

with $L_{20}=S(S+BQ)(S+CQ)(Q^{2}+S^{3}+BQS^{2}+CQS^{2}+BCQ^{2}S)/Q^{4}>0$ since all the parameters are positive and $L_{11}=(4S^{2}(S+BQ)(S+CQ)(Q+1)(3B-C+M+2BQ-2CQ+1)+Q^{2}(Q(M-BQ-CQ+1)(2C+M-BQ+CQ+1)-2(2C+2M-2BQ^{2}+CQ^{2}-2C^{2}Q+3B^{2}Q^{2}+2B^{2}Q^{3}-2C^{2}Q^{2}-C^{2}Q^{3}+2CM-4BQ+CQ+M^{2}-3BCQ^{2}-BCQ^{3}-2BMQ^{2}+CMQ^{2}-4BCQ-4BMQ+CMQ+1)))/(4Q^{2}(Q+1)^{2})$. If $L_{11}\neq 0$, then $P_{3}$ is a cusp point of codimension two by the result presented in [40].

This is also a necessary condition for system (8) to undergo a Bogdanov-Takens bifurcation [40]. One needs to vary two parameters in order to encounter this bifurcation in a structurally stable way and to describe all possible qualitative behaviours nearby [14, 40]. The proof of a Bogdanov-Takens bifurcation can be obtained by following [29] and [51]. In these articles, the authors showed that their system undergoes to a Bogdanov–Takens bifurcation by unfolding the system around the cusp of codimension two. Moreover, by using a series of normal form transformations one can check the non-degeneracy condition. Nowadays, there are several computational methods to find Bogdanov-Takens points. These methods are implemented in software packages such as MATCONT [17]. Figure 8 illustrates the Bogdanov-Takens bifurcation which was detected with MATCONT in the $(Q,C)$-plane with parameter values $(M,B,S)=(0.05,0.1,0.071080895)$ fixed.

3.3 Basins of attraction

In this section, we analyse the impact of the modifications of the parameters $C$ and $Q$ on the basins of attraction of the stable equilibrium points of system (8). Note that the parameter $C=c/(nK)$ of system (8) is equivalent to the alternative food source $c$ in system (7) since the function $\varphi$ is a diffeomorphism preserving the orientation of time. Similarly, the parameter $Q=qnK/r$ of system (8) is equivalent to the predation rate $q$ in system (7). In particular, we consider the system parameters $(B,M,S)=(0.1,0.1,0.157)$111Note that changing $B$ instead of $Q$ has the same qualitative effect on the basin of attraction, see right pane of Figure 8. and vary $Q$ and $C$. For $Q$ and $C$ not too big system (8) has two positive equilibrium points, namely $P_{1}$ and $P_{2}$. The equilibrium points on the axis and $P_{1}$ do not change stability proven in Lemmas 3.1, 3.2, 3.3 and 3.4, while, $P_{2}$ can be stable or unstable.

In order to study the basins of attraction of the equilibrium points $(0,C)$ and $P_{2}$ we use the same notation for the (un)stable manifold of the equilibrium point $P_{1}$ as used in [5, 6]. That is, we defind $W^{u,s}_{\nearrow}(P_{1})$ as the branch of the (un)stable manifold of $P_{1}$ that goes up to the right and $W^{u,s}_{\swarrow}(P_{1})$ as the branch of the (un)stable manifold of $P_{1}$ that goes down to the left. The branch $W^{s}_{\nearrow}(P_{1})$ is connected with $(M,0)$ and $W^{u}_{\swarrow}(P_{1})$ is connected with $(0,C)$ since the nullclines form a bounding box from which trajectories cannot leave. Furthermore, everything in between of these two branches and the $x$-axis also asymptotes to the equilibrium point $(0,C)$. Therefore, the stable manifold of the saddle point $P_{1}$ acts as a separatrix curve between the basins of attraction of $P_{2}$ (when it is stable) and $(0,C)$, see Figure 9.

Considering the invariant region $\Phi$, there are qualitatively six different cases for the boundaries of the basins of the equilibrium points $P_{2}$ and $(0,C)$, then we get:

1. (i)

For $C<C_{HOM}$ such that the equilibrium point $P_{2}$ in system (8) is stable, see Lemma 3.5 (since $C_{HOM}<C_{H}$). For $C$ small enough $W^{s}_{\swarrow}(P_{1})$ intersects the boundary of $\Phi$. Hence, it forms a separatrix curve in $\Phi$, see panel ($i$) in Figure 9. In addition, by increasing $C$ the stable manifold of $P_{1}$ connects first with $(1,0)$ and then with $(M,0)$, again forming the separatrix curve, see panels ($ii$) and ($iii$) in Figure 9. 2. (ii)

For $C=C_{HOM}$, then $W^{s}_{\swarrow}(P_{1})$ connects with $W^{u}_{\nearrow}(P_{1})$, therefore it form a homoclinic curve. Which is the separatrix curve between the basins of attraction of $(0,C)$ and $P_{2}$, see panel ($iv$) in Figure 9. 3. (iii)

For $C_{HOM}<C<C_{H}$, there is an unstable limit cycle surrounding $P_{2}$ which acts as a separatrix curve between the basins of attraction of $P_{2}$ and $(0,C)$. This limit cycle is created around $P_{2}$ via the Hopf bifurcation at $C=C_{H}$ [21] and terminates via a homoclinic bifurcation at $C=C_{HOM}$, see panel ($v$) in Figure 9. 4. (iv)

For $C_{H}<C<C_{SN}$, the equilibrium point $P_{2}$ is unstable, see Lemma 3.5, and $(0,C)$ is globally asymptotically stable. Hence, $\Phi$ is the basin of attraction of $(0,C)$, see panel ($vi$) in Figure 9. 5. (v)

For $C=C_{SN}$, the equilibrium points $P_{1}$ and $P_{2}$ collapse, see Lemma 3.6. Hence, $\Phi$ is the basin of attraction of $(0,C)$, see panel ($vii$) in Figure 9. 6. (vi)

For $C_{SN}<C$, system (8) dose not have positive equilibrium points, see Lemma 3.3. Hence, $\Phi$ is also the basin of attraction of $(0,C)$, see panel ($viii$) in Figure 9.

4 Conclusions

In this manuscript, a modified May–Holling–Tanner predator-prey model with multiple Allee effects for the prey and alternative food sources for the predators was studied. Using a diffeomorphism, we transformed the modified May–Holling–Tanner predator-prey model to a topologically equivalent system, system (8). Subsequently, we analysed system (8) and we proved that the equilibrium points $(0,0)$ and $(1,0)$ are saddle points, $(M,0)$ is a repeller and $(0,C)$ is an attractor for all parameter values, see Lemmas 3.1, 3.2 and 3.3. Additionally, there exist at most two positive equilibrium points, one of them, $P_{1}$, is a saddle point, while the other, $P_{2}$, can be an attractor or a repeller, depending on the trace of its Jacobian matrix. Both equilibrium points can collapse having conditions for a saddle node bifurcations and cusp point [51] (Bogdanov-Takens bifurcation). We also showed the existence of a homoclinic curve, determined by the stable and unstable manifolds of the equilibrium point $P_{1}$ enclosing the second equilibrium point $P_{2}$. When the homoclinic breaks it creates a non-infinitesimal limit cycle, see Lemmas 3.4, 3.5 and Figure 8.

Moreover, by choosing the bifurcation parameters $(C,Q)$, or $(B,C)$, we have obtained significant bifurcation diagrams, see Figure 8. It follows that – for a large nondimensionalised predation $Q$ and a small nondimensionalised proportion of alternative food $C$ – co-existence is expected. Similarly, when the proportion of nondimensionalised alternative food $C$ is bigger than the proportion of nondimensionalised predation $Q$ co-existence is expected. The bifurcation diagrams and associated phase planes, see Figure 9, also shows that there exists complexity for system (8) including the collision of the equilibrium points leading to different type of bifurcation.

Since the function $\varphi$ is a diffeomorphism preserving [36]the orientation of time, the dynamics of system (8) are topologically equivalent to the dynamics of system (7). Hence, the parameters $(C,Q)$ impact the number of equilibrium points of system (8) in the first quadrant and change the behaviour of the system, and, as $C=c/(nK)$ and $Q=qnK/r$, the system parameters $(c,n,k,q,r)$ will thus impact the behaviour of system (7). Therefore, self-regulation depends on the values of these parameters. For instance, keeping all parameters fixed, but increasing the alternative food source $c$, one expects to see a change in behavior and dynamics similar to the one shown in Figure 8 and 9. All these results show that dynamical behavior of system (7) becomes more complex under the modification of the system parameters when compared to the May–Holling–Tanner model with the strong and weak Allee effect (3) studied in [36].

In Figure 1 we showed that the inclusion of a multiple Allee effect changes the shape of the per-capita growth of the prey, and, in particular, reduces the region of depensation. Moreover, we can see in Figure 10 that there exist a critical non-fertile prey population $b_{cr}$ for which the basin of attraction of the equilibrium point $p_{2}$ of system (7) is smaller than the basin of attraction of the related $p_{2}$ of system (1) considering an alternative food source (2) and with a strong Allee effect (3). Note that the non-fertile prey population of $60\%$ is realistic. For instance, Monclus et al. [38] studied the impact of the different population densities on stthe marmot reproduction. This study used the proportion of fertile female adults which fluctuated between $2.13$ and $19.15\%$ of the total population density. Moreover, we can also conclude that the basin of attraction of the stable positive equilibrium point $p_{2}$ increases when we reduce the depensation in the model.

Finally, the techniques used in this manuscript show that there is a strong connection between the analysis of the manifold and the basins of attraction of the equilibrium points. This analysis can be applied in population dynamics in order to predict the behaviour in models where there is variation in the non-fertile population. Moreover, we showed that the combination of different techniques such as numerical simulations and bifurcation analysis can be very useful for showing the temporal dynamics in predation interaction.