Numerical approach for approximating the Caputo fractional-order derivative operator

1.   Introduction

The advance of ecological theory has had an increasing growth due to the use of mathematical models of living systems. In studying ecological models, it is frequent to make a mathematical analysis of simple theoretical models and later a numerical analysis of more realistic and complex models.

Although analytical approaches are exact and general, they can be limited to simple models. However, the formal analysis of these simple models, described by bidimensional nonlinear differential equation systems (ODES) provides the skeleton for most complex theories for the species interplays and it is the main source of specific hypotheses to be tested by numerical simulations.

But, it is well-known that many types of modelling have been formulated for this important interaction, from the seminal Lotka-Volterra model ^[1] in 1925, a model obeying mass-action principle ^[2]. The knowledge detailed of the dynamics of predator-prey models described by ODES is an interesting and basic issue for understanding the behavior of complex food chains; also for the application of other related mathematical tools, as are the delay, impulsive or stochastic differential equations.

One of this type of model was proposed by the English ecologist Patrick H. Leslie in 1948 ^[3], assuming that the equation for predators is a logistic-type growth function, in which the conventional environmental carrying capacity for predators $K_{y}$ is proportional to the prey population size $x = x\left(t\right)$, that is $K_{y} = K(x) = nx$ ^[1,4,5].

Clearly, the model proposed by Leslie does not fit to the Lotka-Volterra scheme ^[1]. It has been strongly criticized by presenting anomalies in their predictions since it vaticinates that even in very low prey population size, when the consumption rate per predator is almost zero, predator population might increase, if the predator/prey ratio is very small ^[1].

Nonetheless, in the case of critical scarcity, some predator species can switch over to other available food. This ability can be modeled by adding a positive constant $k_{2}$ in the carrying capacity $K_{y}(x)$, being described now by $K(x) = nx+k_{2}$, a function of the available prey population ^[1,5].

In that situation is said that the model is represented by a Leslie-Gower scheme and it is also known as modified Leslie-Gower model ^{[6,7,8,9,10,11,12]}; if $x = 0$, then $K(0) = k_{2}$, concluding that the predator is generalist since it searches an alternative food in absent of its favorite prey.

In this work, we describe in detail the dynamics of the model proposed in ^[8], in which was assumed the predator is a generalist and the action of the predators was represented by the hyperbolic functional response. Then, the modified Leslie-Gower here presented can enhance the anomalies of the original model formulated by Leslie ^[3,4].

In spite of this model has been used in several recent papers in other studies ^[13,14,15], using delay, impulsive or stochastic differential equations or incorporating ecological phenomena as Allee effect ^[9,16], or prey refuge use ^[17,18], or it has been used to study an optimal control problem ^[19] but its dynamic yet is not being completely analyzed.

On the other hand, the predator functional response or consumption function refers to the change in attacked prey density per unit of time per predator when the prey density changes ^[5].

Here, we will consider that the predator functional response is depending only on the prey population and expressed by the function $h(x) = \ \frac{q\ x }{x\ +\ k_{1}}$, the hyperbolic functional response ^[20,21]. The parameter $a$ is an abruptness measure of the function ^[22]. If $a\rightarrow 0$, the curve grows quickly, while if $a\rightarrow K$, the curve grows slowly, that is, a bigger amount of prey is necessary to obtain $\frac{q}{2}$. This parameter indicates the quality of the prey as food for the predator.

When $K(0) = k_{2} = 0$, i.e., $K_{y}$ is proportional to prey abundance, and the predator consumption rate is the same hyperbolic functional response, the model is named as May-Holling-Tanner model ^{[1,5,23,24,25,26]}. Then, our analysis also extends the results obtained for that model.

We will describe the behavior of the model using the bifurcation theory ^[20], depending on the parameter values and to classify the different dynamics resulting. Also, we compare the results with those obtained for the May-Holling-Tanner model ^[25] and the model analyzed in ^[7] considering the Allee effect in prey ^[20].

Thus, our outcomes extend and complement the obtained results in ^[8] and ^[27], since we established conditions on the parameter for the existence of two, one or none positive equilibrium points.

Assuming the existence of a unique equilibrium at the interior of the first quadrant we demonstrate the existence of: i) two limit cycles encircling a stable positive equilibrium point; ii) a separatrix curve in the phase plane dividing the behavior of the trajectories, and iii) a homoclinic curve.

We also conclude that local stability of an equilibrium point does not imply the global stability, such as it was also established in ^[25] for the May-Holling-Tanner model.

The case in which two positive equilibria exist in the first quadrant is analyzed in citegonzalez5.

The rest of the paper is organized as follows: In the next section, the modified Leslie-Gower model is presented; in Section 3, the main properties of the model are proved. To reinforce these results, some numerical simulations are shown in Section 3, meanwhile, in Section 4 we discuss the obtained outstanding, explaining the ecological meanings of some of them.

2.   The model

The modified Leslie-Gower type model considering a generalist predator proposed in ^[8], it is described by the following ordinary differential equation system of Kolmogorov type ^[28]

with $x\left(0\right) \geq 0$ and $y\left(0\right) \geq 0$, where $x = x\left(t\right)$ and $y = y\left(t\right)$ indicate the prey and predator population sizes respectively for $t\geq 0$, measured as number of individuals, density or biomass. The parameters are all positives, i.e., $\mu = \left(r_{1}, b_{1}, a_{1}, k_{1}, r_{2}, a_{2}, k_{2}\right) \in \mathbb{R} {}{}_{+}^{7}$.

As system $\left(2.1\right)$ is Kolmogorov type ^[28], the coordinates axis are invariant sets and the model is defined at

   $\Omega = \left\{ \left(x, y\right) \in \mathbb{R}^{2}/\ x\geq 0, \ y\geq 0\right\} = \mathbb{R}_{0}^{+}\times \mathbb{R}_{0}^{+}$.

The equilibrium points of system $\left(2.1\right)$ or vector field $X_{\mu }$ are $\left(0, 0\right)$, $\left(\frac{r_{1}}{b_{1}}, 0\right)$, $\left(0, \frac{r_{2}k_{2}}{a_{2}}\right)$ and $\left(x_{e}, y_{e}\right)$ satisfying the equation of the isoclines $y = \frac{r_{2}}{a_{2}} \left(x+k_{2}\right)$ and $y = \frac{1}{a_{1}}\left(r_{1}-b_{1}x\right) \left(x+k_{1}\right)$. Clearly, $\left(x_{e}, y_{e}\right)$ can be a positive equilibrium point or cannot exists there, depending of the sign of the factor $r_{1}-b_{1}x$.

We note that system $\left(2.1\right)$ can rewrite as

with $K = \frac{r_{1}}{b_{1}}$, $n = \frac{r_{2}}{a_{2}}$ and $c = \frac{r_{2}k_{2} }{a_{2}}$, where $\eta = \left(r_{1}, K, a_{1}, k_{1}, r_{2}, b, c\right) \in \mathbb{R}{}{}_{+}^{7}$ and the parameters have the following meanings

$r_{1}$ is the intrinsic prey growth rate,

$K$ is the prey environmental carrying capacity,

$a_{1}$ is the consuming maximum rate per capita of the predators (satiation rate),

$k_{1}$ is the amount of prey to reach one-half of $q$ (that is, it is half saturation rate),

$r_{2}$ is intrinsic predator growth rate,

$n$ is the food quality and it indicates how the predators turn eaten prey into new predator births,

$c$ is the amount of alternative food available for the predators.

By ecological reason $k_{1} < K$. This parameter $c$ indicates that the predator is a generalist since it does not exist available prey (its favorite food), it searches an alternative resource. If $c = 0$, the system describes the May-Holling-Tanner model ^[23,25,29], which is not defined in $x = 0$.

To simplify the calculations we reduce system $\left(2.2 \right)$ to a normal form; following the methodology used in ^[30,31], which involved a change of variable and a time rescaling ^[32] given by the function $\varphi :\bar{\Omega}\times \mathbb{R} \longrightarrow \Omega \times \mathbb{R}$, so that

