Global stability in a modified Leslie-Gower type predation model assuming mutual interference among generalist predators

1.   Introduction

The dynamic relationship between predators and their prey has been and will continue to be one of the dominant themes in Theoretical Ecology (Mathematical Ecology) and particularly in Population Dynamics. This is due to its universal existence ^[1] and because the more exhaustive knowledge of this interaction allows a better understanding of the behavior of food chains or trophic webs ^[2,3].

It is well known that the first predator-prey model was proposed by the Italian mathematician Vito Volterra ^[4,5] in a well-known monograph published in 1926 ^[6], being described by an autonomous nonlinear ordinary differential equation system (ODEs). This model coincided with a two-dimensional model for biochemical interactions previously proposed by the American mathematician-physicist Alfred J. Lotka; for this, the ODE system is known by the name of Lotka-Volterra model ^[2,3,5,7].

The main dynamic characteristic of this first predation model is that the only point of positive equilibrium is a center ^[8,9], that is, all the paths are concentric closed orbits around that point ^[2,3,5,7]. This implies that the population sizes of predators and their prey would permanently oscillate around that point for any initial condition ^[3,7]. This behavior of the system solutions was strongly questioned after its formulation because in nature no predator-prey interactions with these characteristics are found ^[2,3].

Based on Volterra's work, different proposals came up to face and resolve the various objections made ^[10]. One of the first proposals to solve some of the objections to the Lotka-Volterra model was raised by the Russian biologist Georgii F. Gause (1910–1986). In 1934 he proposed a model that takes into account the intraspecific competition on the prey population ^[11], replacing the Malthusian prey growth incorporated in the Lotka-Volterra model, by the logistic growth function (also called the Verhulst-Pearl equation ^[4]).

Another important subsequent result is the Theorem proposed by the Russian mathematician Andrei N. Kolmogorov in 1936 ^[2,5]. In this, conditions for a general ODE system describing a predator-prey interaction are established, to ensure the existence of a single stable periodic solution ^[12] (mathematically, an attractor limit cycle ^[8,9]). For the phenomenon of oscillations of population sizes of prey and predators, there is sufficient evidence in nature, which are called ecologically stable cycles ^[12].

A different alternative is presented in the model formulated by the British ecologist Patrick Holt Leslie in 1948 ^[13], which does not fit the Lotka-Volterra model scheme ^[3]. Unlike the Gause-type compartmentalized models ^[3,14], based on a principle of mass or energy transfer ^[1], the Leslie model characterizes because the predator growth equation is of the logistic type, just like the prey growth equation ^[3,4].

Leslie assumed that the conventional environmental carrying capacity of the predators $K_{y}$ is proportional to the abundance of prey $x$ ^[2], i.e., it is assumed that $K_{y} = K(x) = nx$ ^{[2,3,15,16,17]}.

However, when the predator is a generalist ^[3] and there is no favorite prey available, predators can switch to a different food source. In this case, the predator environmental carrying capacity can be expressed as $K(x) = nx+c$, where the parameter $c > 0$ indicates the amount of alternative food available to predators ^[18,19]. Thus, a modified Leslie-Gower or a Leslie-Gower scheme is obtained ^[20,21].

On the other hand, it is known that the action of predators in the interaction is called the functional response of predators or consumption function ^[14]; it refers to the change in prey density attacked per unit time per predator, when prey density changes ^[14]. They are classified into various types, depending on the size of the prey population or the size of both populations, and there are various mathematical ways to describe them.

Based on laboratory experiments, the Canadian ecologist Crawford S. Holling (1930–2019) in 1959 ^[22] described three types of saturated functions, which he considered only dependent on the size of the prey population (prey-dependent functional response).

Later, in 1984, Robert J. Taylor ^[23], proposed the non-monotonic functional response, usually used to describe the ecological phenomenon called group defense formation ^[23] and others, as the agammaegation and inhibition for microorganism interactions. Prey-dependent functional responses are classified as Holling-type, Ⅰ, Ⅱ, Ⅲ, and Ⅳ, and there are various mathematical ways to describe each type ^[2,3]. For example, in the so-called May-Holling-Tanner model ^{[2,3,17,20,24,25]}, which is derived from the Leslie model, the hyperbolic functional response is considered ^[18,25]. Functional responses, depending on both populations have also been proposed, such as the Beddington-DeAngelis ^[26,27], the ratio-dependent ^[3], etc.

In the Leslie model, the functional response of the predator is expressed by the linear function $h(x) = \ qx$, also used in the classic Lotka-Volterra model ^[2,3]. Nevertheless, this functional response does not correspond to any of the types proposed by Holling, since it is not bounded.

In this work, the ecological phenomenon called interference or competition among predators (CAP) is incorporated in the Leslie model, which means that two or more predators are bothered when capturing prey and therefore the prey can avoid predation ^[28].

A similar situation can be assumed in open access bioeconomic models when fishermen go to the same place with their boats to harvesting fish ^[29,30]; the boats can collide with each other and even cause agammaessions between the crews (making shots, breaking fishing implements, etc.), as has happened in various places on the planet ^[29,30].

