Bifurcation analysis of a prey-predator model with predator intra-specific interactions and ratio-dependent functional response

TL;DR

This paper analyzes a predator-prey model with predator intraspecific interactions and ratio-dependent response, revealing complex bifurcation behaviors and the effects of parameter changes on system stability.

Contribution

It provides a detailed bifurcation analysis of the Bazykin predator-prey model with novel insights into stability and dynamics.

Findings

Existence and stability of two interior equilibrium points.

Identification of multiple bifurcation types including saddle-node, Hopf, homoclinic, and Bogdanov-Takens.

Numerical simulations showing the impact of predator consumption rate and efficiency on system dynamics.

Abstract

We study the Bazykin predator-prey model with predator intraspecific interactions and ratio-dependent functional response and show the existence and stability of two interior equilibrium points. We prove that the model displays a wide range of different bifurcations, such as saddle-node bifurcations, Hopf bifurcations, homoclinic bifurcations and Bogdanov-Takens bifurcations. We use numerical simulations to further illustrate the impact changing the predator per capita consumption rate has on the basin of attraction of the stable equilibrium points, as well as the impact of changing the efficiency with which predators convert consumed prey into new predators.

Equations

$$\frac{dN}{dt}$$

$$\frac{dP}{dt}$$

$$\phi :\bar{\Omega }\times \mathbb{R}\to \Omega \times \mathbb{R}\text{where}\left(u,v,\tau \right)=\phi \left(N,P,t\right)=\left(\frac{N}{K},\frac{aP}{K},\frac{rKt}{N+aP}\right).$$

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

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

$$u=0\text{and}v\left(Q-1+u\right)=u\left(1-u\right),$$

$$v=0\text{and}v=\frac{-\left(M+Nu\right)+\sqrt{{\left(M+Nu\right)}^2+4Nu\left(C-M\right)}}{2N}.$$

$$\left(C+NQ\right)u^2+{\Sigma }_{1}u+\left(Q-1\right){\Sigma }_2=0$$

$${\Sigma }_1$$

$${\Sigma }_2$$

$$N\left(C+NQ\right)v^2+{\Sigma }_{3}v+\left(C-M\right){\Sigma }_2=0$$

$${\Sigma }_3=2NQ(M-C)+C(M+N)=-2N{\Sigma }_2+C(M-N)=-N{\Sigma }_1+(QN+C)(M-N).$$

$$u_{1,2}=\frac{-{\Sigma }_1\pm Q\sqrt{\Delta }}{2\left(C+NQ\right)}\text{and}{v}_{1,2}=\frac{-{\Sigma }_3\pm C\sqrt{\Delta }}{2N\left(C+NQ\right)}.$$

$$\Delta ={\left(M-N\right)}^2-4N\left(C\left(Q-1\right)-MQ\right)={\left(M-N\right)}^2-4N{\Sigma }_2,$$

$$2N{\Sigma }_1=(M-N)\left(M+N(2Q-1)\right).$$

$$2{\Sigma }_3=(M-N)\left(2C+N-M\right).$$

$$J(u,v)=\left(\begin{array}{cc}2u+v-Qv-2uv-3u^2& -u(Q+u-1)\\ -v(M-C+Nv)& Cu-3N{v}^2-Mu-2Mv-2Nuv\end{array}\right).$$

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

$$(u,v)\to \left(xy,y\right)\text{and the time rescaling}\tau \to \frac{t}{y}.$$

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

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

$$J(O_{xy})=\left(\begin{array}{ccc}1+M-Q& & 0\\ 0& & -M\end{array}\right)$$

$${\lambda }_1(O_{xy})=1+M-Q\text{and}{\lambda }_2(O_{xy})=-M.$$

$$J(I_x)=\left(\begin{array}{ccc}-(1+M-Q)& & \frac{(Q-C)(N(C-1)-M(1+N)+Q-1)(1+M-Q)}{(1+M-C{)}^3}\\ 0& & \frac{{\Sigma }_2}{(1+M-C)}\end{array}\right)$$

$${\lambda }_1(I_x)=-1-M+Q\text{and}{\lambda }_2(I_x)=\frac{{\Sigma }_2}{1+M-C}.$$

$$(u,v)\to \left(X,XY\right)\text{and the time rescaling}\tau \to \frac{t}{X}$$

$$\frac{dX}{dt}$$

$$\frac{dY}{dt}$$

$$J_{(O_{XY})}=\left(\begin{array}{cc}1& 0\\ 0& C-1-M\end{array}\right)$$

$${\lambda }_1(O_{XY})=1\text{and}{\lambda }_2(O_{XY})=C-1-M.$$

$$J(I_Y)=\left(\begin{array}{ccc}\frac{-{\Sigma }_2}{M-Q+1}& & 0\\ \frac{(C-Q)(C-M-1)(M+N-Q-CN+MN+1)}{(M-Q+1{)}^3}& & 1+M-C\end{array}\right)$$

$${\lambda }_1(I_Y)=\frac{-{\Sigma }_2}{M-Q+1}\text{and}{\lambda }_2(I_Y)=1+M-C.$$

$$J(u,v)=\left(\begin{array}{cc}u(1-2u-v)& u(1-Q-u)\\ v(C-M-Nv)& -v(M+N(u+2v))\end{array}\right).$$

Full text

Keywords— Predator-prey model; Bifurcations; Ratio-dependent functional response; Intraspecific interactions.

1 Introduction

The original Lotka-Volterra predator-prey models [1] are relatively straightforward with simple functional forms for species’ growth and interactions. Empirical observations required successive changes to these assumptions, leading, inter alia, to the Bazykin models [2]. Temporal [3, 5, 6, 7] and spatio-temporal [8, 9, 10, 11] ratio-dependent predator-prey models are becoming of more interest in ecology since these models better describe (when compared to the original models) both the theoretical and experimental predator-prey interactions [12, 13]. Additionally, ratio-dependent predator-prey models are more suitable for predator-prey interactions when predators are seeking, capturing or killing other animals [14, 15, 16]. For instance, Jost and Arditi [17] found that ratio-dependent predator-prey models are more suitable for identifying the mathematical predator-prey functional response from real data.

Gause type predator-prey models present an enrichment paradox [16, 18] which causes an increase in the equilibrium density of the predator but not in the prey, destabilising the community equilibrium. Gause type predator-prey models also present a biological control paradox [19, 20] which refers to the fact that these type of models cannot have a stable low prey density equilibrium. Kuang and Beretta [16] showed that ratio-dependent functional responses can resolve these paradoxes in Gause type predator-prey models.