It is easy to see $\det D\varphi \left(u, v, \tau \right) \; = \frac{Kn\left(u\ +\ \frac{k_{1}}{K}\right) \left(u+\frac{c}{Kn}\right) }{r_{1}} > 0$ .

Thus, $\varphi$ is a diffeomorphism ^[32,33] preserving the time orientation; for this reason the vector field $X_{\mu }$ is topologically equivalent to the vector field $Y_{\sigma} = \varphi \circ X_{\mu }$ with $Y_{\sigma } = \; P(u, v)\ \frac{\partial }{\partial u}+Q(u, v)\ \frac{\partial }{\partial v}\;$and the associated differential equation system is given by the polynomial system of fourth degree, being also of Kolmogorov type ^[28]:

where $\sigma = (A, S, C, Q)\in ]0, 1[\times \mathbb{R}_{+}^{3}$ with, $S = \frac{ r_{2}}{r_{1}}$, $Q = \frac{na_{1}}{r_{1}}$, $C = \frac{c}{Kn}$ and $A = \frac{k_{1} }{K} < 1$ by ecological reasons.

System $\left(2.3\right)$ or vector field $Y_{\lambda }$ is defined in

  $\bar{\Omega}$ $= \{(u, v)\in \mathbb{R}^{2}/\ u\geq 0, \ v\geq 0\}$.

2.1.   Determination of the equilibria

The equilibrium points, critical points or singularities of system $\left(2.3\right)$ are $(1, 0)$, $(0, 0)$, $(0, C)$ and the points $\left(u_{e}, v_{e}\right)$ which lie in the intersection of the isoclines

  $v = u+C$ and $v = \frac{1}{Q}\left(1-u\right)(u\ +\ A)$.

Then, the abscissa $u$ of the positive equilibrium points is a solution of the second degree equation:

According to the Descartes' Rule of Signs and to the sign of

   $1-A-Q$, $CQ-A$, and $\Delta = \left(1-A-Q\right) ^{2}-4\left(CQ-A\right)$

equation $\left(2.4\right)$ can have two, one or none positive roots.

1) If $1-A-Q > 0$ and $CQ-A > 0$, the solutions of equation $\left(2.4 \right)$ are:

   $u_{1} = \frac{1}{2}\left(1-A-Q-\sqrt{\Delta }\right)$ and $u_{2} = \frac{1}{2}\left(1-A-Q+\sqrt{\Delta }\right)$.$\qquad$

Then, according to the $\Delta$ sign, we have three possibilities:

1.1 There are two positive equilibrium points $P_{1} = \left(u_{1}, u_{1}+C\right)$ and $P_{2} = \left(u_{2}, u_{2}+C\right)$ with $0 < u_{1} < u_{2}$, if and only if, $C < \frac{1}{4Q}\left(4A+\left(1-A-Q\right) ^{2}\right)$.

1.2 There is a unique equilibrium point at interior of the first quadrant, if and only if, $C = \frac{1}{4Q}\left(4A+\left(1-A-Q\right) ^{2}\right)$, having $\left(u_{1}, u_{1}+C\right) = \left(u_{2}, u_{2}+C\right) = \left(E, E+C\right)$ with $E = \frac{1-A-Q}{2}$.

1.3 There are not equilibrium points at interior of the first quadrant, if and only if, $C > \frac{1}{4Q}\left(4A+\left(1-A-Q\right) ^{2}\right)$.

2) The solutions of equation $\left(3\right)$ are $u_{1} < 0 < u_{2}$, if and only if,

2.1$\qquad 1-A-Q > 0$ and $CQ-A < 0$, or

2.2$\qquad 1-A-Q < 0$ and $CQ-A < 0$.

In both cases, the unique positive equilibrium point is

   $P_{2} = \left(u_{2}, u_{2}+C\right) = \left(L, L+C\right)$, with $L = \frac{1}{2}\left(1-A-Q+\sqrt{\Delta }\right)$.

Note that if $CQ-A < 0$, in term of the original parameters it has

$\qquad k_{2}\frac{a_{1}r_{2}}{a_{2}r_{1}}-k_{1} < 0$, i.e., $k_{2}\frac{ a_{1}r_{2}}{a_{2}r_{1}} < k_{1}$,

then, $\frac{k_{2}r_{2}}{a_{2}} < \frac{k_{1}r_{1}}{a_{1}}$, obtaining the unique case studied in ^[8].

Therefore, the obtained results on ^[8] correspond to these particular cases.

3) If $1-A-Q = 0$ and $CQ-A < 0$, equation $\left(2.4\right)$ has two solutions, one positive and other negative. As $1-A = Q$, the unique equilibrium point at interior of the first quadrant is

   $P_{2} = \left(F, F+C\right)$ with $F = \sqrt{\Delta }$ and $\Delta = A-C\left(1-A\right) > 0$. So, $C < \frac{A}{1-A}$.

4) If $1-A-Q > 0$ and $CQ-A = 0$, equation $\left(2.4\right)$ has two solutions

$\qquad u_{1} = 0$ and $u_{2} = \ G = 1-A-Q = 1-A-\frac{A}{C}$.

Clearly, the point $P_{1}$ coincides with $\left(0, C\right)$. Then,

   $P_{2} = \left(\frac{C-A-AC}{C}, \frac{\left(C-A\right) \left(C+1\right) }{C}\right)$, with $C-A-AC > 0$ and $C-A > 0$.

is the unique positive equilibrium point.

5) There are not equilibrium points at the interior of the first quadrant, if and only if, equation $\left(2.4\right)$ does not have real roots, i.e., if and only if,

5.1   $1-A-Q = 0$ and $CQ-A > 0$, or

5.2    $1-A-Q < 0$ and $CQ-A\geq 0$.

The above classification 1-5 implies the study of different cases in this class of systems, depending on the relations between the parameters $A$, $C$, and $Q$. Furthermore, system $\left(2.1\right)$ has a meaningful difference with May-Holling-Tanner model ^[25], respect to the number of equilibrium points, since it can have up to two positive equilibrium points, besides the new singularity $\left(0, C\right)$.

In Figure 1, we show the relative position among the prey and predator nullclines.

3.   Main results

System $\left(2.3\right)$ or vector field $Y_{\sigma }$ has the following properties:

Lemma 1. Existence of invariant region

The set

   $\bar{\Gamma} = \left\{ \left(u, v\right) \in \bar{\Omega}/\ 0\leq u\leq 1, \; v\geq 0\right\}$

is an invariant region.

Proof. Clearly, the $u-axis$ and the $v-axis$ are invariant sets because the system is a Kolmogorov type.

If $u = 1$, we have

   $\frac{du}{d\tau } = -Qv\left(1+C\right) < 0$

and whatever it is the sign of

   $\frac{dv}{d\tau } = S\left(1+A\right) \left(\ 1+C-\ v\right) \ v$

the trajectories enter and remain in the region $\bar{\Gamma}$.

By the diffeomorphism $\varphi$, the set

   $\Gamma = \left\{ \left(x, y\right) \in\Omega /\ 0\leq x\leq K, \; y\geq 0\right\}$

is an invariant region of the system $\left(2.1\right)$.

Lemma 2. Boundedness of the solutions

The solutions are bounded.

Proof. To determinate that the solutions are bounded we will use the Poincaré compactification ^[35]:

We change the variable and the time rescaling given by the function $\theta :\tilde{\Omega}\times \mathbb{R}\longrightarrow \breve{\Omega}\times \mathbb{R}$, so that:

   $\theta \left(X, Y, T\right) = \left(\frac{X}{Y}, \frac{1}{Y}, Y^{3}T\right) = (u, v, \tau)$