There are three ways to mathematically represent CAP. In this work we will assume the expression formulated by the Canadian mathematician Herbert I. Freedman in 1979 ^[31], modifying the assumption of the usual models to express the functional response. Freedman formulated a change in the predation function in Gause-type models, proposing the expression:

power function, where $\beta$ is the constant of mutual interference, and $h\left(x\right)$ is the functional response of the predator dependent only on the prey population. He also included a second mathematical form for the CAP in the Gause type models by adding the negative term $m\left(y\right) = -ey^{2}$ in the equation that describes the growth of predators ^[5].

The power function appears also proposed in the bioeconomic literature ^[29] as a particular case of a more general function called Cobb-Douglas type production function ^[30,32,33], described by $G\left(x, y\right) = qx^{\alpha }y^{\beta }$, with $0 < \alpha, \beta < 1$ and $q > 0$.

In Figure 1, we see the behavior of the function $g\left(y\right) = y^{\beta }$, with $0 < \beta < 1$, for different values of $\beta$.

It can observe that if $\beta$ $\rightarrow 0$, the graphic of the function grows more abruptly. This means that the competition among predators is more intense, whereas for large values of $\beta$ ($\beta$ $\rightarrow 1$) the competition decreases. When $\beta = 1$, we have the well-known Leslie-Gower model ^[34].

As we will show, the function $g\left(x, y\right)$ has no strong implication on the dynamics of the system, such as happens when is incorporated in the Volterra model ^[31].

Other alternative mathematical ways to model CAP are: the Beddington-DeAngelis functional response ^[3,26,27], dependent on both population sizes, and the addition of a quadratic term in the growth equation of predators, which is incorparated mainly in Gause-type models ^[5,28,31]. Nonetheless, according to our knowledge, the modified Leslie-Gower type model presented here considering a power function describing CAP has not been previously analyzed.

This work is focused on the analysis of two models derived from the Leslie-Gower model, assuming both have an alternative food, but in one of them, we consider a predator compete with each other; the differences with the models in which other ways to describe CAP are considered will be show.

The rest of the article is organized as follows: In Section 2 we will present the modified Leslie-Gower type model considering alternative food and its basic properties. In Section 3 we will show the main properties of the model assuming CAP and an alternative supply to this species and in Section 4 we will discuss the consequences of CAP on the modified Leslie-Gower type model.

2.   The Leslie model with generalist predators

The Leslie-Gower model considering a generalist or alternating food for predators is described by the following Kolmogorov-type system ^[14,36]:

where $x = x\left(t\right)$ and $y = y\left(t\right)$ indicate the population sizes of prey and predator, respectively, at any time $t\geq 0$, subject to $x\left(0\right) > 0$ and $y\left(0\right) > 0$, measured in number of individuals, biomass or density per unit area or volume, with $\sigma = (r, k, q, s, n, c)\in \mathbb{R}_{+}^{6}$. The parameters have the ecological meanings described in the following Table 1:

The parameter $c > 0$ characterizes to generalist predators ^[20,21], because in the absence of his favorite prey (when $x = 0$), it searches for an alternative food source ^[3,19]. If $c = 0$, the original model formulated by Leslie in 1948 is obtained ^[13,34] and the predator is said to be a specialist ^[3]

According to Turchin ^[3], when the favorite prey population size is low, generalist predators should focus on other prey species. When the favorite prey population size is high, predators will switch to hunting it, because it becomes profitable for them to do so ^[3].

In some articles, the parameter $c > 0$ is assumed to measure the extent to which the environment provides protection to predators, but we prefer the annotated meaning in Table 1.

Another form to describe the behavior of generalist predators is using a Holling Ⅲ functional response ^[5,2,3]. In ^[35] a Gause type predator-prey model was studied, considering a sigmoid functional response and complex dynamics were obtained as bi-stability, limit-cycles, and bifurcations.

It should be noted that in model $\left(1\right)$,

(ⅰ) In the absence of favorite prey, when $x = 0$, the growth of the predator population is assured.

(ⅱ) If $y >$ $\frac{r}{q}\left(1-\frac{x}{K}\right)$, then $\frac{dx }{dt} < 0$ and the prey population becomes extinct.

(ⅲ) If $y > nx+c$, then $\frac{dy}{dt} < 0$ and the predator population goes to extinction.

The domain of the system $\left(1\right)$ or vector field $G_{\sigma }\left(x, y\right)$ is the set

Unlike the original Leslie-Gower model, system $\left(1\right)$ is defined at $x = 0$ and has four equilibrium points or singularities. They are: $\left(0, 0\right)$, $\left(K, 0\right)$, $\left(0, c\right)$ and $\left(x_{e}, nx_{e}+c\right)$ with $x_{e} = \frac{\left(r-cq\right) K}{r+nqK}$.

Formally for this last equilibrium point, we have the following lemma

Lemma 1. The equilibrium point $\left(x_{e}, nx_{e}+c\right)$

i) belongs to the interior of the first quadrant and is unique, if and only if, $c < \frac{r}{q}$,

ii) coincides with the equilibrium point $\left(0, c\right)$, if and only if, $c = \frac{r}{q}$,

iii) is outside of the first quadrant, if and only if, $c > \frac{r}{q}$.