Bazykin models [2, 21] with ratio-dependent functional response are examples of modified Gause type predator-prey models. The main difference between a standard Gause type model and a Bazykin model is that the latter considers the effect of adding a linear density dependence in the predator equation. The specific Bazykin model considered here is an extension of the model introduced by Arditi and Ginzburg [12] in which the authors showed that a ratio-dependent interaction is more realistic compared to a standard prey-dependent interaction. The Bazykin model is described by an autonomous two-dimensional system of ordinary differential equations, where the equations for the growth of the prey is a logistic-type function [2, 22, 23] and the growth of the predator is a function of the ratio of prey to predator abundance. The functional response is ratio-dependent, in which the consumption rate of prey depends on the numbers of predators and prey, i.e. $H(N,P)=qN/(N+aP)$ [24]. In other words, the ratio dependence is a type of predator dependence in which the functional response depends on the ratio of prey population size to predator population size, not on the absolute numbers of either species [25]. In particular, it is given by

[TABLE]

Here, $N(t)$ and $P(t)$ represent the proportion of the prey, respectively, predator population at time $t$; $r$ is the intrinsic growth rate of the prey; $K$ is the prey carrying capacity; $q$ is the per capita predation rate; $a$ is the amount of prey by which the predation effect is maximum; $c$ is the efficiency with which predators convert consumed prey into new predators; $\mu_{0}$ is the per capita death rate of predators and $\mu_{1}$ is predator death rate by density.

Haque [26] studied a diffusive version of the Bazykin model (1) and showed that, in absence of the diffusion, ratio-dependent predator-prey models are more appropriate for predator-prey interactions when the predators involve serious hunting processes. However, Haque focused on the case when system (1) has only one positive equilibrium point in the first quadrant (see Section 2.1) and the author showed that there is always coexistence of both populations or the extinction of only the predator population. This manuscript extends the diffusive-free analysis studied by Haque [26] and its main aim is to study the bifurcation dynamics of (1). In particular, we aim to understand the change in dynamics the intraspecific interactions and the ratio-dependent functional response causes. We will focus on the most general parameter setting for which system (1) can have up to two positive equilibrium points in the first quadrant and up to two limit cycles. As a result, system (1) supports complex dynamics such as the extinction, the coexistence, and the oscillation, of both populations over time. Furthermore, we will show that the model (1) supports different types of bifurcations such as saddle-node bifurcations, Hopf bifurcations and homoclinic bifurcations. We will also show the impact that changing the predation rate has on the time series behaviour.

The basic properties of the Bazykin model (1) are briefly described in Section 2. In Section 3 we prove the stability of the equilibrium points and in Section 4 we present the conditions for the different types of bifurcations. In addition, we discuss the impact changing the predation rate, i.e. $q$, or the efficiency with which predators convert consumed prey into new predators, i.e. $c$, has on the basins of attraction of the positive equilibrium points of system (1). We further discuss the results and give the ecological implications in Section 5.

2 Preliminary Results

The ratio-dependent Bazykin predator-prey model is given by (1) and we only consider the model in the domain $\Omega=\{\left(N,P\right)\in\mathbb{R}^{2},~{}N>0,~{}P>0\}$. In order to simplify the analysis we follow the nondimensionalisation approach of several other type of predator-prey models [27, 28, 29] and we introduce the dimensionless variables $(u,v,\tau)$ given by

[TABLE]

Let $C:=c$, $M:=\mu_{0}/r$, $N:=\mu_{1}K/(ar)$ and $Q:=q/(ar)$, then (1) transforms into the nondimensionalised equivalent system

[TABLE]

System (3) is defined in $\bar{\Omega}=\{(u,v)\in\mathbb{R}^{2},~{}u\geq 0,~{}v\geq 0\}$ and $\varphi$ (2) is a diffeomorphism for $(u,v)\neq(0,0)$ which preserve the orientation of time since $\det\left(\varphi\left(N,P,t\right)\right)=K^{2}\left(u^{2}+v^{2}\right)/(ar)>0$ [30, 31]. Moreover, system (3) is of Kolmogorov type [32], that is, the axes $u=0$ and $v=0$ are invariant. The $u$ nullclines are

[TABLE]

while the $v$ nullclines are

[TABLE]

Hence, the equilibrium points on the axes for this system are $(0,0)$ and $(1,0)$111Note that the value $1$ of the equilibrium point $(1,0)$ relates to the rescaled prey carrying capacity $K$.. In order to obtain the equilibrium point(s) inside the first quadrant, which we call positive interior equilibrium point(s), we rewrite equation (4) as $v=u(1-u)/(Q-1+u)$ (assuming $Q-1+u\neq 0$). Then, replacing $v$ into $R(u,v)=0$ (3) we obtain:

[TABLE]

with

[TABLE]

Similarly, we can rewrite equation (5) as $u=v(M+Nv)/(C-M-Nv)$ (assuming $C-M-Nv\neq 0$). Then, replacing $u$ into $W(u,v)=0$ (3) we obtain:

[TABLE]

with

[TABLE]

Therefore, the interior equilibrium point(s) $P_{1,2}=(u_{1,2},v_{1,2})$ of system (3) are determined by the solution(s) of the equations (6) and (8) which are given by:

[TABLE]

Here,

[TABLE]

and $u_{1}<u_{2}<1$ and $v_{1}<v_{2}$.

2.1 Number of positive equilibrium points

We describe the different configurations for the solutions of equations (6) and (8) given in (10). In particular, we determine the number of positive interior equilibrium points.

If $C\leq M$, then $u_{i}$ and $v_{i}$, $i=1,2$, cannot be both positive since $u=v(M+Nv)/(C-M-Nv)$. Therefore, there are no positive interior equilibrium points. 2. 2.

If $C>M$ and

1. (a)

$\Sigma_{2}<0$, then $\Delta>0$ and $(C-M)\Sigma_{2}<0$. Hence, the two solutions of Equation (8) are real, but only $v_{2}$ is positive. Equation (6) also has two real solutions and $u_{2}$ is always positive since

1. i.

for $Q>1$ we have that $(Q-1)\Sigma_{2}<0$; and 2. ii.

for $Q\leq 1$ we have that $\Sigma_{1}<0$.

Therefore, system (3) has one positive interior equilibrium point $P_{2}=(u_{2},v_{2})$ in this case. See region 1 in Figure 1. 2. (b)

$\Sigma_{2}=0$, then one of the equilibrium points $(u_{1,2},v_{1,2})$ merges with $(0,0)$. The remaining equilibrium point is in the interior of the first quadrant if both $\Sigma_{1}$ and $\Sigma_{3}$ are negative. Since $\Sigma_{2}=0$, both of these inequalities simplify to $N>M$. Therefore, in this case system (3) has one positive interior equilibrium point $P_{2}=(u_{2},v_{2})$ if $N>M$ and no positive interior equilibrium points if $N\leq M$. 3. (c)

$\Sigma_{2}>0$ and

1. i.

$\Delta<0$, then then system (3) has no positive equilibrium points in the first quadrant, see region 3 in Figure 1. 2. ii.