Doing a large algebraic work we get the system:

   $\bar{U}_{\eta }:\left\{ \begin{array}{lcl} \frac{dX}{dT} & = & -X\left(X^{3}+\left(AY-Y+CY+SY\right) X^{2}+a_{1}X+a_{2}\right) \\ \frac{dY}{dT} & = & -SY^{2}\left( X^{2}+\left(AY+CY-1\right) X+AY\left(CY-1\right) \right) \end{array} \right.$

with:

   $a_{1} = MY^{2}-CY^{2}-AY^{2}+SY^{2}+ACY^{2}-AMY^{2}-CMY^{2}$

   $a_{2} = QY^{2}-SY^{2}-ACY^{3}+AMY^{3}+CMY^{3}+ASY^{3}+CSY^{3}-ACMY^{3}$

The Jacobian matrix ^[23] in the new system is

$\qquad D\bar{U}_{\eta }(X, Y) = \left(\begin{array}{cc} -\left( 4X^{3}+b_{1}X^{2}+b_{2} X+b_{3}\right) & -DJ\left(1, 2\right) \\ -SY^{2}\left(2X+AY+CY-1\right) & -DJ\left(2, 2\right) \end{array} \right)$

with:

$\qquad DJ\left(1, 2\right) = X \left(\left(A+C+S-1\right) X^{2}+b_{4} X+b_{5}\right)$

$\qquad DJ\left(2, 2\right) = SY\left(2X^{2}-3AY-2X+4ACY^{2}+3AXY+3CXY\right)$

and

$\qquad b_{1} = 3AY-3Y+3CY+3SY$

$\qquad b_{2} = 2QY-2CY^{2}-2AY^{2}-2SY+2ACY^{2}+2ASY^{2}+2CSY^{2}$

$\qquad b_{3} = CQY^{2}-ASY^{2}-ACY^{3}+ACSY^{3}$

$\qquad b_{4} = Q-S-2AY-2CY+2ACY+2ASY+2CSY$

$\qquad b_{5} = 2CQY-2ASY-3ACY^{2}+3ACSY^{2}$

Evaluating the matrix $D\bar{U}_{\eta }(X, Y)$ in the point $(0, 0)$ we obtain

$\qquad D\bar{U}_{\eta }\left(0, 0\right) = \left(\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array} \right)$

To desingularize the origin in the vector field $\bar{U}_{\eta }$ we apply the blowing up method ^[33].

Let $X = rw$ and $Y = w$, replacing and time rescaling $\zeta = w^{2}T$. After an algebraic work we obtain the system:

$\qquad \breve{U}_{\eta }(r, w):\left\{ \begin{array}{lcl} \frac{dr}{d\zeta } & = & r\left(-CQ-Qr+ACw-ACrw\right) \\ \frac{dw}{d\zeta } & = & Sw\left(A+r-rw-ACw-Arw-Crw\right) \end{array} \right.$

The Jacobian matrix of the system which is:

   $D\breve{U}_{\eta }(r, w) = \left(\begin{array}{cc} ACw-2Qr-CQ-2ACrw & -ACr\left(r-1\right) \\ -Sw\left(w+Aw+Cw-1\right) & \breve{U}\left(2, 2\right) \end{array} \right)$

with $\breve{U}\left(2, 2\right) = -S\left(2rw-r-A+2ACw+2Arw+2Crw\right)$

Evaluating the matrix $D\breve{U}_{\eta }(r, w)$ in the point $(0, 0)$ we obtain:

$\qquad D\breve{U}_{\eta }(0, 0) = \left(\begin{array}{cc} -CQ & 0 \\ 0 & AS \end{array} \right)$.

Clearly.

$\qquad \det D\breve{U}_{\eta }(0, 0) = -CAQS < 0$

Therefore, we have that $(0, 0)$ is a saddle point of the vector field $\breve{U}$ and of $\bar{U}$, then the point $(0, \infty)$ is a saddle point in the compactified vector field of $Y_{\sigma }\left(u, v\right)$.

Then, the trajectories are bounded.

To determine the nature local of the hyperbolic equilibrium points we must obtain the Jacobian matrix of system $\left(2.3\right)$, which is:

$\qquad DY_{\sigma }\left(u, v\right) = \left(\begin{array}{cc} a_{11}\left(u, v\right) & -Qu\left(u+C\right) \\ Sv\left(A+C+2u-v\right) & S\left(A+u\right) \left(C+u-2v\right) \end{array} \right)$

with

$a_{11}\left(u, v\right) =$ $-4u^{3}+3\left(1-C-A\right) u^{2}+2\left(\left(A+C\left(1-A\right) -Qv\right) \right) u+C\left(A-Qv\right)$,

Theorem 3. Nature of equilibrium over the axis

a) For all $\sigma = (A, S, C, Q)\in ]0, 1[\times \mathbb{R}_{+}^{3}$

a1. The singularity $\left(1, 0\right)$ is a saddle point.

a2. The equilibrium $\left(0, 0\right)$ is a repeller point.

b) The equilibrium $\left(0, C\right)$ is

(i) a local attractor point if $CQ-A > 0$,

(ii) a hyperbolic saddle point if $CQ-A < 0$,

(iii) a saddle-node if $CQ-A = 0$.

Proof. a1 and a2 are immediate evaluating the Jacobian matrix in each point.

For the point $\left(0, C\right)$ it has

(ⅰ) If $CQ-A > 0$, the $detDY_{\eta }\left(0, C\right)$ is positive and $trDY_{\eta }\left(0, C\right)$ is negative, so the equilibrium point is a local attractor point.

(ⅱ) If $CQ-A < 0$, the $detDY_{\eta }\left(0, C\right)$ is negative, so the equilibrium point is a hyperbolic saddle point.

(ⅲ) If $CQ-A = 0$, the $detDY_{\eta }\left(0, C\right) = 0$ and the equilibrium points $(0, C)$ and $P_1 = (u_1, u_1+C)$ collapse (see Figure 1).

To prove the stability of the singularity $(0, C)$, the center manifold theorem ^[35] will be used.

Setting $(u, v)\rightarrow(X, Y+C)$, system $(2.3)$ is translated to the origin of coordinate; it is given by

The diagonal form of a two–dimensional system can be written by

In addition, the flow on the center manifold is defined by the system of differential equation

Thus, system (3.1) can be written by

So, we have that

Considering the function $h_5(X)$ as the local center manifold defined by

and

In addition, the function $h_5(X)$ satisfies

Thus, replacing (3.3) and (3.4) into equation (3.5), we have

With

Therefore, setting the coefficients $a$, $b$, and $c$ solving the equation (3.6), we have that

Then,

Thus, replacing $h_5(X)$ in equation (3.3); we have that the flow on the center manifold is

With

The series expansion of the function $h_5(X)$, it also is approximate the shape of the local center manifold. Then, the non-hyperbolic equilibrium point $(0, C)$ is an attractor saddle-node.

Remark 4. Assuming $CQ-A < 0$ in the system $\left(2.3\right)$, let us $W^{u}\left(1, 0\right)$, the unstable manifold of the hyperbolic saddle point $\left(1, 0\right)$, and $W^s_\downarrow \left(0, C\right) = \bar{\Sigma}$, the stable manifold of the hyperbolic saddle point $\left(0, C\right)$. Then,

i) the relative position of both manifolds determines a heteroclinic curve, when $W^{u}\left(1, 0\right) \cap \bar{\Sigma}\neq \phi$ (see Figure 2).

ii) when the heteroclinic curve is broken a non-infinitesimal limit cycle is generated, having a heteroclinic bifurcation. This limit cycle could coincide with an infinitesimal limit cycle generated by Hopf bifurcation. The nature (stability) of this non-infinitesimal limit cycle, can be determined evaluating the absolute value of the ratio between the products of the negative and positive eigenvalues of the Jacobian matrix evaluated in both saddle points ^[20]

Theorem 5. The system $\left(2.3\right)$ undergoes a saddle-node bifurcation at $(0, C)$, when $1-A-Q > 0$ and $CQ-A = 0$.

Proof. The Jacobian matrix of the system $(2.3)$ evaluated at the equilibrium point $(0, C)$ is

It is clear to see that $DetDY_{\sigma}(0, C)) = 0$. Moreover, let's represent the dynamical system $(2.3)$ which its vectorial form is given by

.

Let $V = (1, 1)^T$ the eigenvector corresponding to the eigenvalue $\lambda = 0$ of the matrix $DY_{\sigma}(0, C)$. In addition, let $U = (1, 0)^T$ the eigenvector corresponding to the eigenvalue $\lambda = 0$ of the matrix $(DY_{\sigma}(0, C))^T$.

Differentiating the vector function $(2.3)$ with respect to the bifurcation parameter $Q$ we obtain

Therefore,