Proof. The equilibrium point $\left(x_{e}, nx_{e}+c\right)$ is at the intersection of the isoclines and the abscissa $x_{e}$ satisfies the equations

Both are straight lines.

$f_{1}$ is a decreasing function and passes through the points $\left(0, \frac{r}{q}\right)$ and $\left(K, 0\right)$,

$f_{2}$ is an increasing function that passes through the point $\left(0, c\right)$.

Thus, the graphs of $f_{1}\left(x\right)$ and $f_{2}\left(x\right)$ intersect inside the first quadrant $Int\left(\mathbb{R} _{0}^{+}\right) ^{2}$, if and only if, $c < \frac{r}{q}$

Clearly, if $c = \frac{r}{q}$, the equilibrium point $\left(x_{e}, nx_{e}+c\right)$ coincides with the point $\left(0, c\right)$, and is outside the first quadrant, if and only if, $c > \frac{r}{q}$.

To determine the nature of the hyperbolic equilibrium points we require the Jacobian matrix, which is

Thus, it is defined for every point $\left(x, y\right) \in \Omega$.

2.1.   Main properties of the system $\left(1\right)$

The main properties of the modified Leslie-Gower model described by the system $(1)$ or vector field $G_{\sigma }\left(x, y\right)$ are shown below.

Lemma 2. The set

is a positively invariant region.

Proof. Clearly, the coordinate axes are invariant sets

Let's consider $x = K$; then, we get that

Thus, $\frac{dx}{dt} < 0$. Then, the orbits with an initial condition outside the set $\Gamma$ will enter into it; the trajectories with initial conditions within the set $\Gamma$ cannot leave it, whatever the sign of $\frac{dy}{dt}$.

Remark 3. We note that the subregion

is also a positively invariant region for system $\left(1\right)$.

Lemma 4. The solutions are uniformly bounded.

Proof. We will use the comparison theorem for differential inequalities ^[41].

From the first equation we have

We also know that

Let $L = max\left\{ x\left(0\right), K\right\}$; thus $x\left(t\right) \leq L$, for all $t\geq 0$.

Also, $0 < y < nx+c$.

Let's consider the variable $w\left(x, y\right) = \ x+\frac{1}{s}y$. Clearly, $w\left(x, y\right) > 0$, for all $t\geq 0$.

Then,

Now consider the sum of $\frac{dw}{dt}$ and $\sigma w$, with $\sigma > 0$, obtaining

Thus,

Let $M =$ $\frac{rK}{4}+\sigma \ x+\left(\frac{\sigma }{s}+1\right) y$; hence

which is a first order differential inequality.

Applying the comparison theorem for differential inequalities (Page 30 in ^[41]), we get that

For $t = 0$, we obtain that $w\left(0\right) \leq M+H$; that is, $H\geq$ $w\left(0\right) -M$.

Then, there is $m\in \mathbb{N}$ such that $H\leq m$ $\left(w\left(0\right) -M\right)$;

therefore, $w\left(t\right) e^{t}\leq Me^{t}+m\left(w\left(0\right) -M\right)$, i.e., $w\left(t\right) \leq M+m\left(w\left(0\right) -M\right) e^{-t}$

Therefore, when $t\rightarrow \infty$ then $w\left(t\right) \leq M$.

Thus, the solutions are bounded.

Remark 5. The previous result ensures that there is a set

which is the region where all the solutions of the system $\left(1\right)$ with initial conditions in $\Omega$ are confined ^[36].

This property ensures that the model is well-posed ^[1].

Therefore, many predators and a small prey population size may coexist; but this would be a contradictory situation if the predators are specialists.

Lemma 6. The equilibrium point $\left(0, 0\right)$ is a hyperbolic repeller for all parameter values.

Proof. The Jacobian matrix evaluated at point $\left(0, 0\right)$ is

Thus,

Applying the trace and determinant theorem ^[7,37], the thesis is held.

Lemma 7. The singularity $\left(K, 0\right)$ is hyperbolic saddle point for any set of parameters.

Proof. The Jacobian matrix evaluated at point $\left(K, 0\right)$ is

and det$DG_{\sigma }\left(K, 0\right) = -rs\left(c+nK\right) < 0$.

According to the trace and determinant theorem ^[7,37] the thesis is obtained.

Theorem 8. The equilibrium point $\left(0, c\right)$ is

(ⅰ) a hyperbolic saddle point, if and only if, $c < \frac{r}{q}$.

(ⅱ) a non-hyperbolic saddle point, if and only if, $c = \frac{r}{ q}$.

(ⅲ) a hyperbolic attractor, if and only if, $c > \frac{r}{q}$.

Proof. The Jacobian matrix evaluated at point $\left(0, c\right)$ is

Thus, det$DG_{\sigma }\left(0, 0\right) = -s\left( r-cq\right)$, whose sign depends on the factor $b = r-cq$.

(ⅰ) $r-cq > 0$ implies that $\left(0, c\right)$ is a hyperbolic saddle point.

(ⅱ) $r-cq = 0$ implies that $\left(0, c\right)$ is a non-hyperbolic saddle point.