$\Delta=0$, then both roots merge. So, $u_{1}=u_{2}$ and $v_{1}=v_{2}$. This root of order two is in the first quadrant if and only if both $\Sigma_{1}$ and $\Sigma_{3}$ are negative. Along $\Delta=0$ the expression for $2N\Sigma_{1}$ becomes

[TABLE]

So, since $\Sigma_{2}>0$ implies that $Q>1$, we have that $\Sigma_{1}>0$ along $\Delta=0$ if $M>N$ and $\Sigma_{1}<0$ along $\Delta=0$ if $M<N$. Similarly, the expression for $2\Sigma_{3}$ along $\Delta=0$ is

[TABLE]

So, since $C>M$, we have that $\Sigma_{3}>0$ along $\Delta=0$ if $M>N$ and $\Sigma_{3}<0$ along $\Delta=0$ if $M<N$. Therefore, in this case system (3) has one positive interior equilibrium point $P_{2}=(u_{2},v_{2})$ of order two if $N>M$ and no positive interior equilibrium points if $N\leq M$. 3. iii.

$\Delta>0$ and

1. A.

$N>M$, then $\Sigma_{1}<0$ – which follows directly from the above computation in case 2(c)ii upon replacing the =-sign by the ¡-sign – and $\Sigma_{3}=-2N\Sigma_{2}+C(M-N)<0$. Thus, since $\Sigma_{2}>0$ implies that $Q>1$, we have that both the $u$-intercept and $v$-intercept are positive. Consequently, system (3) has two positive interior equilibrium points. See region 2 in Figure 1. 2. B.

$N\leq M$, then $\Sigma_{1}=2\Sigma_{2}+Q(M-N)\geq 0$ and $\Sigma_{3}\geq 0$ – which follows directly from the above computation in case 2(c)ii upon replacing the =-sign by the $\geq$-sign. Furthermore, we still have that both the $u$-intercept and $v$-intercept are positive. Consequently, system (3) now has no positive interior equilibrium points. See region 4 in Figure 1.

3 Main Results

In this section, we discuss the stability of the equilibrium points of system (3) for $Q>1$ which represents the most general case of the number of equilibrium points. Since $Q:=q/(ar)$, the condition $Q>1$ corresponds to the assumption that the predation rate is larger than the product between the prey’s intrinsic growth rate and the number of prey necessary for getting the maximum predation effect. Note that the case when $Q=1$ and $C>M$ (a subcase of case 2a in Section 2.1) was studied by Haque [26].

Lemma 3.1

The set $\Gamma=\{\left(u,v\right)\in\bar{\Omega},~{}0\leq u\leq 1,~{}v\geq 0\}$ is an invariant region and all solutions of (3) which are initiated in the first quadrant are bounded and eventually end up in $\Gamma$.

Proof. This proof follows the proofs of Arancibia et al. [28, 33]. Since the system (3) is of Kolmogorov type the coordinates axes are invariant [34]. Moreover, if $v=0$ then $du/d\tau=u^{2}(1-u)$ and if $u=0$ then $dv/d\tau=-v^{2}(M+Nv)$. So, any trajectory with initial point on the positive vertical $v$-axis tends to zero and any trajectory with initial point on the positive horizontal $u$-axis tend to $u=1$. Next, setting $u=1$ in the first equation of system (3), we have $du/dt=-Qv<0$ and setting $u>1$, we also have $du/dt<0$ since $1-u<0$. Thus for any initial condition initiated in the first quadrant, the corresponding trajectory enters and remains in $\Gamma$ independent of the sign of $dv/dt$.

To finalise the proof we show that no trajectory in the open region $\Gamma$ converges to infinity as $\tau\to\infty$. To show that solutions are bounded it is enough to find a $v^{*}$ such that $dv/d\tau<0$ for $v\geq v^{*}$ and $0\leq u\leq 1$. The equation $dv/d\tau$ can be written as $dv/d\tau=-Cv^{3}+\mathcal{O}(v^{2},uv^{2})$ which implies that we can always find such as $v^{*}$. Therefore, all trajectories end up in $\Gamma_{1}=\{\left(u,v\right)\in\bar{\Omega},~{}0\leq u\leq 1,~{}0\leq v\leq v^{*}\}$.

3.1 The nature of the equilibrium points

To determine the nature of the equilibrium points we compute the Jacobian matrix $J(u,v)$ of (3)

[TABLE]

3.1.1 Equilibrium points on the axes

Lemma 3.2

The equilibrium point $(1,0)$ is a saddle if $C>M$ and a stable node if $C<M$.

Proof. The result follows direct from the Jacobian matrix (12) evaluated at $(1,0)$

[TABLE]

combined with the Hartman-Grobman theorem.

In other words, if $C<M$ then $(1,0)$ is an attractor for positive initial conditions, see the right panel of Figure 2.

To understand the dynamics near the origin we assume, in addition to the assumption $Q>1$, that $C>M$. That is, we assume that the equilibrium point $(1,0)$ is a saddle point, see Lemma 3.2. In other words, we assume that the efficiency with which predators convert consumed prey into new predators is bigger that the ratio between the per capita death rate of predators and the prey intrinsic growth rate. We divide our $(Q,C)$-parameter space in six regions

• •

Region I: $M+1<C<MQ/(Q-1)$;

• •

Region II: $C>MQ/(Q-1)$ and $1<Q<M+1$;

• •

Region III: $C>M+1$ and $Q>M+1$;

• •

Region IV: $M<C<M+1$ and $1<Q<M+1$;

• •

Region V: $M<C<MQ/(Q-1)$ and $Q>M+1$; and

• •

Region VI: $MQ/(Q-1)<C<M+1$,

see also Figure 3.

Theorem 3.1

If we assume that $Q>1$ and $C>M$, then the origin $(0,0)$ in system (3) is a non-hyperbolic degenerate point [35, 36]. Moreover, a neighbourhood of the origin $(0,0)$ presents six types of topologically different structures in the first quadrant of the phase plane:

• •

a saddle sector and a repelling sector in Region I;

• •

an attracting sector and elliptic sector in Region II;

• •

an elliptic sector in Region III;

• •

a saddle sector in Region IV;

• •

an attracting sector and a saddle sector in Regions V; and

• •

an elliptic sector and a repelling sector in Region VI.

We refer to Figure 3 for exemplary sketches of the dynamics near the origin in the six different regions. Proof. First we observe that setting $u=0$ in system (3) gives $dv/d\tau=-v^{2}(M+Nv)\leq 0$ for $v\geq 0$. That is, any trajectory starting along the $v$-axis converges to the origin $(0,0)$. Also, setting $v=0$ in system (3) gives $du/d\tau=u^{2}(1-u)$ and any orbits starting along the $u$-axis near the origin converges to the equilibrium point $(1,0)$. Next, to analyse the dynamics in a neighbourhood of the origin, we consider the vertical blow-up [32] given by the transformation