We will analyze the expression $U[D^2f(u, v, Q)(V, V)]$ where $V = (v_1, v_2)$. The expression inside the parentheses is

Evaluated in the point $(0, C)$ it obtains, $D^2f(0, C, Q)(V, V) = \begin{pmatrix} -2\\ 0 \end{pmatrix}$.

Thus, $U[D^2f(u, v, Q)] = -2\neq0$.

Therefore, by Sotomayor's theorem ^[35], the system $(2.3)$ has a saddle-node bifurcation at $(0, C)$.

Theorem 6. Existence of a homoclinic curve and homoclinic bifurcation

Assuming $CQ-A = 0$ the saddle point $P_{1} = \left(u_{1}, u_{1}+C\right)$ ^[34] coincides with the node attractor $\left(0, C\right)$, obtaining a saddle-node attractor (a non-hyperbolic equilibrium).

a) In the parameter space, there are conditions for which exist a homoclinic curve determined by the central and the unstable manifolds of the attractor saddle-node $\left(0, C\right)$.

b) It exists a non-infinitesimal limit cycle bifurcating from the homoclinic ^[39] surrounding the point $P_{2} = \left(u_{2}, u_{2}+C\right)$.

Proof. Let $W^s_c \left(0, C\right)$, the central stable manifold and $W^u_\nearrow \left(0, C\right)$ the unstable manifold of the non-hyperbolic saddle point $\left(0, C\right)$ (see Figure 3).

a) The trajectory determined by the unstable manifold $W^u_\nearrow \left(0, C\right)$ cannot cross the straight line $u = 1$ towards the right since $\bar{\Gamma}$ is an invariant region (see Lemma 1). This orbit cannot cut or cross the trajectory determined by the central manifold $W^s_c \left(0, C\right)$, by the Existence and Uniqueness Theorem. Moreover, the $\omega -limit$ of the unstable manifold $W^u_\nearrow \left(0, C\right)$ must be:

ⅰ) the point $P_{2}$, when this is an attractor;

ⅱ) a stable limit cycle, if $P_{2}$ is a repeller or if $P_{2}$ is a local attractor, surrounded by two limit cycles (see Figure 6).

On the other hand, the $\alpha -limit$ of the $W^s_c$ $\left(0, C\right)$ can be the sadde point $\left(1, 0\right)$, or lie at infinity in the direction of $u-axis$, or in the point $P_{2}$ if repeller, or in an ustable limit cycle.

Then, there is a subset on the parameter space for which $W_{c}^{s}\left(0, C\right)$ intersects with $W_{+}^{u}\left(0, C\right)$, i.e., $W^s_c\left(0, C\right) \cap W^u_\nearrow \left(0, C\right) \neq \phi$, and a homoclinic curve is obtained. In this case, the same point $\left(0, C\right)$ is the $\alpha$ and the $\omega -limit$ of the right unstable manifold $W^u_\nearrow \left(0, C\right)$.

b) The breaking of the homoclinic curve determined by the intersection of $W^s_c \left(0, C\right)$ and $W^u_\nearrow \left(0, C\right)$ generates a non-infinitesimal limit cycle (existing a homoclinic bifurcation) ^[39], surrounding the point $P_{2} = \left(u_{2}, u_{2}+C\right)$.

Remark 7. 1. To determine the stability of the non-infinitesimal limit cycle, generated by the breaking of the homoclinic curve, it calculates $R$, the absolute value of the ratio between the negative and positive eigenvalues of the Jacobian matrix evaluated in a saddle point, denoted by $\lambda ^{-}$ and $\lambda ^{+}$, respectively, i.e., $R = \left\vert \frac{\lambda ^{-}}{\lambda ^{+}} \right\vert$ ^[20]. If $R = 1$ the limit cycle is neutrally stable ^[20]. Thus, in this case, $R$ depends on the sign of the difference $R_{d} = \left\vert \lambda ^{-}\right\vert -\lambda ^{+}$.

2. When $CQ-A > 0$ the saddle point $P_{1} = \left(u_{1}, u_{1}+C\right)$ exists ^[34]; so, as the changes on the dynamics of the system are smooth, the homoclinic curve determined by this point remains if $(CQ-A) \rightarrow$ 0. Besides, the non-infinitesimal limit cycle generated when $W^s_c \left(0, C\right) \cap W^u_\nearrow \left(0, C\right)$ $= \phi$ has the same stability than the non-infinitesimal limit cycle generated by the breaking of the homoclinic curve determined by $P_{1} = \left(u_{1}, u_{1}+C\right)$. Due to the difficults algebraic resultants when $CQ-A = 0$, we will show the stability of this last periodic orbit using the result obtained in ^[34].

It is easy to see that for $P_{1} = \left(u_{1}, u_{1}+C\right)$ it has

   det$DY_{\lambda }\left(u_{1}, u_{1}+C\right) = Su_{1}\left(C+u_{1}\right) ^{2}\left(A+u\right) \left(2u_{1}-\left(1-A-Q\right) \right)$ and

   tr$DY_{\lambda }\left(u_{1}, u_{1}+C\right) = \left(C+u_{1}\right) \left(u_{1}\left(1-2u_{1}-A\right) -S\left(A+u_{1}\right) \right)$

Then, the eigenvalues of the Jacobian matrix evaluated on $P_{1}$ are:

$\qquad \lambda ^{-} = \frac{1}{2}\left(C+u_{1}\right) \left(\text{tr} DY_{\eta }\left(u_{1}, u_{1}+C\right) -\sqrt{\Delta_{1} } \right)$

$\qquad \lambda ^{+} = \frac{1}{2}\left(C+u_{1}\right) \left(\text{tr} DY_{\eta }\left(u_{1}, u_{1}+C\right) +\sqrt{\Delta_{1}) } \right)$.

with $\lambda ^{-} < 0 < \lambda ^{+}$ and

$\qquad \Delta_{1} = \left(\text{tr}DY_{\eta }\left(u_{1}, u_{1}+C\right) \right) ^{2}-4$det$DY_{\eta }\left(u_{1}, u_{1}+C\right)$.

The following result is obtained

Theorem 8. Stability of the non-infinitesimal limit cycle

The non-infinitesimal limit cycle is:

a) stable, if and only if, $T\left(u_{1}, A.S\right) > 0$, i.e., $S < \frac{ u_{1}\left(1-2u_{1}-A\right) }{A+u_{1}}$.

b) unstable, if and only if, $T\left(u_{1}, A.S\right) < 0$, i.e., $S > \frac{ u_{1}\left(1-2u_{1}-A\right) }{A+u_{1}}$.

c) neutrally stable, if and only if, $T\left(u_{1}, A.S\right) = 0$, i.e., $S = \frac{u_{1}\left(1-2u_{1}-A\right) }{A+u_{1}}$.

Proof. It has that $R_{d} =$ $\left\vert \lambda ^{-}\right\vert -\lambda ^{+} =$ $\left(C+u_{1}\right)$tr$DY_{\eta }\left(u_{1}, u_{1}+C\right)$.

The sign of tr$DY_{\eta }\left(u_{1}, u_{1}+C\right)$ depends on the sign of the factor

$\qquad T\left(u_{1}, A.S, \right) = u_{1}\left(1-2u_{1}-A\right) -S\left(A+u_{1}\right)$.

Considering $R = 1$ or $R_{d} = 0$ we have $S = \frac{u_{1}\left(1-2u_{1}-A\right) }{A+u_{1}}$ and the non-infinitesimal limit cycle is neutrally stable.

The other cases a) and b) are obtained with $R > 1$ and $R < 1$, respectively.

3.1.   Nature of a unique positive equilibrium point

In the following analysis, we consider only the cases 2, 3, 4, 5, recalling being in these cases $CQ-A\leq 0$; the equilibrium $\left(0, C\right)$ is a hyperbolic or a non-hyperbolic saddle point and $P_{1} = \left(u_{1}, u_{1}+C\right)$ lies in the second quadrant or coincide with $\left(0, C\right)$, when $u_{1} = 0$.

We note that the unique positive equilibrium lies in the region