(ⅲ) $r-cq < 0$ it has that $\left(0, c\right)$ is a hyperbolic attractor.

Theorem 9. When the equilibrium point or singularity $\left(x_{e}, nx_{e}+c\right)$

(i) belongs to the first quadrant, then is locally asymptotically stable (l.a.s).

(ii) coincides with the point $\left(0, c\right)$ is a non-hyperbolic attractor.

Proof. The Jacobian matrix evaluated at point $\left(x, nx+c\right)$ is

Then,

Replacing $x$ by $x_{e} = \frac{\left(r-cq\right) K}{r+nqK}$ it has

which can be positive, negative or zero.

Besides, tr$DG_{\sigma }\left(x, nx+c\right) =$ $\frac{ Kr-Ks-2rx-Kcq-Knqx}{K}$.

It has that, tr$DG_{\sigma }\left(x_{e}, nx_{e}+c\right) = -\frac{ rs+r^{2}-cqr+Knqs}{r+Knq} < 0$, if and only if, $c < \frac{r}{q}$.

Hence, $\left(x_{e}, nx_{e}+c\right)$ is locally asymptotically stable, for $c < \frac{r}{q}$.

If $r-cq = 0$ implies that $\left(x_{e}, nx_{e}+c\right)$ conincides with $\left(0, c\right)$ is a non-hyperbolic attractor.

Remark 10. Considering the Poincaré-Bendixson Theorem ^[8,39], if the equilibrium point $\left(x_{e}, nx_{e}+c\right)$ is locally asymptotically stable (l.a.s.) in the bidimensional system $\left(1\right)$, then it is globally asymptotically stable (g.a.s.). This result is reinforced by building a suitable Lyapunov function ^[36,40] and also proving that there are no cycles ^[12], which is shown in the following theorems.

Theorem 11. The equilibrium point $\left(x_{e}, nx_{e}+c\right)$ is globally asymptotically stable (g.a.s.).

Proof. Because $\left(1\right)$ is a system of Kolmogorov type describing a continuous-time predator-prey model, we will consider the Lyapunov function ^[38] described by

which is similar to a proposal in ^[36] and in ^[40], with $c_{1}$ and $c_{2} > 0$, to be determined.

Clearly $V\left(x, y\right) > 0$, for all $\left(x, y\right)$ $\in Int\left(\mathbb{R} _{0}^{+}\right)$ and $x\neq x_{e}$ e $y\neq y_{e}$ ^[36].

Deriving the function $V$ with respect to $t$, we have:

Factoring for the respective derivatives, we get:

Replacing it has

or else,

Remembering that$\left(x_{e}, y_{e}\right)$ satisfies the equations of both isoclines we have

Reordered we become

Expanding the products we have

Choosing $c_{2} = \frac{Kq}{sn}c_{1}$, finally it is obtained that

Thus $V\left(x, y\right)$ is a suitable Lyapunov function ^[36,38,40], since it is defined negative.

Then the equilibrium point $\left(x_{e}, nx_{e}+c\right)$ is g.a.s, by Lyapunov's stability Theorem ^[38].

Theorem 12. There are not periodic solutions.

Proof. We will use the Bendixson-Dulac criterion ^[29] considering the function $g\left(x, y\right) = \frac{1}{xy} > 0$ for all $\left(x, y\right)$ at inside the first quadrant, $Int\left(\mathbb{R} _{0}^{+}\right) ^{2}$.

Considering the sum

where $M\left(x, y\right) = \left(r\left(1-\frac{x}{K}\right) -qy\right) x$ and $N\left(x, y\right) = s\left(1-\frac{y}{nx+c}\right) y$,

it gets that

Then, $S\left(x, y\right) < 0$ for all $\left(x, y\right)$ at inside the first quadrant.

Therefore, by the Bendixson-Dulac criterion ^[8,29,39], there are nonconstant periodic solutions (or limit cycles) within the first quadrant.

3.   Model with CAP and generalist predator

The second model to study is described by the Kolmogorov-type system ^[14,36]:

where $x = x\left(t\right)$, $y = y\left(t\right)$, are the population sizes. The parameter $0 < \beta < 1$, indicates the level of interference between predators and the other parameters have the same meanings given above in Table 1. Therefore, $\lambda =$ $\left(r, K, q, s, n, c, \beta \right) \in \mathbb{R}_{+}^{6}\times \left] 0, 1\right[$.

The equilibrium points are $\left(0, 0\right)$, $\left(K, 0\right)$, $\left(0, c\right)$ and $\left(x_{e}, nx_{e}+c\right)$, where $x_{e}$ satisfies the isoclines equation

The graphic of $g_{1}$ represents a decreasing curve that passes through the points $\left(0, \left(\frac{r}{q}\right) ^{\frac{1}{\beta }}\right)$ and $\left(K, 0\right)$.

The graph of $g_{2}$ represents an increasing straight line through the point $\left(0, c\right)$.

Then the isoclines intersect, if and only if, $c < \left(\frac{r}{q}\right) ^{ \frac{1}{\beta }}$.

To determine the nature of the hyperbolic equilibrium points we require the Jacobian matrix, which is

but it is not defined for points on the $x-axis$, i.e., when $y = 0$. Therefore, we must use an alternative method to obtain the type of stability of the points $\left(0, 0\right)$ and $\left(K, 0\right)$.

In order to simplify the calculations, it is convenient to reduce system $\left(2\right)$ to a normal form; we will follow the methodology used in ^{[15,17,42,43]} making a change of variable and a time rescaling, obtaining the following:

Theorem 13. System $\left(2\right)$ is topologically equivalent to the Kolmogorov-type system ^[14,36]:

where $\eta = \left(Q, S, C, \beta \right) \in \mathbb{R}_{+}^{3}\times ]0, 1[$. Its domain is

with $Q = \frac{q\left(nK\right) ^{\beta }}{r}$, $S = \frac{s}{r}$ y $C = \frac{c }{nK}$.

Proof. Let us $x = Ku$ and $y = nKv$; hence $\frac{dx}{dt} = K\frac{du}{dt}$ and $\frac{dy }{dt} = nK\frac{dv}{dt}$.

Replacing and factoring in system $\left(2\right)$ we have:

Making a time rescaling (independent variable), given by $\tau = \frac{r}{u}t$; by the chain rule we get

Replacing in the vector field$V_{\lambda }\left(x, y\right)$ it has

Defining $Q = \frac{q\left(nK\right) ^{\beta }}{r}$, $S = \frac{s}{r}$ and $C = \frac{c}{nK}$, we get the system $(3)$.

Remark 14. We have built the function $\varphi :\overline{\Omega }\times \mathbb{R} \longrightarrow {\Omega \times \mathbb{R}}$ given by

The Jacobian matrix of transformation $\varphi$ is

and det$D\varphi (u, v, \tau) = \frac{nK^{2}}{r}\left(u+\frac{c}{nK}\right) > 0$.

Then, $\varphi$ is a diffeomorphism ^[8,9] preserving the time orientation, whereby the vector field $X_{\lambda }$ is topologically equivalent to the vector field $U_{\eta }\left(u, v\right) = \varphi \circ X_{\lambda }\left(x, y\right)$.

In this way, system $\left(3\right)$ is a continuous extension of system $\left(2\right)$.

The equilibrium points of the system $\left(3\right)$ are $\left(0, 0\right)$, $\left(1, 0\right)$, $\left(0, C\right)$ and $\left(u_{e}, u_{e}+C\right)$, with $u_{e}$ determined by the intersection of the isoclines

The Jacobian matrix of system $\left(3\right)$ is

with,

according to the derivative of the triple product. We note that $DU_{\eta }\left(u, v\right)$ is not defined at the points $\left(0, 0\right)$ and $\left(1, 0\right)$, due to the entry $DU_{\eta }\left(u, v\right) _{12} = -Q\beta uv^{\beta -1}\left(u+C\right)$ is not defined for $v = 0$.

3.1.   Main results

The main properties of the system $(3)$ or vector field $U_{\eta }\left(u, v\right)$ and therefore of the modified Leslie-Gower model described by the system $\left(2\right)$ are:

Lemma 15. The set

is a positively invariant region.

Proof. Clearly, coordinate axes are invariant sets

Considering $u = 1$, we obtain that $\frac{dx}{dt} = -Q\left(1+C\right) v^{\beta }$.

Therefore, solutions with an initial condition outside the set $\bar{\Gamma}$ will enter to the set $\bar{\Gamma}$ and the solutions within the set will not leave it, whatever the sign of the equation $\frac{dv}{d\tau } = S\left(1+C-v\right) v$.

Remark 16. The set $\Gamma = \left\{ \left(x, y\right) \in \mathbb{R} ^{2}:0\leq x\leq K\text{, }y\geq 0\right\}$ is a positively invariant region $\left(2\right)$.

The subregion

is also a positively invariant region for the system $\left(3\right)$.

Lemma 17. The solutions are bounded.

Proof. Using the Poincaré compactification ^[39].

Let us $X = \frac{u}{v}$ and $Y = \frac{1}{v}$; clearly, if $v\rightarrow \infty$, the point $\left(0, \infty \right)$ is associated with the point $\left(0, 0\right)$.

Thus, $u = \frac{X}{Y}$ and $v = \frac{1}{Y}$, obtaining

Replacing, it get

Factoring by $\frac{1}{Y}$ in each equation becomes

The latter system is not defined for $Y = 0.$ To solve this difficulty, a rescaling of the time given by $T = \frac{1}{Y^{2}}\tau$ is considered. So by chain rule, it has

Applying these relationships to the vector field $\bar{U}_{\eta }\left(X, Y\right)$ we obtain

Clearly, the Jacobian matrix evaluated at point $\left(0, 0\right)$ is not defined since the derivative of the term $QXY^{1-\beta }$ respect to Y is $\frac{d}{dY}\left(QXY^{1-\beta }\right) = Q\frac{X}{Y^{\beta }} \left(1-\beta \right)$.

To overcome this difficulty, we will use the method of the blowing-up directional with the change of variables $X = pq^{2}$, $Y = q$. Hence, $\frac{dX }{dT} = q^{2}\frac{dp}{dT}+2pq\frac{dq}{dT}$, $\frac{dY}{dT} = \frac{dq}{dT}$; thus, $\frac{dp}{dT} = \frac{1}{q^{2}}\left(\frac{dX}{dT}-2pq\frac{dq}{dT} \right)$.

By replacing and factoring the new system is obtained

whose Jacobian matrix evaluated at point $\left(0, 0\right)$ is the null matrix.

Rescaling now by changing $\kappa = qT$, we get

The Jacobian matrix of the vector field $\check{U}_{\eta }\left(p, q\right)$ is

with

and

Then, $D\check{U}_{\eta }(0, 0) = \left(\begin{array}{cc} -S\ \ & \ \ 0 \\ 0\ \ & \ \ S \end{array} \right)$.

Therefore, the point $(0, 0)$ of the vector field $\check{U}_{\eta }\left(p, q\right)$ is a hyperbolic saddle, attractor by the $x-axis$ and repeller by the $y-axis$. Then there is a non-hyperbolic chair of the vector fields $\hat{U}_{\eta }\left(p, q\right)$ and $\bar{U}_{\eta }\left(X, Y\right)$.

So the point $\left(0, \infty \right)$ is a non-hyperbolic saddle of the compacted vector field of $U_{\eta }\left(u, v\right)$.

Thus, the trajectories are bounded.

Remark 18. As in the previous model, this property ensures that the model is well-posed. ^[1].

As in the model without CAP, for the existence of the equilibrium point $\left(u_{e}, u_{e}+C\right)$ we have:

Lemma 19. The equilibrium point $\left(u_{e}, u_{e}+C\right)$

(i) belongs to the first quadrant $Int\left(\mathbb{R} _{0}^{+}\right)$ and is unique, if and only if, $CQ^{\beta } < 1$,

(ii) coincides with the point $\left(0, C\right)$, if and only if, $CQ^{\beta } = 1$,

(iii) is outside the first quadrant, if and only if, $CQ^{\beta } > 1$.

Proof. The equilibrium point $\left(u_{e}, u_{e}+C\right)$ is at the intersection of the isoclines and the abscissa $x_{e}$ satisfies the equations

We have that:

$\qquad g_{1}$ is a decreasing function and passes through the points $\left(0, \left(\frac{1}{Q}\right) ^{\beta }\right)$ y $\left(1, 0\right)$,

$g_{2}$ is an increasing linear function that passes through the point $\left(0, C\right)$.

Therefore, the graphics of $g_{1}\left(u\right)$ and $g_{2}\left(u\right)$ intersect inside the first quadrant $Int\left(\mathbb{R} _{0}^{+}\right) ^{2}$, if and only if, $C < \left(\frac{1}{Q}\right) ^{\beta }$.

Since the Jacobian matrix is not defined when $v = 0$, it considers the change of variable given by $w = v^{\beta }$.

Thus, $\frac{dv}{d\tau } =$ $\frac{1}{\beta }w^{\frac{1}{\beta }-1}\frac{dw}{ d\tau }$. Replacing the second equation of the system $\left(2\right)$ we obtain the new system

with $\eta = \left(Q, S, C, \beta \right) \in \mathbb{R}_{+}^{3}\times ]0, 1[$. The equilibrium points are $\left(0, 0\right)$, $\left(1, 0\right)$, $\left(0, C^{\beta }\right)$ and $\left(u_{e}^{\ast }, \left(u_{e}^{\ast }+C\right) ^{\beta }\right)$ with $u_{e}$ determined by the intersection of the isoclines

or else, $u_{e}$ satisfies the equation

The isocline $h_{2}\left(u\right)$ is a straight line that passes through the points $\left(0, \frac{1}{Q}\right)$ and $\left(1.0\right)$, while the isocline $h_{1}\left(u\right)$ is an increasing curve concave downwards, passing through the point $\left(0, C^{\beta }\right)$. See Figure 2.

Clearly, there is an intersection between both curves, if and only if, $C^{\beta } < \frac{1}{Q}$, i.e., if and only if, $QC^{\beta } < 1$.

The Jacobian matrix of system $\left(4\right)$ is:

with $DW_{\sigma }\left(u, w\right) _{11} = u\left(u+C\right) \left( -1\right) +\left(1-u-Qw\right) \left(u+C\right) +\left(1-u-Qw\right) u$.

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

Proof. The Jacobian matrix evaluated in $\left(0, 0\right)$ is:

Then, det$DW_{\eta }\left(0, 0\right) = \beta C^{2}S > 0$ and tr$DW_{\eta }\left(0, 0\right) = C\left(1+\beta S\right) > 0.$

Therefore, keeping in mind the trace and the determinant Theorem ^[7,37], the equilibrium point $\left(0, 0\right)$ is repeller.

Lemma 21. The equilibrium point $\left(1, 0\right)$ is a hyperbolic saddle.

Proof. The Jacobian matrix evaluated in $\left(1, 0\right)$ is:

hence, det$DW_{\eta }\left(1, 0\right) = -\beta \left(1+C\right) ^{2}S < 0$.

Then, considering the trace and the determinant Theorem ^[7,37], the equilibrium point $\left(1, 0\right)$ is a hyperbolic saddle.

Remark 22. With the two previous lemmas and bearing in mind the change of variable made to obtain the system $\left(4\right)$ or vector field $W_{\eta }\left(u, w\right)$, the equilibrium points $\left(0, 0\right)$ and $\left(1, 0\right)$, of the system $\left(3\right)$ or vector field $U_{\eta }\left(u, v\right)$, are non-hyperbolic repeller and non-hyperbolic saddle, respectively. Therefore, the equilibrium points $\left(0, 0\right)$ and $\left(K, 0\right)$ of the system $\left(2\right)$ are also non-hyperbolic repeller and non-hyperbolic saddle, respectively.

To determine the nature of the equilibrium points $\left(0, C\right)$ and $\left(u_{e}, u_{e}+C\right)$ we will again consider the system $\left(3\right)$.

Lemma 23. The equilibrium point $\left(0, C\right)$ is:

i) a hyperbolic saddle, if and only if, $QC^{\beta } < 1$,