[TABLE]

The transformation (13) blows up the origin of system (3) in the entire $x$-axis [32]. The goal is to analyse the equilibrium points on the nonnegative half axis $x\geq 0$, $y=0$, in the new coordinate system

[TABLE]

System (14) has up to two equilibrium points on the nonnegative horizontal $x$-axis: the origin $O_{xy}$ and, for $(C,Q)$ in regions I, II, V or VI, a second equilibrium point $I_{x}=(\mu,0)$ with $\mu=(1+M-Q)/(C-M-1)$. The corresponding Jacobian matrix at $O_{xy}$ is

[TABLE]

with eigenvalues

[TABLE]

Hence, $O_{xy}$ is

• •

a saddle if $1<Q<M+1$, see Regions I, II and IV in Figure 3; and

• •

a stable node if $Q>M+1$, see Regions III, V and VI in Figure 3.

At $I_{x}$ we get the Jacobian matrix

[TABLE]

with eigenvalues

[TABLE]

Hence, $I_{x}$ is

• •

a saddle in regions I and V, see Figure 3.

• •

an stable node in region II, see Figure 3.

• •

an unstable node in region VI, see Figure 3.

We repeat the blow-up procedure to analyse the behaviour of system (3) near the $v$-axis. Thus we consider the horizontal blow-up [32] given by the transformation

[TABLE]

The goal now is to analyse the equilibrium point at the nonnegative half axis $X=0$, $Y\geq 0$, in the new system coordinates which is given by setting (16) in system (3)

[TABLE]

System (17) has again up to two equilibrium points in the nonnegative vertical $Y$-axis. The equilibrium point $(0,0)$ and, for $(C,Q)$ in regions I, II, V or VI, a second equilibrium point $I_{Y}=(0,\mu^{*})$ with $\mu^{*}=(C-M-1)/(1+M-Q)$.

At the origin $O_{XY}=(0,0)$ we get the Jacobian matrix

[TABLE]

with eigenvalues

[TABLE]

Hence, the stability of $O_{XY}$ depends on the parameters $C$ and $M$. In particular, $O_{XY}$ is

• •

a saddle if $M<C<M+1$, see Regions IV, V and VI in Figure 3.

• •

an unstable node if $C>M+1$, see Regions I, II and III in Figure 3.

The Jacobian matrix at $I_{Y}$ is given by

[TABLE]

with eigenvalues

[TABLE]

Upon comparing (15) with (18) we note that the stability of $I_{Y}$ is the same as the stability of $I_{x}$.

Taking the inverse of (13) and (16) the equilibrium points $I_{x}$ and $I_{Y}$ collapse to the origin $(0,0)$ of (3), and the stable and unstable eigenvectors are mapped to the curves $\Gamma^{s}$ and $\Gamma^{u}$ respectively, see the blow-down pictures in the right column of Figure 3. Therefore, it follows that the origin is a non-hyperbolic degenerate point with the properties as stated in Theorem 3.1.

3.1.2 The positive equilibrium points

Next, we consider the stability of the two positive equilibrium points $P_{1,2}$ of system (3) in the interior of $\Gamma$. These equilibrium points lie on the curve $v=u\left(1-u\right)/\left(Q-1+u\right)$ such that $W(u,v)=0$ (3), and they only exist if the system parameters are such that $\Delta>0$ (11). The Jacobian matrix of system (3) at these equilibrium points becomes

[TABLE]

The determinant and the trace of the Jacobian matrix (19) are given by

[TABLE]

This gives the following results.

Theorem 3.2

Let the system parameters of (3) be such that $C>M$, $\Sigma_{2}>0$, $\Delta>0$ (11) and $N>M$ (case 2iiiA of Section 2.1), then there are two interior equilibrium points and the equilibrium point $P_{1}=\left(u_{1},v_{1}\right)$ defined in (10) is a saddle point.

Proof. Since $P_{1}$ is in the first quadrant, we know that $0<u_{1}+v_{1}=\left(N-M-\sqrt{\Delta}\right)/\left(2N\right)$, with $\Delta$ defined in (11). Evaluating the determinant (20) at $P_{1}=\left(u_{1},v_{1}\right)$ gives

[TABLE]

Hence, the equilibrium point $P_{1}$ is a saddle point.

Theorem 3.3

Let the system parameters of (3) be such that we are in cases 2iiiA or 2a of Section 2.1, then system  (3) has at least one positive interior equilibrium point $P_{2}=\left(u_{2},v_{2}\right)$. Therefore, the stability of equilibrium point $P_{2}$ defined in (10) is:

1. (i)

stable if $T\left(u_{2},v_{2}\right)<0$; 2. (ii)

unstable if $T\left(u_{2},v_{2}\right)>0$; and 3. (iii)

a weak-focus if $T\left(u_{2},v_{2}\right)=0$,

where $T$ is defined in (20).

Proof. Evaluating the determinant (20) at $P_{2}=\left(u_{2},v_{2}\right)$ gives

[TABLE]

since $u_{2}+v_{2}=(N-M+\sqrt{\Delta})/(2N)>0.$ The results now follow by analysing the trace (20) of the Jacobian matrix (19) evaluated at $P_{2}$. The function $T\left(u,v\right)$ (20) at the equilibrium point $P_{2}$ is given by

[TABLE]

with

[TABLE]

and

[TABLE]

Therefore, the sign of the trace $T\left(u_{2},v_{2}\right)$ (20), and thus the behaviour of $P_{2}$, depends on the parity of $T_{1}\sqrt{\Delta}+T_{2}$. Evaluating $T_{1}$, $T_{2}$ and $\Delta$ at $(C,M,N,Q)=(0.363,0.16,0.25,1.6)$ gives $T_{1}\sqrt{\Delta}+T_{2}=-0.036091<0$, while $T_{1}\sqrt{\Delta}+T_{2}=0.025983>0$ at $(C,M,N,Q)=(0.363,0.16,0.25,1.8)$. By the continuous dependence of the trace on the parameter $Q$, we get that there exist $Q\in(1.6,1.8)$ for which the trace must be zero.