and its nature depends on the relation among $v_{u}$ and $v_{s}$, the ordinate of the points $\left(u^{\ast }, v_{u}\right) \in W^{u}\left(1, 0\right)$ and $\left(u^{\ast }, v_{s}\right) \in \bar{\Sigma}$, respectively.

The Jacobian matrix in the unique positive equilibrium point $\left(u, u+C\right)$ is:

   $DY_{\sigma }\left(u, u+C\right) = \left(C+u\right) \left(\begin{array}{cc} -u\left(A+2u-1\right) & -Qu \\ S\left(A+u\right) & -S\left(A+u\right) \end{array} \right)$.

Thus,

   det$DY_{\sigma }\left(u, u+C\right) = S\left(C+u\right) ^{2}\left(A+u\right) u\left(2u-\left(1-A-Q\right) \right)$, and

   tr$DY_{\sigma }\left(u, u+C\right) = \left(C+u\right) \left(u\left(1-2u-A\right) -S\left(A+u\right) \right)$,

where

$\qquad \qquad u = L = \frac{1}{2}\left(1-A-Q+\sqrt{\Delta }\right)$, with

$\qquad \Delta =$ $\left(1-A-Q\right) ^{2}-4\left(CQ-A\right)$.

Case 2: First, we analize the cases

2ⅰ)   $1-A-Q > 0$ and $CQ-A < 0$, or

2ⅱ)   $1-A-Q < 0$ and $CQ-A < 0$.

In both cases, $\Delta = \left(1-A-Q\right) ^{2}-4\left(CQ-A\right) > 0$.

Then, there exists the unique positive $\left(L, L+C\right)$ in the first quadrant, since $u_{1} < 0$. The equilibrium $\left(0, C\right)$ is a hyperbolic saddle point.

Theorem 9. Nature of the positive equilibrium point Case 2

The equilibrium point $\left(L, L+C\right)$ is

i) an attractor point, if and only if, $S > \frac{L\left(1-2L-A\right) }{L+A}$,

ii) a repeller, if and only if, $S < \frac{L\left(1-2L-A\right) }{L+A}$; furthermore, a Hopf bifurcation occurs and it exists, at least an infinitesimal limit cycle ^[32] surrounding this equilibrium point.

iii) a weak focus, if and only if, $\ S = \frac{L\left(1-2L-A\right) }{L+A}$.

Proof. The Jacobian matrix is

   $DY_{\sigma }\left(L, L+C\right) = \left(L+C\right) \left(\begin{array}{cc} \left(1-2L-A\right) L & -QL \\ S\left(L+A\right) & -S\left(L+A\right) \end{array} \right)$,

and

   det$DY_{\sigma }\left(L, L+C\right) = SL\left(L+C\right) ^{2}\left(L+A\right) \left( 2L+A-1+Q\right) > 0$,

As the factor $2L-\left(1-A-Q\right) > 0$, then, det$DY_{\eta }(L, L+C) > 0$.

So, the nature of $\left(L, L+C\right)$ depends on the sign of

$\qquad$tr$DY_{\sigma }(L, L+C) = \left(C+L\right) \left(L\left(1-2L-A\right) -S\left(L+A\right) \right)$

i.e., this sign depends on the factor

   $T\left(L, A.S\right) = L\left(1-2L-A\right) -S\left(L+A\right)$.

Thus, we have,

ⅰ) tr$DY_{\sigma }\left(L, L+C\right) < 0$, if and only if, $S > \frac{L\left(1-2L-A\right) }{L+A}$; therefore, the point $\left(L, L+C\right)$ is an attractor.

ⅱ) tr$DY_{\sigma }\left(L, L+C\right) > 0$, if and only if, $S < \frac{L\left(1-2L-A\right) }{L+A}$; then, the point $\left(L, L+C\right)$ is a repeller.

As the trace changes its sign, a Hopf bifurcation occurs ^[32,33] at the equilibrium point $\left(L, L+C\right)$. Then, this point is surrounded by a stable limit cycle.

Furthermore, the transversality condition ^[32] is verified, since it has

   $\frac{\partial }{\partial S}\left(\text{tr}DY_{\sigma }\left(L, L+C\right) \right) = -\left(L+A\right)$.

ⅲ) tr$DY_{\sigma }\left(L, L+C\right) = 0$, if and only if, $S = \frac{L\left(1-2L-A\right) }{L+A}$;

thus, the point $\left(L, L+C\right)$ is a fine focus, whose weakness must be determined.

The equilibrium $\left(L, L+C\right)$ can be node or focus depending on the quantity

$\qquad H = \left(\text{tr}DY_{\sigma }(L, L+C)\right) ^{2}-4\left(\det DY_{\sigma }(L, L+C)\right)$

$= \left(C+u_{2}\right) ^{2}\left(\begin{array}{c} \left(A+u_{2}\right) ^{2}S^{2}+ 2u_{2}\left(A+u_{2}\right) \left(A+2\left(1-A-Q\right) -2u_{2}-1\right) S \\ + u_{2}^{2}\left(A+2u_{2}-1\right) ^{2} \end{array} \right)$

The sign of $H$ depends on the factor

$H_{1} = \left(A+u_{2}\right) ^{2}S^{2}+ 2u_{2}\left(A+u_{2}\right) \left(A+2\left(1-A-Q\right) -2u_{2}-1\right) S+ u_{2}^{2}\left(A+2u_{2}-1\right) ^{2}$

As it is known, the point $\left(L, L+C\right)$ is

   ⅰ) a focus, if and only if, $H_{1} < 0$, and

   ⅱ) a node, if and only if, $H_{1} > 0$.

Remark 10. 1. As the trajectories of system $\left(2. 3\right)$ are bounded, and as a consequence of Poincaré-Bendixson Theorem, their $\omega -limit$ can be the unique positive equilibrium point or a stable limit cycle around this point, or else, both a positive equilibrium and stable limit cycle.

The weakness of the fine focus $\left(L, L+C\right)$ must be determined, but we will not make the calculations of the Lyapunov numbers for this case ^[32,36]. Nonetheless, the simulations show the existence of two limit cycles encircling the unique positive equilibrium point, the innermost unstable and the outermost stable.

2. When $S > \frac{\left(1-2L-A\right) L}{A+L}$, in the article by Aziz et al ^[8], a Lyapunov function is constructed proving the point $\left(L, L+C\right)$ is globally asymptotically stable, but this is not always true.

Case 3: Assuming $1-A-Q = 0$ and $CQ-A < 0$.

The equilibrium $\left(0, C\right)$ is saddle point; in the following theorem we establish the nature of the $\left(u_{2}, u_{2}+C\right) = \left(F, F+C\right)$ where $F = \sqrt{ \Delta }$ with $\Delta = -4\left(CQ-A\right) > 0$ and $1-A = Q$.

Hence, it also must fulfill that $A-C\left(1-A\right) = \left(C+1\right) A-C > 0$.

Theorem 11. Nature of the positive equilibrium Case 3

The equilibrium point $\left(F, F+C\right)$ is:

i) a local attractor point, if and only if, $S > \frac{\left(1-2F-A\right) F}{ A+F}$,

ii) a repeller point, if and only if, $S < \frac{\left(1-2F-A\right) F}{A+F}$; furthermore, there exists at least a limit cycle, surrounding $\left(F, F+C\right)$,

iii) a weak focus, if and only if, $S = \frac{\left(1-2F-A\right) F}{A+F}$.

Proof. It is immediate; evaluating the Jacobian matrix it shows that the properties are similar to the above case 2. As $CQ-A < 0$, the equilibrium $\left(0, C\right)$ is a saddle point; then, the $\omega -limit$ of trajectories is the unique positive equilibrium point $\left(F, F+C\right)$ or a stable limit cycle surrounds this point.

Using the Lyapunov function proposed in ^[8], it is possible to prove that if $S > \frac{\left(1-2F-A\right) F }{A+F}$, the point $\left(F, F+C\right)$ is globally asymptotically stable,

Case 4: Assuming $1-A-Q > 0$ and $CQ-A = 0$

In this case, the equation $\left(2.4\right)$ has two solutions

$\qquad u_{1} = 0$ and $u_{2} = \ G = 1-A-Q = 1-A-\frac{A}{C} = \frac{C-A-AC}{C}$.