ii) a non-hyperbolic saddle, if and only if, $QC^{\beta } = 1$,

iii) an attractor node, if and only if, $QC^{\beta } > 1$,

Proof. The Jacobian matrix evaluated in $\left(0, C\right)$ is:

As det$DW_{\eta }\left(0, C\right) = -SC^{2}\left(1-QC^{\beta }\right)$ we have that

(ⅰ) if $QC^{\beta } < 1$, the equilibrium point $\left(0, C\right)$ is saddle,

(ⅱ) if $QC^{\beta } = 1$, the equilibrium point $\left(0, C\right)$ is non-hyperbolic saddle,

(ⅲ) If $QC^{\beta } > 1$, thus

So $\left(0, C\right)$ is an attractor node.

Remembering that If $QC^{\beta } > 1$, the positive equilibrium is out of the first quadrant, by the Poincaré-Bendixson Theorem, the point $\left(0, C\right)$ is globally asymptotically stable (g.a.s.).

Remark 24. When $C = 0$, the equilibria $\left(0, C\right)$ and $\left(0, 0\right)$ collapse. Then, for certain conditions on the parameter values, the point is a non-hyperbolic attractor. So, there exist trajectories in the system $\left(1\right)$, with the point $\left(0, 0\right)$ as their $\omega -limit$ (See Figure 3).