Next, we discuss the stable manifold $W^{s}(P_{1})$ of the saddle point $P_{1}$. This manifold often acts as a separatrix curve between the basins of attraction of the equilibrium points $(0,0)$ and $P_{2}$, see, for example, Figure 4. Let $W^{u,s}_{\nearrow,\swarrow}(P_{1})$ denote the branch of the (un)stable manifold of $P_{1}$ whose trajectory flows according to the direction of the arrow, see Figure 5 for an example. Following Flores and González‐-Olivares [37] and Aguirre et al. [38], we get from the orientation of the nullclines and the eigenvectors of $P_{1}$ that $W^{s}_{\nearrow}(P_{1})$ and $W^{u}_{\swarrow}(P_{1})$ are connected with $(0,0)$. By continuity of the vector field in $(Q,C)$, there are parameter values for which the trajectory $W^{s}_{\swarrow}(P_{1})$ intersects the boundary of $\Gamma$ (defined in Lemma 3.1). Therefore, the stable manifold $W^{s}(P_{1})$ of $P_{1}$ act as a separatrix between the basin of attraction of the equilibrium point $P_{2}$ (when this is an attractor) and the equilibrium point $(0,0)$, see the top panels of Figures 4 and 5. We observe that by changing the parameter $Q$ the stable manifold $W^{s}_{\swarrow}(P_{1})$ moves down and the basin of attraction of $P_{2}$ shrinks (see light blue region of Figure 4), while the basin of attraction of $(0,0)$ increases (see orange region of Figure 4). Additionally, by the continuous dependence of the vector field on the parameter $Q$, there exist a $Q$ value for which $W^{u}_{\nearrow}(P_{1})$ intersects $W^{s}_{\swarrow}(P_{1})$ to form a homoclinic curve, see Figures 4 and 5. Upon further decreasing $Q$, and if the equilibrium point $P_{2}$ is an attractor and the $\omega$-limit of the unstable manifold $W^{u}_{\nearrow}(P_{1})$ is the equilibrium point $(0,0)$, the homoclinic curve breaks and there exists an unstable limit cycle which acts as a separatrix between the basins of attraction of these equilibrium points. Note that when the homoclinic curve breaks it generates a non-infinitesimal limit cycle (originating from a homoclinic bifurcation), which could coexist with the other limit cycle obtained via a Hopf bifurcation (infinitesimal limit cycle) when $P_{2}$ is a center-focus, see the upcoming Section 4.2 and Figure 9.

We summarise the above discussion in the following lemma:

Lemma 3.3

There exist conditions on the parameter values for which there is

1. (i)

a homoclinic curve determined by the stable and unstable manifold of equilibrium point $P_{1}=\left(u_{1},v_{1}\right)$; and 2. (ii)

a limit cycle that bifurcates from the homoclinic curve which surrounds the equilibrium point $P_{2}=\left(u_{2},v_{2}\right)$.

3.1.3 The collapse of the positive equilibrium points

Next, we study the case 2(c)ii in Section 2.1 where $C>M$, $\Sigma_{2}>0$, $\Delta=0$ and $N>M$ in equation (11). The equilibrium points $P_{1}$ and $P_{2}$ collapse such that $u_{1}=u_{2}=u_{E}=(2\Sigma_{2}+Q(M-N))/(2(C+NQ))$ and $v_{1}=v_{2}=v_{E}=(2N\Sigma_{2}-C(M-N))/(2N(C+NQ))$. That is system (3) has one equilibrium point of order two in the first quadrant given by $E=(u_{E},v_{E})$, see Figure 6.

Theorem 3.4

Let the system parameter be such that $C>M$, $\Sigma_{2}>0$, $\Delta=0$ and $N>M$ (case 2(c)ii of Section 2.1), then the stability of equilibrium point $E=(u_{E},v_{E})$ is as follows:

1. (i)

$E$* is a saddle-node repeller if $0<C<\dfrac{-A+\sqrt{16MN^{3}(M-N)^{2}+A^{2}})}{8N^{2}}$ ,* 2. (ii)

$E$* is a saddle-node attractor if $C>\dfrac{-A+\sqrt{16MN^{3}(M-N)^{2}+A^{2}})}{8N^{2}}$ ,*

where $A:=-8MN^{2}-(1-N)(M+N)^{2}$.

Proof. Since $u_{E}+v_{E}=(N-M)/(2N)$ and $\Delta=0$ we have $Q=\left((M-N)^{2}+4CN\right)/\left(4N(C-M)\right)$. Therefore, the Jacobian matrix of (3) at the equilibrium point $E$ becomes

[TABLE]

Thus, $\det{J(E)}=0$ and the trace of the Jacobian matrix (21) is given by

[TABLE]

with $\beta=N\left(M^{2}-6MN+N^{2}\right)-\left(M+N\right)^{2}$ and $\gamma=MN\left(M-N\right)^{2}$. The trace ${\rm{tr}}\left(J(E)\right)=0$ if and only if

[TABLE]

with $A$ as given in the statement of this theorem. For $0<C<C^{*}$ (23), we have that ${\rm{tr}}\left(J(E)\right)>0$ (22) and the equilibrium point $E$ is a saddle-node repeller. For $C>C^{*}$ (23), we have that ${\rm{tr}}\left(J(E)\right)<0$ (22) and the equilibrium point $E$ is a saddle-node attractor.

3.1.4 The collapse of $P_{1}$ with $(0,0)$