Then, $\left(G, G+C\right) = \left(\frac{C-A-AC}{C}, \frac{\left(C-A\right) \left(C+1\right) }{C}\right)$ is the unique equilibrium point at interior of the first quadrant, if and only if,

   $C-A-AC = C-\left(C+1\right)A > 0$ and $C-A > 0$.

The point $P_{1} = \left(u_{1}, u_{1}+C\right)$ coincides with $\left(0, C\right)$, being a saddle-node equilibrium (a non-hyperbolic stable equilibrium point). The system $\left(2.3\right)$ can be rewrite as:

where $\rho = \left(A, C, S\right) \in \left] 0, 1\right[ \times \mathbb{R} ^{2}_{+}$.

Theorem 12. Nature of the positive equilibrium Case 4

Let us $\left(u^{\ast}, v_{s}\right)$ $\in W^s_\downarrow \left(0, C\right)$ = $\bar{\Sigma}$, the stable manifold of the saddle point $\left(0, C\right)$ and $\left(u^{\ast}, v_{u}\right)$ $\in W^{u}\left(1, 0\right)$, the unstable manifold of the hyperbolic saddle point $\left(1, 0\right)$ (see Figure 3).

a) Assuming $v_{u} < v_{s}$, the positive equilibrium $\left(G, G+C\right) = \left(\frac{C-A-AC}{C}, \frac{\left(C-A\right) \left(C+1\right) }{C}\right)$ is:

i) a local attractor point, if and only if, $S > \frac{\left(\left(C+2\right) A-C\right) \left(C-A-AC\right) }{C\left(C-A\right) }$.

ii) a repeller point, if and only if, $S < \frac{\left(\left(C+2\right) A-C\right) \left(C-A-AC\right) }{C\left(C-A\right) }$, with $\left(C+2\right) A-C > 0$; moreover, there exists at least a limit cycle surrounding $\left(G, G+C\right)$,

iii) a weak or fine focus, if and only if, $S = \frac{\left(\left(C+2\right) A-C\right) \left(C-A-AC\right) }{C\left(C-A\right) }$, with $\left(C+2\right) A-C > 0$; the weakness must be determined.

b) Assuming $v_{u} > v_{s}$, the positive equilibrium $\left(G, G+C\right)$ is a repeller (focus or node), and $\left(0, 0\right)$ is a almost global stable point ^[37,38].

Proof. a) It is also immediate; evaluating the Jacobian matrix it shows that the properties are similar to the above case.

As in the cases 2 and 3, if $S > \frac{\left(\left(C+2\right) A-C\right) \left(C-A-AC\right) }{C\left(C-A\right) }$, the point $\left(G, G+C\right)$ is a local attractor.

b) When $S < \frac{\left(\left(C+2\right) A-C\right) \left(C-A-AC\right) }{C\left(C-A\right) }$, there exists a limit cycle, whose diameter increases until to coincide with the heteroclinic curve joining the equilibrium points $\left(0, C\right)$ and $\left(1, 0\right)$.

Laterly, this heteroclinic breaks up and the singularity $\left(0, C\right)$ is an attractor for almost all trajectories ^[37,38].

Remark 13. When $v_{u} > v_{s}$, the equilibrium point $\left(0, C\right)$ is an almost global stable point ^[37,38], since there exists the unique singularity $\left(G, G+C\right)$ at the interior of the first quadrant. Thus, no necessary is globally asymptotically stable

Theorem 14. Weakness of the positive equilibrium

The singularity $\left(G, G+C\right) = \left(\frac{C-A-AC}{C}, \frac{\left(C-A\right) \left(C+1\right) }{C}\right)$ of the vector field $Z_{\rho }$ is at least a two order weak focus, if and only if, $S = \frac{\left(2A-C+AC\right) \left(C-A-AC\right) }{C\left(C-A\right) }$.

Proof. As $S = \frac{\left(2A-C+AC\right) \left(C-A-AC\right) }{C\left(C-A\right) }$ and $Q = \frac{A}{C}$, system $\left(2.3\right)$ can be expressed by