Theorem 25. When the equilibrium point $\left(u_{e}, u_{e}+C\right)$ belongs to $Int\left(\mathbb{R} _{0}^{2}\right)$, is a hyperbolic node locally asymptotically stable (l.a.e.).

Proof. The Jacobian matrix of system $\left(3\right)$ evaluated in $\left(u, u+C\right)$ is:

Then,

and

Remark 26. As the unique positive equilibrium point $\left(u_{e}, u_{e}+C\right) \in Int\left(\mathbb{R} _{0}^{2}\right)$ is l.a.e., if and only if, $QC^{\beta } < 1$, Considering the Poincaré-Bendixson Theorem, we have that it is globally asymptotically stable (g.a.s.). This result is reinforced by the following theorems.

Theorem 27. The equilibrium point $\left(u_{e}, u_{e}+C\right)$ is globally asymptotically stable (g.a.s).

Proof. Consider the function described by

with $c_{1}$ and $c_{2} > 0$, to be determined.

Clearly $V\left(u, v\right) > 0$, for all $\left(u, v\right)$ $\in Int\left(\mathbb{R} _{0}^{+}\right) ^{2}$ and $u\neq u_{e}$, $v\neq v_{e}$ ^[36].

Deriving the function $V$ with respect to $\tau$, it has:

Replacing $\frac{1}{u}\frac{du}{d\tau }$ and $\frac{1}{v}\frac{dv}{d\tau }$ is obtained

i.e.,

As $\left(u_{e}, u_{e}+C\right)$ satisfies the equations of both isoclines we have

Reordering it becomes

Expanding the products we have

Choosing, $c_{2} = \frac{Q}{\beta S}c_{1}$ we obtain that

but,

with $p\left(v, v_{e}\right)$ a positive function.

Then,

Therefore, $V\left(u, v\right)$ is a suitable Lyapunov function ^[36,38,40], since it is defined negative.

Then the point $\left(u_{e}, u_{e}+C\right)$ is g.a.s, by the Lyapunov stability Theorem ^[36,38].

Theorem 28. There are no cycles.

Proof. Using the Bendixson-Dulac criterion with the function $g\left(u, v\right) = \frac{1}{u\left(u+C\right) v} > 0$, for all $\left(u, v\right)$ inside the first quadrant, $Int\left(\mathbb{R} _{0}^{+}\right)$.

Considering the sum

where $M\left(u, v\right) = \left(1-u-Qv^{\beta }\right) u\left(u+C\right)$ y $N\left(u, v\right) = S\left(u+C-v\right) v$,

it is obtained that

namely

Hence $S\left(u, v\right) < 0$ for all $\left(u, v\right)$ inside the first quadrant.

Therefore, by the Bendixson-Dulac criterion ^[8,39], there are no cycles (or limit cycles) within the first quadrant.

4.   Some numerical simulations