Next, we study the case when $C>M$, $\Sigma_{2}=0$ and $N>M$ such that $P_{1}$ merges with $(0,0)$. That is, we study case 2b of Section 2.1. Equation (11) has one positive solution and system (3) thus has one positive equilibrium point $P_{2}$ which is given by $P_{2}=(Q(N-M)/(C+NQ)$, $C(N-M)(N(C+NQ))$.

Theorem 3.5

If $C>M$, $\Sigma_{2}=0$ and $N>M$ (case 2b of Section 2.1), then the equilibrium point $P_{2}$ is a stable node.

Proof. Evaluating the determinant (20) at $P_{2}$ gives

[TABLE]

since $N>M$. The trace ${\rm{tr}}\left(J\left(u_{2},v_{2}\right)\right)$ (20) gives

[TABLE]

Finally, system (3) has no positive equilibrium points in the first quadrant in the remaining cases of Section 2.1, that is, in cases 2(c)i and 2iiiB. Therefore, the equilibrium point $(0,0)$ is a global attractor in these cases since $C>M$.

4 Bifurcation Analysis

In this section we discuss some of the possible bifurcation scenarios of system (3). The stability of $(1,0)$ and $P_{1}$ does not change for $\Delta>0$ (11), $Q>1$ and $C>M$, but whenever $\Delta=0$, $P_{1}$ undergoes a bifurcation as it collides with $P_{2}$. The stability of the equilibrium point, which is called $E$, at the bifurcation depends in a continuous fashion on the system parameters $C$ and $Q$, see Theorem 3.4.

Therefore, $C$ and $Q$ can be selected as natural candidates to act as bifurcation parameters.

4.1 Saddle-node Bifurcation

Theorem 4.1

Let the system parameters be such that $\Delta=0$ (11), $Q>1$, $C>M$, and

[TABLE]

then system (3) undergoes a saddle-node bifurcation at the equilibrium point $E$ (for changing $C$).

Proof. The proof of this theorem is based on Sotomayor’s Theorem [39]. For $\Delta=0$, there is only one equilibrium point $E=(u_{E},v_{E})$ in the first quadrant, with $u_{E}=(2C(1-Q)+Q(M+N))/(2(C+NQ))$ and $v_{E}=(CN(2Q-1)-2MNQ-CM)/(2N(C+NQ))$. From the proof of Theorem 3.4 we know that $\det\left(J(E)\right)=0$ if $\Delta=0$. Additionally, let $U=(1,1)^{T}$ be the eigenvector corresponding to the eigenvalue $\lambda=0$ of the Jacobian matrix $J(E)$, and let

[TABLE]

be the eigenvector corresponding to the eigenvalue $\lambda=0$ of the transposed Jacobian matrix $J(E)^{T}$.

If we represent (3) by its vector form

[TABLE]

then differentiating $F$ at $E$ with respect to the bifurcation parameter $C$ gives

[TABLE]

Therefore,

[TABLE]

Note that $W\cdot F_{C}(u_{E},v_{E};C)\neq 0$ under certain conditions on the parameters $(C,Q,N,M)$. Next, we analyse the expression $W\cdot[D^{2}F_{C}(u_{E},v_{E};C)(U,U)]$ and we first compute the quadratic form associated to the Hessian matrix

[TABLE]

At the equilibrium point $E$ and $U$, this becomes

[TABLE]

Therefore,

[TABLE]

and $W\cdot[D^{2}F(u_{E},v_{E};C)(U,U)]\neq 0$ by our assumptions on the parameters $(C,Q,N,M)$ (24). Therefore, by Sotomayor’s Theorem [39] it now follows that system (3) has a saddle-node bifurcation at the equilibrium point $E$ if the conditions on the parameters are met.

4.2 Hopf Bifurcation

Let $(U,V)$ be the coordinates of an equilibrium of (3) in the first quadrant. Then the set of equations $W(U,V)=0$ and $R(U,V)=0$ define implicitly a locally invertible transformation given by

[TABLE]

in $\Lambda:=\{(C,M,U,V)\in\mathbb{R}^{4}_{+}:\,CU-M(U+V)>0,\,\,U-1<0\}.$ Considering the time rescaling $\tau\mapsto\tau/(V(U+V))$, system (3) in parameter space $\Lambda$ has the form:

[TABLE]

where, for convenience, we again use the notation $(x,y)$ to name the state variables. System (25) is $C^{\infty}$-equivalent to (3) in parameter space $\Lambda$. Moreover, the (positive) equilibrium coordinates appear now as the explicit parameters $(U,V)$. In what follows, we will derive conditions such that (25) undergoes a generic Hopf bifurcation at $(U,V)$, see [40, 41] for more details.

The Jacobian matrix of (25) at $(U,V)$ is

[TABLE]

The trace and determinant of $J(U,V)$ are given by

[TABLE]

respectively, where

[TABLE]

Whenever $T=0$ and $D>0$, the eigenvalues of $J(U,V)$ are purely imaginary and non-trivial. Moreover, since we have

[TABLE]

the Hopf bifurcation in (25) is generically unfolded by parameter $M$ [41]. In particular, it follows from (28) that equation $T(C,M,U,V)=0$ implicitly defines the function

[TABLE]

We now calculate the first Lyapunov quantity [40, 41] in order to determine genericity conditions. We follow the derivation in [40] and move the equilibrium point $(U,V)$ of the system (25) to the origin via the translation $x\mapsto x+U$, $y\mapsto y+V$ to obtain the equivalent system

[TABLE]

In particular, the Jacobian matrix of (30) at the equilibrium point $(0,0)$ coincides with $J(U,V)$ in (26). Substitution of (29) into $J(U,V)$ gives

[TABLE]

with ${\rm{tr}}\left(J_{H}\left(U,V\right)\right)\equiv 0$ and $\det\left(J_{H}(U,V)\right)=UV^{2}(U+V)^{2}D_{H}$ (27) where

[TABLE]

If $D_{H}>0$, then \mathbf{v}_{1}=\left(\begin{array}[]{c}\dfrac{U(U-2U^{2}+V-3UV-V^{2})}{CV^{2}}\\ 1\end{array}\right) and \mathbf{v}_{2}=\left(\begin{array}[]{c}-w\\ 0\end{array}\right) are the generalised eigenvectors of $J_{H}(U,V)$, where

[TABLE]

The change of coordinates \left(\begin{array}[]{c}x\\ y\end{array}\right)\mapsto[\mathbf{v}_{1}\,\mathbf{v}_{2}]\left(\begin{array}[]{c}x\\ y\end{array}\right) allows us to express system (30) with $T(C,M,U,V)=0$ in the form

[TABLE]

where

[TABLE]

[TABLE]

and G=G\big{(}U,V\big{)}=U-2U^{2}+V-3UV-V^{2}.

System (32) and equations (31), (33) and (34) allow us to use the derivation in [40] for the direct calculation of the first Lyapunov quantity $L_{1}$. In this way we obtain the following expression:

[TABLE]

where

[TABLE]

Thus we have obtained the following result.

Theorem 4.2

Let $(C,M,U,V)\in\Lambda$ be such that $T(C,M,U,V)=0$, $D_{H}>0$ and $l_{1}\neq 0$ (35). Then (25) undergoes a codimension-one Hopf bifurcation at the equilibrium point $(U,V)$. In particular, if $l_{1}<0$ (resp. $l_{1}>0$), the Hopf bifurcation is supercritical (respectively subcritical), and a stable (respespectively unstable) limit cycle bifurcates from $(U,V)$ under suitable parameter variation.

Figure 7 shows the relevant bifurcation sets in the $(U,V)$-plane for increasing values of $C$. The black curve defined by $T(C,M,U,V)=0$ corresponds to a Hopf bifurcation only inside the shaded region labelled as $\Lambda_{h}=\Lambda\cap\{D_{H}>0\}$. The Hopf bifurcation curve is separated into two segments: the upper one corresponds to $l_{1}<0$ (supercritical), and the lower one to $l_{1}>0$ (subcritical). The division occurs at the point $B=L_{1}^{-1}(0)\cap T^{-1}(0)\cap\Lambda_{H}$ –also known as Bautin bifurcation– where the Hopf bifurcation curve intersects the level set $\{L_{1}=l_{1}=0\}$; at such point the Hopf bifurcation of (25) at $(U,V)$ is degenerate. The actual codimension of this singularity – and the stability of further limit cycles that bifurcate – is determined by the sign of the so-called second Lyapunov quantity $L_{2}$ [41]. In order to check the sign of $L_{2}$ for representative parameter values, we performed a bifurcation analysis with Auto-07p [42]. Figure 8 shows the corresponding bifurcation diagrams in the $(N,Q)$-plane in a neighbourhood of the Bautin point $B$ when $C=0.363$ and $C=0.6$ in rows (a) and (b), respectively, and $M=0.16$ remains fixed. The right column shows the same curves near $B$ after the change of parameter $\widetilde{Q}=100(Q-H(N))$, where $Q=H(N)$ defines (locally) the curve of Hopf bifurcation in the $(N,Q)$-plane. The advantage of this rescaling is that the curves in the $(N,\widetilde{Q})$-plane are now well separated. In Figure 8, a curve of saddle-node (or limit point) bifurcations of limit cycles ($LPC$) emanates from $B$ and eventually moves towards a curve of heteroclinic bifurcation (labeled $het$). At the curve $het$, a heteroclinic cycle is formed which joins $(0,0)$ with $(1,0)$ along the horizontal axis, and then goes back to the origin in the interior of the phase space. Moreover, for parameter values in the small region bounded by the curves $LPC$ and $het$, we have two concentric limit cycles surrounding the point $P_{2}$ in phase space. The continuation procedure confirms that for $C=0.363$ we have $L_{2}>0$ and the largest limit cycle is unstable (similar to Figure 4); and for $C=0.6$, we have $L_{2}<0$, so the largest limit cycle is stable (similar to Figure 9). The theory of generalised Hopf bifurcations is only local, meaning that it is valid for parameter values in a neighbourhood of the special Bautin point $B$ [40, 41]. Hence the two phase portraits in Figures 4 and 9 do not contradict the theory nor one another: they can be found well away from each other in parameter space. We conjecture that there is a higher codimension bifurcation point organising the overall dynamics in our model (which is yet to be found), and the two phase portraits are just two possible “perturbations” away from that point.

4.3 Bogdanov-Takens bifurcation

We now give conditions such that our model undergoes a Bogdanov-Takens bifurcation under suitable parameter variation at the equilibrium point $E$. We refer to [41] and the references therein for the derivation of the genericity and transversality conditions that need to be verified during this proof.

For the sake of clarity, it is convenient to state the dependence of the vector field (3) on parameters $C$ and $Q$ explicitly. Hence, throughout this section we denote

[TABLE]

where we use notation $(u,v)$ for the state variables. Also, let us denote the Jacobian matrix of $X$ with respect to the variables $(u,v)$ as

[TABLE]

Note that this notation for a Jacobian is standard in bifurcation theory as it explicitly highlights the dependence on all the dependent and independent variables [41].

Step 1. We verify that the singularity has a double zero eigenvalue with geometric multiplicity one.

From the proof of Theorem 3.4, conditions $\Delta=0$ in (11) and $tr(J(E))=0$ in (22) determine a bifurcation point $(C,Q)=(C^{*},Q^{*})$ implicitly in the form:

[TABLE]

where $A=-8MN^{2}-(1-N)(M+N)^{2}.$ Moreover, at $(C,Q)=(C^{*},Q^{*})$, the equilibrium point $E=(u_{E},v_{E})$ has a Jacobian matrix given by

[TABLE]

with $J_{12}=M^{2}-2MN+N(4C^{*}+N)$ and a double zero eigenvalue. However, note that (36) is not the null matrix. The corresponding generalised eigenvectors of (36) are given by

[TABLE]

where

[TABLE]

and

[TABLE]

It follows that (36) is nilpotent and that the double zero eigenvalue has geometric multiplicity one.

Step 2. The next goal is to state the transversality condition of a Bogdanov-Takens bifurcation, namely, that the map

[TABLE]

is regular at $(u,v,C,Q)=(u_{E},v_{E},C^{\ast},Q^{\ast})$, where $\mathcal{T}$ and $\mathcal{D}$ are the trace and determinant of the Jacobian matrix $\partial X/\partial(u,v)\left(u,v;C,Q\right)$, respectively.

After some calculations, the determinant of the $4\times 4$ Jacobian matrix $D\Psi$ of the map $\Psi$ can be written as

[TABLE]

with:

[TABLE]

In particular, straightforward substitution and algebraic simplification leads to $F(u_{E},v_{E},C^{\ast},Q^{\ast})=G_{1}$, where

[TABLE]

and $Z=(M+N)\sqrt{M^{2}(N-1)^{2}+(N-1)^{2}N^{2}+2MN(1+6N+N^{2})}.$

It follows that if

[TABLE]

then we have $\det D\Psi(u_{E},v_{E},C^{\ast},Q^{\ast})\neq 0,$ which ensures that the map $\Psi$ is regular at $(u,v,C,Q)=(u_{E},v_{E},C^{\ast},Q^{\ast})$.

Step 3. We now construct a change of coordinates —and give sufficient conditions to do so— which transforms $X(u,v;C,Q)$ into a normal form of the Bogdanov-Takens bifurcation; we refer to [41] again.

Let us first define the following auxiliary expressions:

[TABLE]

where $z_{1}$ and $z_{2}$ are given by (38) and (39), respectively.

Let us now move the equilibrium point $(u_{E},v_{E})$ of (3) to the origin via the translation $u\mapsto u+u_{E}$, $v\mapsto v+v_{E}$ to obtain the equivalent system

[TABLE]

In particular, the Jacobian matrix of (45) at the equilibrium point $(0,0)$ at the bifurcation point $(C^{*},Q^{*})$ coincides with $\partial{X}/\partial(u,v)\left(u_{E},v_{E};C^{*},Q^{*}\right)$ in  (36).

Let $\mathbf{P}=[\mathbf{v}_{1},\mathbf{v}_{2}]$ be the matrix whose columns are $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$; see (37)–(39). Next, consider the following change of coordinates:

[TABLE]

Then, the vector field given by

[TABLE]

is $\mathcal{C}^{\infty}$-conjugated to system (45).

Taking a Taylor expansion of $\mathbf{J}(x,y;C,Q)$ with respect to $(x,y)$ around $(x,y)=(0,0)$ and evaluating at $(C,Q)=(C^{\ast},Q^{\ast})$, one obtains

[TABLE]

provided $z_{1}z_{2}-1\neq 0$, where

[TABLE]

If $b_{20}\neq 0$ and $a_{20}+b_{11}\neq 0$, then the theory of normal forms for bifurcations [41] ensures that our system fulfils the necessary genericity conditions to undergo a codimension two Bogdanov-Takens bifurcation. In particular, since $z_{1}\neq 0$ in (38), condition

[TABLE]

ensures that $b_{20}\neq 0$. Furthermore, after some algebraic manipulation one obtains $a_{20}+b_{11}=G_{3}G_{4}.$ Hence, condition $a_{20}+b_{11}\neq 0$ is equivalent to

[TABLE]

In summary, step 1 and inequalities (41), (46) and (47) ensure that the genericity and transversality conditions of a codimension two Bogdanov-Takens normal form are satisfied. Hence, there exists a smooth, invertible transformation of coordinates, an orientation-preserving time rescaling, and a reparametrisation such that, in a sufficiently small neighbourhood of $(u,v,C,Q)=(u_{E},v_{E},C^{\ast},Q^{\ast})$, the system (3) is topologically equivalent to one of the following normal forms of a Bogdanov-Takens bifurcation:

[TABLE]

where the sign of the term $\xi_{1}\xi_{2}$ in (48) is determined by the sign of $b_{20}(a_{20}+b_{11})$.

We are now in a position to state the corresponding result.

Theorem 4.3

Let the system parameter be such that $\Delta=0$ (11), $tr(J(E))=0$ (22) and $G_{1,2,3,4}\neq 0$ (40), (42) – (44), then the equilibrium point $E=(u_{E},v_{E})$ in (3) undergoes a codimension-two Bogdanov-Takens bifurcation.

4.4 Bifurcation diagram

In order to obtain the bifurcation diagram of system (3) for $Q>1$, $C>M$ and the parameters $M$ and $N$ fixed we follow [28, 29] and use the numerical bifurcation package MATCONT [43]222Note that the Matlab package ode45 was used to generate the data for the simulations and then the PGF package (or tikz) was used to generate the graphics format.. The bifurcation curves obtained from Sections  4.1 – 4.3 divide the $(Q,C)$-parameter space into four parts, see Figure 10. We observe that system (3) has only one positive equilibrium point if $(Q,C)$ are located on the saddle-node curve, see Theorem 3.4, while system (3) does not have any equilibrium points in the first quadrant, and therefore $(0,0)$ is global attractor, if $(Q,C)$ are located above the saddle-node curve (green region of Figure 10).

System (3) has two positive equilibrium points, namely $P_{1}=\left(u_{1},v_{1}\right)$ and $P_{2}=\left(u_{2},v_{2}\right)$ with $0<u_{1}\leq u_{2}$ and $0<v_{1}\leq v_{2}$, if $(Q,C)$ are located below the saddle-node curve (blue, yellow and brown regions of Figure 10). In these regions the equilibrium point $P_{1}$ is always a saddle point. The equilibrium point $P_{2}$ is stable when $(Q,C)$ are located below the Hopf curve (yellow and brown regions of Figure 10) and it is surrounded by an unstable limit cycle when $(Q,C)$ is, in addition, above the homoclinic curve (yellow region of Figure 10). The equilibrium point $P_{2}$ is unstable when $(Q,C)$ are located between the Hopf curve and the saddle-node curve (blue region Figure 10). See Figure 11 for typical phase planes each of the four regions.

5 Conclusions

In this manuscript, a Bazykin predator-prey model with predator intraspecific interactions and the ratio-dependent functional response was studied. Using a diffeomorphism we transformed the Bazykin predator-prey model (1) to a topologically equivalent system (3) and subsequently analysed this nondimensionalised system.

This system has four system parameters which determine the number of the equilibrium points and their stability. We showed that the equilibrium point $P_{1}$ is always a saddle point and $(1,0)$, which corresponds to the rescaled carrying capacity, is a saddle point if $C>M$ and a stable node if $C<M$. We showed in Theorem 3.1 that the origin has complex dynamics and by using vertical and horizontal blow-ups we determined the dynamic in the neighbourhood of the origin. Furthermore, for some sets of parameters values the stable manifold of $P_{1}$ determines a separatrix curve which divides the basins of attraction of $(0,0)$ and $P_{2}$. As a result, the equilibrium point $P_{2}$ can be stable, stable surrounded by unstable limit cycle or unstable, depending on the trace of its Jacobian matrix, see Theorem 3.3. Additionally, we proved that system (3) undergoes a codimension-one Hopf bifurcation at the equilibrium $P_{2}$ and we also provided evidence for a degenerate codimension-two Hopf bifurcation at the Bautin point, see Theorem 4.2. Moreover, the equilibrium points $P_{1}$ and $P_{2}$ collapse for $\Delta=0$ (11), i.e. $P_{1}=P_{2}=E$, and system (3) undergoes to a saddle-node bifurcation [39], see Theorem 4.1, or a Bogdanov-Takens bifurcation, see Theorem 4.3.

Since the function $\varphi$ (2) is a diffeomorphism preserving the orientation of time, the dynamics of system (3) is topologically equivalent to the dynamics of Bazykin predator-prey model (1). Therefore, we can conclude that there exists self-regulation in system (1) for certain population sizes, that is, the species can coexist. However, system (3) is sensitive to changes in the parameters and also disturbances of the population size. We can see this impact in the size of the basins of attraction of the equilibrium points $(0,0)$ and $P_{2}$ in Figures 4 and 11. In these figures, the orange region represents population sizes that lead to the extinction of both populations and the grey region represents population sizes that lead to the stabilisation of both populations over time. We showed that the stabilisation of the predator and the prey populations depend on the values of the parameters $Q$ and $C$ by taking the parameters $M$ and $N$ fixed333Note that the parameter $Q$ correspond to the rescaled per capita predation rate $q$ and the parameter $C$ correspond to the rescaled efficiency with which predators convert consumed prey into new predators $c$.. For example, small predation rate $q$, or small efficiency with which predators convert consumed prey into new predators $c$, increases the area of coexistence which is related to basins of attraction of $P_{2}$, see Figure 11. Moreover, we show in Figure 11 that the stabilisation or extinction of both the predator and the prey populations depends also on the species initial density. For example, for low predation rate the stable manifold of the equilibrium point $P_{1}$ act as a separatrix and initial conditions with large predator density lead to the extinction of both populations, see bottom panel of Figure 11. In contrast, if the parameters $(Q,C)$ are located in the yellow region of the bifurcation diagram, see Figures 10 and 11, the stabilisation of both populations are bounded by the unstable limit cycle. In this case, the stabilisation or extinction depends on the initial conditions as is showed in Figure 11. If the initial condition is in the region interior to the limit cycle then the predator and prey population will coexist. However, if the initial condition is in the region outside of the limit cycle then both populations become extinct.

Finally, we observe that intraspecific interaction occurs in the predator population. Also, in most cases considered in this manuscript either both species go extinct or a coexistence state will be reached. For instance, in the left panel of Figure 2 we observe the extinction of both populationswe observe the extinction of both populations (while in the right panel only the predator population goes extinct). The dynamics observed in the left panel can be connected to the intraspecific interaction in the predator population. As observed in most of our results presented, this type of interaction improves the species’ adaption (see Figure 4). For example, when the prey is abundant, which may limit the predator population density, the intraspecies competition among predators becomes stronger. Hence, the predator density is insufficient to control the prey population and thus the prey is more likely to escape the predator [4, 21].

Acknowledgment

Pablo Aguirre was funded by Proyecto Interno UTFSM PI-LI-19-06. Pablo Aguirre also thanks Proyecto Basal CMM Universidad de Chile.