where $\nu = \left(A, C\right) \in \left] 0, 1\right[ \times \mathbb{R} ^{+}$. The Jacobian matrix is

   $DZ_{\nu }(u, u+C) = \left(\begin{array}{cc} u\left(u+C\right) \left(1-2u-A\right) & -\frac{A}{C}u\left(u+C\right) \\ \frac{\left(2A-C+AC\right) \left(C-A-AC\right) \left(u+C\right) \left(A+u\right) }{C\left(C-A\right) } & -\frac{\left(2A-C+AC\right) \left(C-A-AC\right) \left(A+u\right) \left(C+u\right) }{C\left(C-A\right) } \end{array} \right)$

and evaluated in the equilibrium $\left(G, G+C\right) = \left(\frac{C-A-AC}{C }, \frac{\left(C-A\right) \left(C+1\right) }{C}\right)$, it has

   $DZ_{\nu }(G, G+C) = \frac{\left(C-A\right) \left(C+1\right) \left(C-A-AC\right) }{C^{3}}\left(\begin{array}{cc} 2A-C+AC & -A \\ 2A-C+AC & -\left(2A-C+AC\right) \end{array} \right)$.

Setting $u = U+G$ and $v = V+G+C$; then, the new system translated to origin of coordinates system $\left(U, V\right)$ after algebraic manipulations is

The Jacobian matrix of the vector field $\bar{Z}_{\nu }\left(U, V\right)$ at the point $\left(0, 0\right)$ is the same than the Jacobian matrix of the vector field $Z_{\nu }\left(u, v\right)$ at the point $(G, G+C)$.

We denote

$\qquad W^{2} =$ det$D\bar{Z}_{\nu }(0, 0) =$ $R^{2}\left(C-A-AC\right) \left(2A-C+AC\right)$,

where $R = \frac{\left(C-A\right) \left(C+1\right) \left(C-A-AC\right) }{C^{3}}$, with $C-A-AC > 0$.

The first Lyapunov quantity ^[32] is $\eta _{1} =$ tr$DZ_{\nu }\left(G, G+C\right) =$ tr$D\bar{Z}_{\nu }(0, 0) =$ $0$. To determine the normal form of vector field $\bar{Z}_{\nu }$, we use the procedure described in ^[23].

The Jordan matrix associated to vector field $\bar{Z}_{\nu }$ is $J = \left(\begin{array}{rc} 0 & - W \\ W & 0 \end{array} \right)$.

The matrix change of basis ^[23] is

$\qquad M = \left(\begin{array}{rc} \bar{Z}_{11}-\text{tr}D\bar{Z}_{\nu } & -\text{det}D\bar{Z}_{\nu } \\ \bar{Z}_{21} & 0 \end{array} \right) = \left(\begin{array}{rc} R^{2}\left(2A-C+AC\right) & -W \\ R^{2}\left(2A-C+AC\right) & 0 \end{array} \right)$.

Now consider the change of variables given by

   $\left(\begin{array}{c} U \\ V \end{array} \right) = M\left(\begin{array}{c} x \\ y \end{array} \right)$

that is,

$\qquad U = R^{2}\left(2A-C+AC\right) x-Wy$

$\qquad V = R^{2}\left(2A-C+AC\right) x$

or

$\qquad x =$ $\frac{1}{R^{2}\left(2A-C+AC\right) }V\qquad \qquad$

$\qquad y = \frac{-U+V}{W}$

Then the new system is

$\qquad \tilde{Z}_{\nu }:\left\{ \begin{array}{lcl} \frac{dx}{d\tau } & = & \frac{1}{S\left(A+G\right) }\frac{dV}{d\tau } \\ \frac{dy}{d\tau } & = & \frac{1}{W}\left(-\frac{dU}{d\tau }+\frac{dV}{d\tau }\right) \end{array} \right.$

After a large algebraic calculations we obtain

   $\breve{Z}_{\nu }:\left\{ \begin{array}{lcl} \frac{dx}{d\tau } & = & \begin{array}{l} -W\left(C+G\right) y-WS\left(A+C+2G\right) xy+ \frac{ W^{2}\left(C+G\right) }{A+G}y^{2} \\ -S^{2}W\left(A+G\right) x^{2}y+S^{2}W^{2}xy^{2} \end{array} \\ \frac{dy}{d\tau } & = & \begin{array}{l} \frac{W}{C+G}x+ \frac{W^{3}}{G^{3}}\frac{ 2C+3G}{\left(C+G\right) ^{4}}x^{2}-S\left(A+G\right) \left(AS+CS+2GS+2\right) xy \\ +W\left(CS+GS+1\right) y^{2}+ S^{3}\left(A+G\right) ^{3}\frac{A-C+AC+4CG+C^{2}}{C}x^{3} \\ -S^{3}\left(A+G\right) ^{2}x^{2}y +S^{2}W\left(A+G\right) xy^{2}-W^{2}\left(A+C+4G-1\right) y^{3} \\ +\frac{S^{4}}{W}\left(A+G\right) ^{4}x^{4}-4S^{3}\left(A+G\right) ^{3}x^{3}y+ GS^{2}\left(2C+3G\right) \frac{\left(A+G\right) ^{2} }{W}x^{2}y^{2} \\ -4SW^{2}\left(A+G\right) xy^{3}+W^{3}y^{4} \end{array} \end{array} \right.$

To obtain the normal form we make the change of variables given by

$\begin{array}{lcl} \qquad w & = & Hx \end{array}$ and $\begin{array}{lcl} z & = & Jy \end{array}$

with H and J parameter to detetermine.

Thus,

$\begin{array}{lcl} \qquad \frac{dw}{d\tau } & = & H\frac{dx}{d\tau } \end{array}$ and $\begin{array}{lcl} \frac{dz}{d\tau } & = & J\frac{dy}{d\tau } \end{array}$.

So,

$\begin{array}{lcl} \qquad \frac{dw}{d\tau } & = & -HW\left(C+G\right) \frac{1}{J}z+HOT \\ \qquad \frac{dz}{dT} & = & J\frac{W}{\left(C+G\right) }\frac{1}{H}w+HOT \end{array}$

As $HW\left(C+G\right) \frac{1}{J} = J\frac{W}{\left(C+G\right) }\frac{1}{H} = 1$, then $J = H\left(C+G\right)$, obtaining

$\begin{array}{lcl} \qquad \frac{dw}{d\tau } & = & -Wz+HOT \\ \qquad \frac{dz}{d\tau } & = & Ww+HOT. \end{array}$

Considering the time rescaling given by $T = W\tau$

   $\frac{dw}{d\tau } = \frac{dw}{dT}\frac{dT}{d\tau } = W\frac{dw}{dT}$ and $\frac{ dz}{d\tau } = \frac{dz}{dT}\frac{dT}{d\tau }$

and finally the normal form ^[32] of vector field $Z_{\nu }$

   $\tilde{Z}_{\nu }:\left\{ \begin{array}{lcl} \frac{dw}{dT} & = & -z+HOT \\ \frac{dz}{dT} & = & w+HOT \end{array} \right.$

Using the Mathematica package ^[41] for the simbolic calculus, we obtain that the second Lyapunov quantity ^[32] given by $\eta _{2} = \frac{G^{4}\left(-1+A+2G\right) ^{2}\left(C+G\right) ^{2}}{ 8\left(A+G\right) ^{2}W^{3}}f_{1}\left(A, C, G\right)$

with $f_{1}\left(A, C, G\right) = f_{10}\left(A, G\right) +Cf_{11}\left(A, G\right) +C^{2}f_{12}\left(A, G\right) +C^{3}f_{13}\left(A, G\right)$

where

$f_{10}\left(A, G\right) = -8G^{7}+ 4\left(1-2A\right) G^{6}-A\left(A+11\right) G^{5}+ A\left(A^{2}-8A+5\right) G^{4}+ A^{2}\left(7A-34\right) G^{3}+A^{2}\left(7A^{2}+4\right) G^{2}+A^{3}\left(A^{2}-10A+3\right) G+A^{4}\left(1-A\right)$.

$f_{11}\left(A, G\right) = -36G^{6}+16\left(2-3A\right) G^{5}-\left(13A^{2}+15A+6\right) G^{4}+A\left(11A^{2}-64A+37\right) G^{3}+A\left(7A^{3}-46A^{2}+39A-8\right) G^{2}+A^{2}\left(1-A\right) \left(11A-A^{2}-4\right) G+\left(A-A^{2}-2\right) A^{3}$.

$f_{12}\left(A, G\right) = -50G^{5}+22\left(3-4A\right) G^{4}+\left(75A-61A^{2}-28\right) G^{3}+ \left(33A^{2}-19A^{3}-12A+4\right) G^{2}+ A\left(1-A\right) \left(-4A+2A^{2}-1\right) G+A^{2}\left(1-A\right)$.

$f_{13}\left(A, G\right) = \left(-16G^{3}+24\left(1-A\right) G^{2}-12\left(1-A\right) ^{2}G+ 2\left(1-A\right) ^{3}\right) G$.

Evaluating the factor $f_{1}\left(A, C, G\right)$ for $G = 0, 6$, $C = 0.1$ and after for $G = 0.5$, $C = 0.25$ considering differents values for $A$, it can see this factor can change its sign; then, $\eta _{2}$ changes its sign for some values of the parameters $A$, $C$ and $G$.

Therefore, the weakness of $\left(G, G+C\right)$ is two; thus, it can exist up two limit cycles around of the unique positive equilibrium point, the innermost unstable and the outermost stable.

Remark 15. The above calculations of the weakness of the focus $\left(G, G+C\right)$ requires two strict relationships in the parameter space; the variation of one of them originates the existence of one or two limit cycles. So, these two relations determine a subset of measure non-zero in the parameter space where the existence of two limit cycles can be assured.

Applying new bifurcation methods and geometric approaches developed in ^[40], it can be completed the qualitative analysis effected in the current work.

Nevertheless, any tiny change in some of the involved parameters will imply inequality rather than equality; hence, we have the structural instability of the system.

This is ecologically important, inasmuch as in reality, none of the equalities given in the above Theorem will be possible to maintain by a long time.

Case 5: Non-existence of positive equilibrium points

The condition on the parameter values are:

   5.1 $1-A-Q = 0$ and $CQ-A > 0$, or

   5.2 $1-A-Q < 0$ and $CQ-A\geq 0$.

   5.3 $1-A-Q > 0$ and $CQ-A > 0$ and $\Delta < 0$.

Theorem 16. No positive equilibrium

When there exist no positive equilibria, the point $\left(0, C\right)$ is globally asymptotically stable.

Proof. By lemma $5,$ we know the solutions are bounded; by lemma 3, $\bar{\Gamma}$ is invariant region. Reminding besides the equilibrium $\left(1, 0\right)$ is a saddle point, then the Poincaré-Bendixson Theorem applies and the unique $\omega -limit$ of the trajectories is the point $\left(0, C\right)$.

4.   Numerical simulations and a bifurcation diagram

In order to reinforce the obtained results we show some numerical simulations. The cases 2, 3, 4 and 5 have been considered, recalling in those cases the point $\left(0, C\right)$ is a hyperbolic or non-hyperbolic saddle point.

Case 2. Assuming $1-A-Q > 0$ and $CQ-A < 0$.

This relations are fulfill in the Figures 4 to 9, for $A = 0.2000005$, $C = 0.39999925$, $Q = 0.5$. By Lemma 3.3 the singularity $\left(0, C\right)$ is a hyperbolic saddle point. However, in some simulations the equilibrium $\left(0, C\right)$ behaves like a local attractor for some trajectories, without being.

The Hopf bifurcation is originated when $S = \frac{\left(\left(C+2\right) A-C\right) \left(C-A-AC\right) }{ C\left(C-A\right) }$

Case 3. Assuming $1-A-Q = 0$ and $QC-A < 0$.

Remerbering that in this case the equilibrium $\left(0, C\right)$ is also a saddle point.

This is the case in which the dynamic of system is most simple (see Figures 10 to 12).

This relations on the parameters are fulfill choosing the following values $C = 0.2$; $Q = 0.825$; $A = 0.175$, having that $CQ-A = \left(0.825\right) \left(0.2\right) -0.175 = -0.01 < 0$ and $1-0.175-0.825 = 0$.

Case 4. Assuming $1-A-Q > 0$ and $QC-A = 0$.

The point $P_{1}$ coincides with $\left(0, C\right)$, being a non-hyperbolic saddle point and it has a saddle-node bifurcation.

We consider the values $A = 0.2$; $C = 0.4$; $Q = 0.5$, for the following simulations (Figures 13 to 18); thus, $1-A-Q = 1-0.2-$ $0.5 = 0.3 > 0$ and $CQ-A = \left(0.4\right) \left(0.5\right) -0.2 = 0$.

The Hopf bifurcation is originated when

   $S = \frac{\left(\left(\left(0.4\right) +2\right) \left(0.2\right) -\left(0.4\right) \right) \left(\left(0.4\right) -\left(0.2\right) -\left(0.2\right) \left(0.4\right) \right) }{\left(0.4\right) \left(\left(0.4\right) -\left(0.2\right) \right) } = 0.12$,

then, the positive equilibrium is a weak focus.

If $S > 0.12$, the positive equilibrium is an attractor.

If $S < 0.12$, the positive equilibrium is a repeller.

Case 5. The point $\left(0, C\right)$ is a global attractor, and there no exists positive equilibrium points.

4.1.   Bifurcation diagram

To get the bifurcation diagram of the system $\left(2.3\right)$ for $A$ and $Q$ fixed we use the numerical bifurcation package ^[42], following to ^[43].

Furthermore, the bifurcation curves obtained from Lemma $3$ and Theorems $5-13$ divide the $\left(C, S\right)$ parameter space into seven parts. Modifying the parameter $C$ impacts the number of positive equilibrium points and the stability of the equilibrium point $\left(0, C\right)$ of system $\left(2.3\right)$.

The modification of the parameter $S$ changes the stability of the positive equilibrium point $P_{2}$ of system $\left(2.3\right)$, while the equilibrium points $\left(0, 0\right)$ and $\left(1, 0\right)$ do not change their behavior.

When the parameters $C, S$ are located in Region $I$ (red area), there are no positive equilibrium points in system $\left(2.3\right)$. If the parameters $C$ and $S$ are located in this area, where $\Delta < 0$, the point $\left(0, C\right)$ is globally asymptotically stable.

When $C = C^{\ast \ast }$ the equilibrium point $P_{1}$ and $P_{2}$ collapse ^[34], see $SN_{1}$ in Figure 21. So that, system $\left(2.3\right)$ experiences a Bogdanov-Takens bifurcation ^[34].

If the parameter $C^{\ast } < C < C^{\ast \ast }$ and $S$ are located in the Region $II$ (dark green area), Region $IV$ (dark grey area) or Region $VI$ (dark blue area), system $\left(2.3\right)$ has two equilibrium points $P_{1}$ and $P_{2}$. The equilibrium point $P_{1}$ is always a saddle point ^[34], while $P_{2}$ can be unstable (Region $VI$) or stable (Region $II$).

For $C$ and $S$ in Region $IV_{1}$ the equilibrium point $P_{2}$ is stable surrounded by two limit cycles, while when $C$ and $S$ are located in Region $IV_{2}$ the equilibrium point $P_{2}$ is unstable surrounded by a stable limit cycle.

Similarly, when $C$ and $S$ are in Region $V_{1}$ the equilibrium point $P_{1}$ is in the second or third quadrant and the equilibrium point $P_{2}$ is stable surrounded by two limit cycle, while when $C$ and $S$ are located in Region $V_{2}$ the equilibrium point $P_{1}$ is also in the second or third quadrant and the equilibrium point $P_{2}$ is unstable surrounded by a stable limit cycle.

Additionally, when $C = C^{\ast }$ the equilibrium point $P_{1}$ and $\left(0, C\right)$ collapse, then system $\left(2.3\right)$ experience a saddle-node bifurcation, see $SN_{2}$ in Figure 21.

If $C < C^{\ast }$ then the equilibrium point $P_{1}$ crosses to the second or third quadrant, so that the equilibrium point $\left(0, C\right)$ is a saddle point and $P_{2}$ is the only positive equilibrium point.

Furthermore, $P_{2}$ can be unstable (Region $VII$) or stable (Region $III$). For $C$ and $S$ in Region $V$ the unstable equilibrium point $P_{2}$ is surrounded by a stable limit cycle.

5.   Conclusions

In this work, a modified Leslie-Gower model ^[8], particularly a modified May-Holling-Tanner predator-prey model was studied, considering that predators can eat other prey in the case of severe scarcity of its most favorite food. This situation was taken account by adding a positive constant $c$ in the function $K_{y}$ representing the environmental carrying capacity for predators. This implies the existence of a new equilibrium point $\left(0, c\right)$ in the $y-axis$.

By means of a diffeomorphism, we obtained the topologically equivalent system $\left(2.3\right)$, dependent only of four parameters listed as $A$, $C$, $Q$ and $S$. We show that in this model (and in the original one) there are five different dynamical situations according to the relation among the parameters $A$, $C$, and $Q$, specifically on the term $QC-A$.

Although the results obtained in the commented paper ^[8] are correct and valid, they are only for an important and special case. in this work, we analyzed the case named 2, 3, 4 and 5, not considered in the mentioned paper. In all these cases there exist only one positive equilibrium points in the system $\left(2.3\right)$ Furthermore, it exists the point $\left(0, C\right)$ over the $v-axis$ and the equilibrium $\left(0, 0\right)$ and $\left(1, 0\right)$.

We show that the model has a rich dynamic since this model can exhibit various kinds of bifurcations (e.g. saddle-node, Hopf-Andronov, homoclinic, Hopf multiple bifurcations ^[23,32,33]) as likewise infinitesimal and non-infinitesimal limit cycles, generated by Hopf and homoclinic bifurcation, respectively.

Conditions for the existence of equilibrium points and their nature was established. We proved that the equilibrium point $\left(0, 0\right)$ is always a repeller for all parameter values, which means that there is no extinction of both populations simultaneously. Moreover, $\left(1, 0\right)$ is a saddle point, implying that the predator population can go to depletion; meanwhile, the prey attains its maximum population size in the common environmental.

We also prove the existence of a homoclinic curve determined by the stable and unstable manifolds of the positive saddle point $\left(0, C\right)$, encircling the unique positive equilibrium point $P_{2}$; when it breaks up originates a non-infinitesimal limit cycle.

The dynamics of the studied model, in which the predators have as an alternative food to low densities of prey, differs from the May-Holling-Tanner model ^[23,25], since in system $\left(2.1\right)$:

ⅰ) can have one, two or none positive equilibrium points at the interior of the first quadrant with a more varied dynamic; whereas, the May-Holling-Tanner model has a unique positive equilibrium point, which never can be a cusp point (Bogdanov-Takens bifurcation).

ⅱ) There exist parameter constraints for which a homoclinic curve exists, something that not appears in the May-Holling-Tanner model.

ⅲ) It has a separatrix curve dividing the behavior of the trajectories, which are originated by, the non-hyperbolic saddle point $\left(0, 0\right)$, in the May-Holling-Tanner model, and in the modified model is created by the hyperbolic attractor point $\left(0, C\right)$.

ⅳ) There exist the bi-stability phenomenon since two limit cycles can bifurcate of a weak focus, surrounding an attractor positive equilibrium point, being the innermost unstable (frontier of the attraction basin) and the outermost stable.

However, taking account of the complexity of the model analyzed, one may ask if, in nature, there exists a predator-prey interaction with these characteristics. Moreover, these complex dynamics are a prominent issue to be considered by the ecologists and agencies responsible for conservation and management of renewable resources, especially in those populations sensitive to disturbances of the environment

In short, in this paper we extend the dynamical properties of the model proposed in ^[8] and the partial results obtained in previous papers ^{[11,12,14,15,16]}; we also complement the outcomes obtained in ^[25] for the Holling-Tanner model, showing that the modified model has interesting and rich mathematical and ecological behaviors.

Acknowledgments

EGO was partially financed by DIEA-PUCV 124.730/2012.

Conflict of interest

All authors declare no conflicts of interest in this paper.