In Figure 3, considering values for system $\left(3\right)$, we can see that the point $\left(0, 0\right)$ of system $\left(1\right)$ is a global attractor for the trajectories.

The following three figures show the global stability of the positive equilibrium point (Figures 4, 5 and 6), when exists at interior of the first quadrant $Int\left(\mathbb{R} _{0}^{+}\right)$ and the global stability of the equilibrium point $\left(0, C\right)$.

5.   Conclusion

In this work, the study of the dynamics of a predator model prey of the Leslie-Gower type was carried out considering the competition (or interference) among predators, CAP, described by the system $\left(1\right)$.

We have also detailed the main properties of the modified Leslie-Gower model considering a generalist predator in order to carry out a comparative study of the influence that the CAP has on the original model.

To simplify the calculations in the second model, a reparameterization and a time rescaling were carried out and the system $\left(3\right)$ topologically equivalent to $\left(2\right)$ was obtained.

The main features of the modified Leslie-Gower type models with CAP are:

● There is a region of positive invariance

● The solutions are bounded, that is, the system is well-posed ^[1].

● The model with CAP is not defined in $\left(0, 0\right)$ and neither the Jacobian matrix of the topologically equivalent system $\left(2\right)$ is defined for $y = 0$. The nature of the equilibrium points on the $x-axis$ was determined using system $\left(3\right)$, topologically equivalent to system $\left(2\right)$.

● The equilibrium point $\left(K, 0\right)$ is a non-hyperbolic saddle.

● The equilibrium point $\left(0, 0\right)$ is a non-hyperbolic repeller.

● There are no cycles (or limit cycles) inside the first quadrant $Int\left(\mathbb{R} _{0}^{+}\right) ^{2}$.

Therefore, both models, with or without CAP, have similar dynamics, since in both the common equilibrium points have the same characteristics (nature), except for the hyperbolicity of the equilibrium points $\left(0, 0\right)$ and $\left(K, 0\right)$ in each of the models.

The main difference for the model with generalist predators is the existence of the equilibrium $\left(0, C\right)$ over the vertical axis. Then, there are conditions ensuring that this point is globally asymptotically stable (g.a.s.). This fact is obtained for $CQ^{\beta } > 1$, with $\left(Q, S, C, \beta \right) \in \mathbb{R}_{+}^{2}\times ]0, 1[$, implying that the CAP can produce the extinction of the prey population and the permanence of the predator populations in its maximum environmental carrying capacity.

Both in system $\left(2\right)$ and in the Leslie-Gower model without interference among predators described by system $\left(1\right)$, there is a unique point of coexistence of species, which is globally asymptotically stable (gas) when it exists.

This means that whatever the initial size of each of the populations, as time passes, the population sizes stabilize at that point of equilibrium. This property of both models is highly desirable by some ecologists or managers of exploitation of renewable resources. For example, if we consider that a fish population is caught by many stakeholders (artisanal or industrial fishing, that is, humans acting as predators.

These conclusions are deduced in our analysis, assuming the CAP described by the function $B\left(x, y\right) = qxy^{\beta }$, with $0 < \beta < 1$, but other outcomes would be obtained considering distinct prey-dependent functional response, i.e., assuming $B\left(x, y\right) = h\left(x\right) y^{\beta }$, with $0 < \beta < 1$ and $h\left(x\right)$ a Holling type functional response.

Although similar results have been obtained for the Leslie-Gower model with or without CAP ^[19], we cannot conclude that CAP induces the stability of the predator-prey interaction. We know that in many of the real-world interactions, there are extinctions of one or both populations ^[30], or there exist periodic oscillations of the population sizes, which can be explained by the existence of limit cycles ^[12,15,43].

These different situations are obtained as possible dynamics appearing in models that consider the CAP described by other mathematical forms. Such cases are obtained, for example, by assuming a Beddington-DeAngelis functional response ^[3,26] in the Leslie-Gower model ^[27] or else, assuming the addition of a quadratic term in the predator growth equation ^[5,28], although this second form is mainly included in Gause type predation models ^[5,28].

In the model studied in ^[27], the positive equilibrium point is not necessarily g.a.s., when it exists, and there are also limit cycles around a positive equilibrium point, which implies oscillations of population sizes ^[27]. For the same set of parameters, different dynamic behaviors coexist, due to the appearance of a separatrix curve in the phase plane ^[27]. Therefore, the trajectories or solutions of the system are highly sensitive to the initial conditions, since there exist trajectories with very near initial conditions, but to distinct side of that separatrix curve having far $\omega -limit$.

In this way, an important dilemma arises, which is: What is the most efficient model to represent a certain interaction of predation in nature?. We believe that the same interaction can be represented by a given model (whatever its mathematical nature) at a defining moment or time, but its validity will depend strongly on the underlying hypotheses, which are usually not declared.

Bearing in mind that the dynamic complexity of an ODE system does not imply the ecological complexity of the model described by that system, we could not reliably affirm that one model is better than another. However, more complex ecological situations (for example, considering gestation time, non-homogeneous distribution of populations or classified by age and sex, etc.) require the use of more sophisticated mathematical tools than ODEs.

Conflict of interest

The authors declare no conflicts of interest in this paper.