Purinergic P2X7R expressed on regulatory T cells potentially links molecular mimicry to autoimmune responses

1.   Introduction

The predator-prey model is a model in biomathematics used, and modified, by many researchers to describe the dynamics between a prey and predator species in order to determine conditions in their parameters that allow for either the conservation or extinction of one or both species, where such models are generally posed by systems of ordinary ^[1,2,3,4,5] or partial differential equations ^[6,7,8,9,10].

In particular, the Leslie-Gower models assumes that the intrinsic growth of both species is described by logistic functions, whose carrying capacity between prey and predators is subject to the capacity of the environment or directly proportional to the population size of the prey, respectively ^[11]. Similarly, this model assumes a functional response of the predator that is generally described by Holling-type functions ^{[12,13,14,15,16]}.

For example, assuming that the predator response functional is given by a Holling Ⅱ function $f(x) = \frac{ax}{x+d}$ ^[17,18,19], where $a > 0$ is the maximum per capita consumption rate of predators, and $0 < d < K$ is the semi-saturating rate of capture, where predators can be harvested ^[20,21], the dynamics of population sizes of prey $x(t)\geq 0$ and predators $y(t)\geq 0$ is given by

where $K > 0$ is the carrying capacity of prey in the environment, which normally determines the resources available for prey survival, $n > 0$ is the quantity of prey consumed by the predator and converted into new births, $q > 0$ is the catchability coefficient, $E > 0$ is the harvesting efforty, and $r, s > 0$ are the intrinsic birth rates of prey and predators, respectively. In particular, if the diet of predators is based exclusively on prey consumption, then $c = 0$ ^[22,23,24]; otherwise, $c > 0$ describes the alternative food for predators, so their diet is not based exclusively on prey consumption ^[25,26].

However, given the existence of prey that take refuge from the predator to subsist in the environment ^{[1,29,30,31,32]}, and whose intrinsic growth is subject to an auxiliary parameter $b > 0$ (half-saturation) that affects the hyperbolic function of the per capita birth rate curve of the prey ^[27,28], in addition to the maximum birth rate $r > 0$, then model (1.1) is modified by

where $P < K$ is the population size of refugia in prey. In particular, the intrinsic growth of prey approximates logistic growth for a large $x$; however, if $x$ is small, the parameter $b > 0$ influences its dynamics, of which it may represent examples such as inbreeding, ease of reproduction, difficulty in finding mates, or mating resistance at low temperatures.

Observe that the model (1.2) is used if the initial prey population size is above the critical refuge population size, that is, if $P < x(0)$. For the case where $x(0) < P$, the entire prey population size is refugeed from the predator whenever $x < P$, so there is no interaction between the two species and, hence, the functional response of the predator is deactivated when $x < P$. In this case, the dynamics of the prey is described by a logistic differential equation $\dot{x} = rx\left(1-\frac{x}{K}\right)$ for the case in which there are either no external factors inhibiting its growth, or by an autonomous differential equation $\dot{x} = \frac{rx}{x+b}\left(1-\frac{x}{K}\right)$, until they exceed their critical refuge population size and the quantity $x-P$ becomes susceptible to predator consumption.

On the other hand, if the predator's diet is based only on the consumption of prey, they are forced to change their diet due to the deficit of prey in order to survive in the environment ^[33,34,35], so the model (1.2) should be modified by

with

However, if the predator's diet is not based exclusively on prey consumption, but has equal probability of consumption with the alternative food provided in the environment, the predator's diet is not affected and depends only on the alternative food provided by the environment in the absence of prey. In this case, the model (1.2) is modified by

with $\epsilon$ as described in (1.4).

Consequently, and given the possible scenarios for predator growth and survival at low and high prey densities subject to their refuge, the objective of this paper is to perform a qualitative and bifurcation analyses to models (1.3) and (1.5), as shown in Sections 2 and 3, respectively, to determine similarities or difference in both models.

2.   Prey refuge without alternative food for predators above threshold value

Let $x(t)\geq 0$ and $y(t)\geq 0$ be considered as the population sizes of prey and predators, respectively, whose dynamics are given by the model (1.3), defined in the biological sense region

with $M = \mbox{max}\{ n(K-P), c \}$. In particular, the vector fields of the model (1.3) are given by

where $X$ and $Y$ are defined in $\Sigma^{+} = \left\lbrace (x, y)\in \Omega: P < x\leq K \right\rbrace$ and $\Sigma^{-} = \{(x, y)\in \Omega: 0 < x < P \}$, respectively.

Since the model (2.2) is equivalent to a Filippov system ^{[36,37,38,39]}, if $f(x, y) = x-P$, $Xf(p) = \langle X(p), \mathrm{grad}f(p)\rangle$ and $Yf(p) = \langle Y(p), \mathrm{grad}f(p)\rangle$, then $Xf(p) = Yf(p) = \frac{rP}{P+b}\left(1-\frac{P}{K}\right) > 0$ for all $p \in \Sigma = \{(x, y)\in \Omega: x = P \}$. Therefore, since the vector field $Y$ is composed of a system of autonomous and independent differential equations, the trajectories $\varphi_{Z_{0}}$ with initial condition $(x(0), y(0)) \in \Sigma^{-}$ must cross $\Sigma$ and remain in $\Sigma^{+}$. Additionally, $\varphi_{Z_{0}}$ with the initial condition $p = (P, y(0)) \in \Sigma$ is defined as $\varphi_{Z_{0}}(p, t) = \varphi_{X}(p, t)$ for all $t \in I \cap \{t > 0\}$ and $\varphi_{Z_{0}}(p, t) = \varphi_{Y}(p, t)$ for $t \in I \cap \{t \leq 0\}$.

Lemma 1. For an arbitrary initial condition $(x(0), y(0)) \in \Omega$, the model (2.2) has a unique $\varphi_{Z_{0}}$ trajectory and remains in $\Omega$.

Proof. Since the vector fields $X$ and $Y$ are continuously differential and the trajectories $\varphi_{Z_{0}}$ cross $\Sigma$ or remain $\Sigma^{+}$ for all initial conditions $(x(0), y(0)) \in \Sigma^{-}$ or $(x(0), y(0)) \in \Sigma^{+}$, by the existence and uniqueness theorem ^[40] for every vector field the uniqueness of the trajectories are guaranteed. Additionally, for all $0\leq x < P$, if $x = 0$ then $\dot{x} = 0$ for all $y\geq 0$. Similarly, if $y = 0$, then $\dot{y} = 0$, and if $y = M$, we have that $\dot{y} \leq 0$. Thus, the trajectories $\varphi_{Z_{0}}$ do not cross $\Sigma^{-}$. Moreover, for all $P < x\leq K$, if $x = K$, then $\dot{x} = 0$; if $y = 0$, then $\dot{y} = 0$; and if $y = n(K-L)+c$, then $\dot{y}\leq 0$, so $\varphi_{Z_{0}}$ does not cross $\Sigma^{+}$. The dynamics of the model at points $(P, 0)$ and $(P, M)$ do not cross $\Omega$ and remain to be analyzed. In this case, since the trajectories $\varphi_{Z_{0}}$ with the initial condition $(x(0), y(0)) \in \Sigma^{-} \cup \Sigma^{+}$ cross $\Sigma$, we have that the trajectories $\varphi_{Z_{0}}$ do not escape from the points $(P, 0)$ and $(P, M)$. □

2.1.   Existence and local stability at equilibrium points

The equilibria in the model (2.2) over the coordinate axes are given by $P_{0} = (0, 0) \in \Sigma^{-}$, $P_{1} = (K, 0) \in \Sigma^{+}$, and $P_{2} = \left(0, \frac{c(s-qE)}{s}\right) \in \Sigma^{-}$ if $s-qE > 0$. The local stability of the equilibria is summarized in the following result.

Lemma 2. If $s-qE < 0$, then $P_{0}$ and $P_{1}$ are saddle points. If $s-qE > 0$, then $P_{0}$ is a locally unstable node, $P_{1}$ a locally stable node and $P_{2}$ is a saddle point.

Proof. The local stability of $P_{1} \in \Sigma^{+}$ and $P_{0}, P_{2} \in \Sigma^{-}$ are guaranteed from the computation of the eigenvalues of the Jacobian matrix $DX$ of the vector field $X$ computed in $P_{1} \in \Sigma^{+}$ and of the Jacobian matrix $DY$ of the vector field $Y$ computed in $P_{0}, P_{2} \in \Sigma^{-}$. □

On the other hand, since the vector field $Y$ has no interior equilibria, the intersections in $\Sigma^{+}$ of nullclines

are defined in the vector field $X$, and determine the number of interior equilibria in the model (2.2). In particular, note that the model (2.2) has no interior equilibria if $s-qE < 0$ since the nullclines (2.3) do not intersect in $\Sigma^{+}$.

For $s-qE > 0$, the possible positive roots of the polynomial

with

determine the number of interior equilibria in the model (2.2).

The following result shows the existence of a unique interior equilibrium $P_{3} = (x_{*}, y_{*})$ in the model (2.2), with $x_{*}$ as the positive root of the polynomial (2.5) and $y_{*} = f_{1}(x_{*})$, as observed in Figure 1.

Figure 1.  Nullclines (2.3) with $s-qE > 0$. The green and margenta curves represent the nullclines $y = f_{1}(x)$ and $y = f_{2}(x)$, respectively. The brown and black dotted lines are $x = P-d$ and $x = P$, respectively.

DownLoad: Full-Size Img PowerPoint

Lemma 3. If $s-qE > 0$, the model (2.2) has a unique interior equilibrium $P_{3}$.

Proof. By Descartes' sign criterion, we have that the polynomial $Q(x)$ (2.4) must have two positive roots: $x_{*}$ and $x_{+}$. However, since $y = f_{1}(x)$ is a straight line with a positive slope and intercept at $(P, 0)$, $f_{2}(x) \to -\infty$ when $x \to \pm \infty$ and $f_{2}(x) \to \pm\infty$ for $x \to P^{\pm}$, and by constraint, we have that $x_{+} < P$ and $P < x_{*}$. Therefore, the model (2.2) has a unique interior equilibrium. □

The local dynamics of the interior equilibrium $P_{3}$ remain to be analyzed. In this case, since $\Omega$ is invariant and the equilibria $P_{0}$, $P_{1}$, and $P_{2}$ are unstable, by the Poincare - Bendixson Theorem ^[41,42] we have that $P_{3}$ cannot be a local saddle point, that is, $\mbox{det }DX(P_{3}) > 0$, with $DX$ the Jacobian matrix of the vector field $X$ computed in $P_{3}$.

On the other hand, the sign of the trace in the Jacobian matrix $DX$ of the vector field $X$ computed in $P_{3}$, with

which determines the local stability of the equilibrium $P_{3}$, and is summarized in the following result. In particular, if $P_{3}$ is locally unstable, by the Poincare -Bendixson Theorem it guarantees the existence of at least one stable $\Gamma_{X}$ limit cycle in $\Sigma^{+}$, which does not collide with $\Sigma$.

Lemma 4. For $s-qE > 0$ and $\triangle: = [{{tr}}DX(P_{3})]^{2}-4{det}DX(P_{3}) < 0$, then $P_{3}$ is either a locally stable or unstable focus if ${{tr}}DX(P_{3}) < 0$ or $0 < {{tr}}DX(P_{3})$, respectively. For $s-qE > 0$ and $\triangle: = [{{tr}}DX(P_{3})]^{2}-4{det}DX(P_{3}) > 0$, then $P_{3}$ is either a locally stable or unstable node if ${{tr}}DX(P_{3}) < 0$ or $0 < {{tr}}DX(P_{3})$, respectively.

2.2.   Global stability at equilibrium points

In the first instance, the following result shows conditions on the model parameters (2.2) for which guarantees the global stability of $P_{1}$.

Theorem 1. If $s-qE < 0$, then $P_{1}$ is globally asymptotically stable in $\Omega$.

Proof. Since $P_{0}$ is a saddle point, $P_{2}, P_{3} \in \Omega$ and $\dot{y} = (s-qE)y-\frac{sy^{2}}{n(x-P)} < 0$ for all $(x, y) \in \Sigma^{+}$, then the model (2.2) has no limit cycles in $\Sigma^{+}$. Furthermore, since the model (2.2) has no interior equilibria and no limit cycles in $\Sigma^{-}$, and the differential equations describing the vector field $Y$ are autonomous and independent, all trajectories $\varphi_{Z_{0}}$ with initial condition $(x(0), y(0)) \in \Sigma^{-}$ cross $\Sigma$ and remain in $\Sigma^{+}$, so the trajectories $\varphi_{Z_{0}}$ converge to $P_{1}$. □

On the other hand, by using the Bendixon-Dilac ^[40] criterion on the model (2.2), the following result shows a first sufficient condition on the model parameters (2.2) for which guarantees the global stability of $P_{3}$.

Theorem 2. If $s-qE > 0$ and $b(b+P) < (b+K)(d-P)$, then $P_{3}$ is globally asymptotically stable in $\Omega$.

Proof. If

are functions describing the vector field $X$ in $\Sigma^{+}$, and

is a Dirac function, then

with $\alpha = b(d-K-P)+K(d-P)-b^2 > 0$ for all $(x, y) \in \Sigma^{+}$, so the model (2.2) has no limit cycles in $\Sigma^{+}$. Therefore, and analogous to the proof of Theorem 1, we have that $P_{3}$ is globally asymptotically stable in $\Omega$. □

From the Lyapunov criterion ^[40], the following result shows a second sufficient condition on the parameters of the model (2.2) for which the global stability of $P_{3}$ is guaranteed in $\Omega$.

Theorem 3. If $s-qE > 0$ and $d+2P > K$, the equilibrium $P_{3}$ is globally asymptotically stable in $\Omega$.

Proof. Similarly to what we proved in Theorem 1, it is sufficient to prove that $P_{3}$ is globally asymptotically stable in $\Sigma^{+}$. Indeed, if

for all $(x, y) \in \Sigma^{+}$ and $\rho(x) = \frac{a(x-P)}{(x-P)+d}$ is a Lyapunov function, $g(x) = \frac{\frac{rx}{x+b}\left(1-\frac{x}{K}\right)}{\rho(x)}$ and (2.3), we have that

Since $P < x$, a sufficient condition for $\dot{V} < 0$ is to prove that $g(x)$ is non-increasing. Indeed, by substituting $\bar{x} = x-P$, then

and so

with $\beta = d-K+2P > 0$. □

2.3.   Bifurcation analysis

In the first instance, the model (2.2) has a Hopf bifurcation if $\mbox{tr }DX(P_{3}) = 0$, as observed in the following result. In particular, and as observed in Lemma 4, if $\mbox{tr }DX(P_{3}) < 0$, then $P_{3}$ is unstable with at least one stable limit cycle $\Gamma_{X}$ around $P_{3}$, where $\Gamma_{X}$ collides with $P_{3}$ if $\mbox{tr }DX(P_{3}) = 0$. However, if $\mbox{tr }DX(P_{3}) > 0$, then $P_{3}$ is stable.

Theorem 4. The model (2.2) has a Hopf bifurcation around $P_{3}$ if ${{tr}}DX(P_{3}) = 0$.

Proof. Since

the transversality condition for existence of the Hopf bifurcation is guaranteed ^[43]. □

Moreover, the model (2.2) has a transcritical bifurcation when $s-qE = 0$, as summarized by the following result. In this case, if $s-qE > 0$, then $P_{3} \in \Sigma^{+}$ and $P_{1}$ is a saddle point. For $s-qE = 0$, $P_{3}$ collides with $P_{1}$, so $P_{1}$ is a stable saddle-node equilibrium since $\mbox{det }DX(P_{1}) = 0$ and $\mbox{tr }DX(P_{1}) = -r < 0$. However, if $s-qE < 0$, then $P_{3} \notin \Sigma^{+}$, and $P_{1}$ is locally stable.

Theorem 5. The model (2.2) has a transcritical bifurcation around $P_{1}$ if $s-qE = 0$.

Proof. Let $X_{s}(x, y, s) = \left(0, y\left(1-\frac{y}{n(x-P)}\right)\right)^{T}$ be the derivative of the vector field $X(x, y, s)$ with respect to $s$, $U = \left(\frac{a(K-P)(b+K)}{r(d+K-P)}, -1\right)^{T}$ is the eigenvector associated with the null eigenvalue in the Jacobian matrix $DX$ calculated in $P_{1}$, and $W = (0, -1)^{T}$ is the eigenvector associated with the null eigenvalue in the Jacobian matrix $DX^{T}$ calculated in $P_{1}$.

If $s = s_{0}: = qE$, by the transcritical bifurcation theorem ^[44], three conditions must be guaranteed: $W^{T}X_{s}(P_{1}, s_{0}) = 0$, $W^{T}[DX_{s}(P_{1}, s_{0})U]\neq 0$, and $W^{T}[D^{2}X(P_{1}, s_{0})(U, U)]\neq 0$. Indeed,

● Since $X_{s}(P_{1}, s_{0}) = (0, 0)^{T}$, then $W^{T}X_{s}(P_{1}, s_{0}) = 0$.

● Since $DX_{s}(P_{1}, s_{0})U = (0, -1)^{T}$, then $W^{T}[DX_{s}(P_{1}, s_{0})U] = 1 > 0$.

● If $X = (g_{1}(x, y, s), g_{2}(x, y, s))^{T}$, with $g_{1}$ and $g_{2}$ as described in (2.6), and $U = (u_{1}, u_{2})$, then

and so $W^{T}[D^{2}f(P_{1}, s_{0})(U, U)] = \dfrac{2s_{0}}{n(K-P)} > 0.$ □

To conclude, and as shown in Lemmas 2 and 4, Figure 2 shows the bifurcation diagram of the model (2.2) in the plane $(P, E)$, realized by the numerical bifurcation package SlideCont ^[45], whose selected fixed parameters were chosen in order to show all possible dynamics of the model (2.2) in the bifurcation plane $(P, E)$.

Figure 2.  Bifurcation diagram of the model (2.2) in the plane $(P, E)$ with fixed parameters: $a = 0.3$, $b = 3$, $c = r = 0.5$, $d = 0.001$, $K = 2$, $n = 2.8$, $q = 1$ and $s = 0.1$. Dotted line: $\triangle = 0$. The Normal and dash dotted lines represent the Hodf and Transcritical bifurcation, respectively. Region Ⅰ: $P_{2}, P_{3} \notin \Sigma^{+}$ and $P_{1}$ globally asymptotically stable. Region Ⅱ$_{a}$: $P_{2}, P_{3} \in \Sigma^{+}$ with $P_{3}$ globally asymptotically stable node. Region Ⅱ$_{b}$: $P_{2}, P_{3} \in \Sigma^{+}$ with $P_{3}$ globally asymptotically stable focus. Region Ⅲ: $P_{2}, P_{3} \in \Sigma^{+}$ with $P_{3}$ unstable focus and $\Gamma_{X}$ a globally asymptotically stable limit cycle.

DownLoad: Full-Size Img PowerPoint

In particular, since the trajectories $\varphi_{X}$ and $\varphi_{Y}$ in the model (2.2) do not collide in $\Sigma$, that is, the trajectories $\varphi_{Z_{0}}$ with initial condition $(x(0), y(0)) \in \Sigma^{-}$ must cross $\Sigma$ and remain in $\Sigma^{+}$, the phase portraits characterizing each bifurcation region, shown in Figures 3, are drawn by the two vector fields $X$ and $Y$, separated by the straight line $x = P$, and simulated in Matlab ^[46], where $Y$ can be defined for all $x\leq P$ because the trajectory $\varphi_{Z_{0}}$ with initial condition $p = (P, y(0)) \in \Sigma$ is defined as $\varphi_{Z_{0}}(p, t) = \varphi_{X}(p, t)$ for all $t \in I \cap \{t > 0\}$ and $\varphi_{Z_{0}}(p, t) = \varphi_{Y}(p, t)$ for $t \in I \cap \{t \leq 0\}$. The bifurcation regions are described below.

Figure 3.  Dynamics of the model (2.2) characterizing each bifurcation region shown in Figure 3. The vector fields $X$ and $Y$ are represented by blue and orange colors, respectively, and $\Sigma$ is the black dashed line. Margenta curve: stable limit cycle $\Gamma_{X}$ around $P_{3}$. Red point: $P_{0}$. Pink point: $P_{1}$. Blue point: $P_{2}$. Black point: $P_{3}$.

DownLoad: Full-Size Img PowerPoint

● If $s-qE < 0$, represented in region Ⅰ, the model (2.2) has no interior equilibria, so the $\varphi_{Z_{0}}$ trajectories converge to equilibrium $P_{1}$, as observed in Figure 3(b). However, if $s-qE > 0$, then $P_{3} \in \Omega$, as observed in regions Ⅱ$_{a}$, II$_{b}$ and Ⅲ. Note that the transcritical bifurcation separates region Ⅰ from regions Ⅱ$_{a}$, II$_{b}$, and Ⅲ.

● In regions Ⅱ$_{a}$ and II$_{b}$ we have that $P_{3}$ is a stable node or focus, as observed in Figures 3(c, d), respectively. In both regions, the trajectories $\varphi_{Z_{0}}$ converge to the equilibrium $P_{3}$. Moreover, the curve $\triangle = 0$ separates regions Ⅱ$_{a}$ and II$_{b}$.

● The model (2.2) has a stable limit cycle $\Gamma_{X}$ in $\Sigma^{+}$, characterized in region Ⅲ, and which will not collide with $\Sigma$. Moreover, $P_{3}$ is unstable, so the $\varphi_{Z_{0}}$ trajectories converge to $\Gamma_{X}$, as observed in Figure 3(e). In particular, the Hopf bifurcation separates region Ⅲ from regions Ⅱ$_{a}$ and II$_{b}$.

3.   Prey refuge with alternative food for predators above threshold value

Let $x(t)\geq 0$ and $y(t)\geq 0$ be the population size of prey and predators, respectively, whose dynamics are given by the model (1.5) in the biological sense region:

In this case, the vector field of the model (1.5) is given by

where the vector fields $\tilde{X}$ and $Y$ are defined in $\tilde{\Sigma}^{+} = \left\lbrace (x, y)\in \tilde{\Omega}: P < x\leq K \right\rbrace$ and $\tilde{\Sigma}^{-} = \{(x, y)\in \tilde{\Omega}: 0\leq x < P \}$, respectively.

If $\tilde{\Sigma} = \{(x, y)\in \tilde{\Omega}: x = P \}$ and similarly to the dynamics of the model (2.2) in $\Sigma$, the trajectories $\varphi_{Z_{c}}$ with initial condition $p = (P, y) \in \tilde{\Sigma}$ is defined as $\varphi_{Z_{c}}(p, t) = \varphi_{\tilde{X}}(p, t)$ for $t \in I \cap \{t > 0\}$ and $\varphi_{Z_{c}}(p, t) = \varphi_{Y}(p, t)$ for $t \in I \cap \{t \leq 0\}$ because $\tilde{X}f(p) = Yf(p) = rP\left(1-\frac{P}{K}\right) > 0$.

Analogous to the proof shown in Lemma 1, some results are guaranteed in the model (3.2).

Lemma 5. $\tilde{\Omega}$ is an invariant set and the trajectory $(x(0), y(0)) \in \tilde{\Omega}$ exists and is unique.

3.1.   Existence and local stability at equilibrium points

The equilibria of the model (3.2) on the coordinate axes $(x, y)$ is given by $P_{0}, P_{2} \in \tilde{\Sigma}^{-}$, and $\tilde{P}_{1} = (K, 0) \in \tilde{\Sigma}^{+}$, where $P_{0}$ and $P_{2}$ are equivalent to the equilibria of the model (2.2), so that their local stability is described in Lemma 2. In addition, and equivalently to what is shown in Lemma 2, $\tilde{P}_{1}$ is either a locally stable node or saddle point if $s-qE < 0$ or $s-qE > 0$, respectively.

On the other hand, the intersections in the nullclines of the vector field $\tilde{X}$ in $\tilde{\Sigma}^{+}$

with $y = f_{2}(x)$ as represented in (2.3), and $y = f_{3}(x)$ a line with slope $\frac{nx(s-qE)}{s}$ and intercept at $\left(0, -\frac{(s-qE)(nP-c)}{s}\right)$ and $\left(\frac{nP-c}{n}, 0\right)$, determine the number of interior equilibria of the model (3.2), equivalent to determining the positive roots of the polynomial

with

In this case, it is evident that the nullclines (3.3) do not intersect at $\tilde{\Sigma}^{+}$ if $s-qE < 0$. Similarly, if $s-qE > 0$ and $nP-c > 0$, it follows that there is a unique point of the intersection at the nullclines (3.3) in $\tilde{\Sigma}^{+}$, which corresponds to a unique interior equilibrium $\tilde{P}_{*} = (x_{1}^{*}, y_{1}^{*}) \in \tilde{\Sigma}^{+}$ in the model (3.2), with $y_{1}^{*} = f_{3}(x_{1}^{*})$, as observed in Figure 4, and whose proof is analogous to what is shown in Lemma 3.

Figure 4.  Nullclines (3.3) with $s-qE > 0$. Margenta and turquoise curves are the nullclines $y = f_{2}(x)$ and $y = f_{3}(x)$, respectively. Brown and dark cyan dotted lines are $x = P-d$ and $x = \frac{nP-c}{n}$, respectivaly. Dotted black line is $x = P$.

DownLoad: Full-Size Img PowerPoint

Lemma 6. If $s-qE < 0$, the model (3.2) has no interior equilibria. If $s-qE > 0$ and $nP-c > 0$, the model (3.2) has only one equilibrium.

On the other hand, if $s-qE > 0$ and $nP-c < 0$, then $y_{1}^{*} = f_{3}(x_{1}^{*}) > 0$ for all $x_{1}^{*} > 0$, so the positive roots of the polynomial (3.4) are interior equilibria of the model (3.2). In this case, assuming that $\tilde{P}_{*} = (x_{1}^{*}, y_{1}^{*})$ is an interior equilibrium that always exists in $\tilde{\Omega}$ with $s-qE > 0$ and $nP-c < 0$, a division between the polynomial $\tilde{Q}(x)$ (3.4) and $x_{1}^{*}$ leads to the polynomial

as a factor of $\tilde{Q}(x)$ and the rest of the division is $\tilde{Q}(x_{1}^{*}) = 0$.

By clearing $\tilde{B}+\tilde{A}x_{1}^{*}$ from (3.4) and substituting it in (3.6), we have

If $s-qE > 0$ and $nP-c < 0$, then the following result shows conditions on the coefficients of the polynomial $\bar{Q}(x)$ that determine the number of interior equilibria in the model (3.2).

Lemma 7. Let $s-qE > 0$, $nP-c < 0$ and $\tilde{\triangle} = \left(\tilde{B}+\tilde{A}x_{1}^{*}\right)^{2}+\frac{4\tilde{A}\tilde{D}}{x_{1}^{*}}$.

1) If $\tilde{\triangle} < 0$, the model (3.2) has only one interior equilibrium.

2) If $\tilde{\triangle} = 0$, the model (3.2) has two interior equilibria.

3) For $\tilde{\triangle} > 0$,

(a) If $\tilde{B}+\tilde{A}x_{1}^{*} < 0$, the model (3.2) has three interior equilibria.

(b) If $\tilde{B}+\tilde{A}x_{1}^{*}\geq 0$, the model (3.2) has only one interior equilibrium.

In particular, if the model (3.2) has a unique interior equilibrium point $\tilde{P}_{*}$, and since $\tilde{\Omega}$ is an invariant with equilibria on the unstable coordinate axes, by constraint, we have that $\tilde{P}_{*}$ must be locally either a node or a focus, stable or unstable, so its local stability result is equivalent to what is shown in Lemma 4, where the trace in the Jacobian matrix $D\tilde{X}$ computed on $\tilde{P}_{*}$ is given by:

On the other hand, if the model (3.2) has three different interior equilibria $\tilde{P}_{3} = (x_{1}^{*}, g_{3}(x_{1}^{*}))$, $\tilde{P}_{4} = (x_{-}^{*}, g_{3}(x_{-}^{*}))$ and $\tilde{P}_{5} = (x_{+}^{*}, g_{3}(x_{+}^{*}))$, with $x_{1}^{*} < x_{-}^{*} < x_{+}^{*}$, the eigenvalues of the Jacobian matrix $D\tilde{X}$ computed at each interior equilibrium $\tilde{P}_{i}, i = 3, 4, 5$, shows that $\tilde{P}_{4}$ is a saddle point; therefore, $\tilde{P}_{3}$ and $\tilde{P}_{5}$ must be locally stable or unstable nodes and/or focus, as summarized by the following result.

Lemma 8. If the model (3.2) has three different interior equilibriums $\tilde{P}_{3} = (x_{1}^{*}, f_{3}(x_{1}^{*}))$, $\tilde{P}_{4} = (x_{-}^{*}, f_{3}(x_{-}^{*}))$, and $\tilde{P}_{5} = (x_{+}^{*}, f_{3}(x_{+}^{*}))$, with $x_{1}^{*} < x_{-}^{*} < x_{+}^{*}$, then $\tilde{P}_{4}$ is a saddle point.

For the case of local stability of the equilibria $\tilde{P}_{3}$ and $\tilde{P}_{5}$, the following section generalizes the cases of local stability for both equilibria.

3.2.   Existence of homoclinic curve

Let us consider the case where the model (3.2) has three interior equilibria $\tilde{P}_{3}$, $\tilde{P}_{4}$, and $\tilde{P}_{5}$, with $\tilde{P}_{3}$ stable and $W^{s}_{+}(\tilde{P}_{4})$, $W^{u}_{+}(\tilde{P}_{4})$ the upper stable and unstable varieties to the right in $\tilde{\Sigma}^{+}$ of the saddle point $\tilde{P}_{4}$, respectively, and whose dynamics is observed in Figure 5.

Figure 5.  Relative positions of the varieties $W_{u}^{+}(\tilde{P}_{4})$ and $W_{s}^{+}(\tilde{P}_{4})$, with red and blue curves are $W_{u}^{+}(\tilde{P}_{4})$ and $W_{s}^{+}(\tilde{P}_{4})$, respectively. Margenta and turquoise curves are the nullclines (3.3) as shown in Figure 4.

DownLoad: Full-Size Img PowerPoint

In this case, since $\omega$-limit of $W^{u}_{+}(\tilde{P}_{4})$ could be the equilibrium $\tilde{P}_{5}$, when locally stable, or a stable limit cycle surrounding $\tilde{P}_{5}$ when locally unstable, or the equilibrium $\tilde{P}_{3}$ as observed in Figure 5, if $(\tilde{x}_{0}, \tilde{y}_{a}) \in W^{u}_{+}(\tilde{P}_{4})$ and $(\tilde{x}_{0}, \tilde{y}_{b}) \in W^{s}_{+}(\tilde{P}_{4})$, with $P < \tilde{x}_{0} < K$, $\tilde{y}_{a} = \tilde{f}_{1}(a, b, c, d, E, K, P, q, r, s)$ and $\tilde{y}_{b} = \tilde{f}_{2}(a, b, c, d, E, K, P, q, r, s)$ continuous functions in $X$, then there exist $(\tilde{x}_{0}, \tilde{y}_{a}^{*}), (\tilde{x}_{o}, \tilde{y}_{b}^{*}) \in \tilde{\Sigma}^{+}$ such that $\tilde{y}_{a}^{*} = \tilde{y}_{b}^{*}$; therefore, the varieties $W^{u}_{+}(\tilde{P}_{4})$ and $W^{s}_{+}(\tilde{P}_{4})$ intersect and form a homoclinic curve $\tilde{\gamma}$ created by $\tilde{P}_{4}$ and surrounding $\tilde{P}_{5}$ as shown in the following result.

Theorem 6. If $\tilde{P}_{3}$ is stable and $W^{s}_{+}(\tilde{P}_{4}) \cap W^{u}_{+}(\tilde{P}_{4}) \neq \emptyset$, then the model (3.2) has a homoclinic curve $\tilde{\gamma}$ surrounding $\tilde{P}_{5}$.

In particular, since the trace of $DX$ is calculated in the equilibria $\tilde{P}_{3}$, $\tilde{P}_{5}$ is equivalent to that shown in (3.8), replacing $x_{1}^{*}$ by $x_{-}^{*}$ or $x_{+}^{*}$, respectively, the following result shows a generalization about the local stability of the equilibria $\tilde{P}_{3}$ and $\tilde{P}_{5}$.

Theorem 7. If $\tilde{P}_{3}$ is unstable, that is, ${{tr}}D\tilde{X}(\tilde{P}_{3}) > 0$, then

1) $\tilde{P}_{5}$ is stable if ${{tr}}D\tilde{X}(\tilde{P}_{5}) < 0$.

2) $\tilde{P}_{5}$ is unstable if ${{tr}}D\tilde{X}(\tilde{P}_{5}) > 0$. In this case, at least one stable cycle is formed surrounding $\tilde{P}_{3}$, $\tilde{P}_{4}$, and $\tilde{P}_{5}$.

For $\tilde{P}_{3}$ stable, that is, ${{tr}}D\tilde{X}(\tilde{P}_{3}) < 0$,

1) If $\tilde{y}_{b} < \tilde{y}_{a}$,

(a) $\tilde{P}_{5}$ is unstable if ${{tr}}D\tilde{X}(\tilde{P}_{5}) > 0$

b) $\tilde{P}_{5}$ is stable if ${{tr}}D\tilde{X}(\tilde{P}_{5}) < 0$. In this case, at least one unstable cycle is formed, surrounding $\tilde{P}_{5}$.

2) If $\tilde{y}_{b} > \tilde{y}_{a}$, $\tilde{P}_{5}$ is a stable node or focus.

In addition, if $\bar{\triangle}(\tilde{P}_{i}): = [{{tr}}D\tilde{X}(\tilde{P}_{i})]^{2}-4{det}D\tilde{X}(\tilde{P}_{i})$, $i = 3, 5$, then $\tilde{P}_{i}$ is locally a node or a focus if $\bar{\triangle}(\tilde{P}_{i}) > 0$ or $\bar{\triangle}(\tilde{P}_{i}) < 0$, respectively.

Proof. The Theorem is immediate in view of the results shown in Lemma 4 and the Poincare - Bendixson Theorem ^[40]. Moreover, since there is a stable limit cycle surrounding the equilibrium $\tilde{P}_{3}$, this cycle increases until it coincides with the homoclinic curve $\tilde{\gamma}$ that joins the equilibrium $\tilde{P}_{2}$. Upon breaking $\tilde{\gamma}$, that is, if $\tilde{y}_{b} < \tilde{y}_{a}$, the local stability of $\tilde{P}_{3}$ is not altered. □

3.3.   Bifurcation analysis

In the first instance, if $\tilde{P}$ is a locally non saddle point equilibrium, then the model (3.2) has a Hopf bifurcation if $\mbox{tr }D\tilde{X}(\tilde{P}) < 0$, the proof of which is analogous to what is shown in Theorem 4. Similarly, the model (3.2) has a transcritical bifurcation if $s-qE = 0$, that is, when the equilibrium $\tilde{P}_{*}$ collides with $\tilde{P}_{1}$, and whose proof is similar to that shown in Theorem 4, so it is omitted for brevity. On the other hand, since the model (3.2) has a homoclinic curve $\tilde{\gamma}$ if $\tilde{P}_{3}$ is unstable and the varieties $W^{s}_{+}(\tilde{P}_{4})$ and $W^{u}_{+}(\tilde{P}_{4})$ intersect, then such a model has a homoclinic bifurcation as shown in the following result, the proof of which is deduced from Theorems 6 and 7.

Theorem 8. If the model (3.2) has three interior equilibria $\tilde{P}_{3}$, $\tilde{P}_{4}$ and $\tilde{P}_{5}$, with $\tilde{P}_{3}$ stable and $\tilde{y}_{a}-\tilde{y}_{b} = 0$, the model (3.2) has a homoclinic bifurcation.

Finally, if $\tilde{P}_{3}, \tilde{P}_{4}, \tilde{P}_{5} \in \tilde{\Omega}$, the model (3.2) has a saddle-node bifurcation if the saddle point $\tilde{P}_{4}$ collides with $\tilde{P}_{3}$ or $\tilde{P}_{5}$, forming the saddle-node equilibria $\tilde{P}_{+} = (\tilde{x}_{+}, f_{2}(\tilde{x}_{+}))$ with $\tilde{x}_{+} = -\frac{\tilde{B}+\tilde{A}x_{1}^{*}}{2\tilde{A}}$ or $\tilde{P}_{-} = (\tilde{x}_{-}, f_{3}(\tilde{x}_{-}))$ with $\tilde{x}_{-} = x_{-}^{*}$, as observed in the following result. In particular, given that the Jacobian matrix $D\tilde{X}$ in the vector field $\tilde{X}$ of the model (3.2) calculated in $\tilde{P}_{+}$ is

with $\mbox{det } D\tilde{X}(\tilde{P}_{+}) = 0$ and

then $\tilde{P}_{+}$ is a stable or unstable saddle-node equilibrium if $\mbox{tr } DX(\tilde{P}_{+}) < 0$ or $\mbox{tr } DX(\tilde{P}_{+}) > 0$, respectively. Analogously for $\tilde{P}_{-}$ if $\tilde{x}_{+}$ is replaced by $\tilde{x}_{-}$.

Theorem 9. If $\tilde{P}_{3}$ and $\tilde{P}_{4}$, or $\tilde{P}_{5}$ and $\tilde{P}_{4}$, collide, the model (3.2) has a saddle-node bifurcation around $\tilde{P}_{+}$ o $\tilde{P}_{-}$, respectively.

Proof. The case in which $\tilde{P}_{4}$ collides with $\tilde{P}_{5}$ will be proved, similarly for the case in which $\tilde{P}_{3}$ and $\tilde{P}_{4}$ collide. Let $\tilde{X}_{c}(x, y, c)$ be the derivative of the vector field $\tilde{X}(x, y, c)$ with respect to $c$, that is, $\tilde{X}_{c}(x, y, c) = \left(0, \frac{sy^{2}}{[n(x-P)+c]^{2}}\right)^{T}$, and $U = \left(\frac{s}{n(s-qE)}, 1\right)^{T}$, $W = \left(\frac{(s-qE)[(\tilde{x}^{+}-P)+d]}{a(\tilde{x}^{+}-P)}, -1\right)$ the eigenvectors associated with the null eigenvalue in $D\tilde{X}$ and $D\tilde{X}^{T}$ calculated in $\tilde{P}_{+}$, respectively.

If $c_{0}$ satisface $\tilde{\triangle}(c_{0}) = 0$, by the saddle-node bifurcation Theorem ^[44], two conditions must be guaranteed: $W^{T}\tilde{X}_{c}(\tilde{P}_{+}, c_{0})\neq 0$ and $W^{T}[D^{2}\tilde{X}(\tilde{P}_{+}, c_{0})(U, U)]\neq 0$. Indeed,

1) Since $\tilde{X}_{c}(\tilde{P}_{+}, c_{0}) = \left(0, \frac{(s-qE)^{2}}{s}\right) ^{T}$ then $W^{T}\tilde{X}_{c}(\tilde{P}_{+}, c_{0})\frac{(s-qE)^{2}}{s} > 0$.

2) If $\tilde{X} = (\tilde{g}_{1}(x, y, c), \tilde{g}_{2}(x, y, c))^{T}$, with $\tilde{g}_{1} = \frac{rx}{x+b}\left(1-\frac{x}{K}\right)-\frac{a(x-P)y}{(x-P)+d}$ and $\tilde{g}_{2} = sy\left[1-\frac{y}{n(x-P)}\right]-qEy$, and $U = (u_{1}, u_{2})$, then

and $W^{T}[D^{2}\tilde{X}(\tilde{P}_{+}, c_{0})(U, U)] = \dfrac{2s(3\tilde{A}\tilde{x}_{+}+\tilde{B})}{n^2K[(x^{+}-P)+b](s-qE)^{2}} \neq 0.$ □

Figure 6 shows the bifurcation curves in the $(P, E)$ plane, and realized by the numerical bifurcation package SlideCont ^[45], which the transcritical bifurcation separates regions Ⅰ and Ⅱ, the homoclinic bifurcation separates regions Ⅵ and Ⅶ and the Hopf bifurcation divides regions Ⅱ and Ⅲ, Ⅴ and Ⅵ, and Ⅶ and Ⅷ. In this case, the stability of $\tilde{P}_{*}$, $\tilde{P}_{3}$ and $\tilde{P}_{5}$ are altered in the curves separating regions Ⅱ-Ⅲ, Ⅴ-Ⅵ, and Ⅶ-Ⅷ, respectively.

Figure 6.  Bifurcation diagram of the model (3.2) in the plane $(P, E)$ with fixed parameters: $a = 0.3$, $b = 3$, $c = 5$, $d = 0.001$, $K = 2$, $n = 2.8$, $q = 1$, $r = 0.5$ and $s = 0.1$. Dash dotted line: Transcritical bifurcation. Normal line: Hopf Bifurcation. Dash dot dotted line: Saddle-node bifurcation. Dashed line: Homoclinic bifurcation. Region Ⅰ: $P_{2}$, $\tilde{P}_{*}$, $\tilde{P}_{i} \notin \tilde{\Omega}$, $i = 3, 4, 5$, and $\tilde{P}_{1}$ globally asymptotically stable. Region Ⅱ: $P_{2} \in \tilde{\Omega}$ and $\tilde{P}_{*} \in \tilde{\Omega}$ globally asymptotically stable. Region Ⅲ: $P_{2} \in \tilde{\Omega}$, $\tilde{P}_{*} \in \tilde{\Omega}$ unstable and existence of a globally asymptotically stable limit cycle $\Gamma_{\tilde{X}}$ around $\tilde{P}_{*}$. Region Ⅳ: $P_{2} \in \tilde{\Omega}$, $\tilde{P}_{i} \in \tilde{\Omega}$ unstables, $i = 3, 4, 5$, and existence of a globally asymptotically stable limit cycle. $\Gamma_{\tilde{X}}$ around $\tilde{P}_{i}$. Region Ⅴ: $P_{2} \in \tilde{\Omega}$, $\tilde{P}_{i} \in \tilde{\Omega}$, $i = 3, 4, 5$, with $\tilde{P}_{3}$ unstable and $\tilde{P}_{5}$ globally asymptotically stable. Region Ⅵ: $P_{2}$, $\tilde{P}_{i} \in \tilde{\Omega}$, $i = 3, 4, 5$, with $\tilde{P}_{3}$ and $\tilde{P}_{5}$ locally asymptotically stable. Region Ⅶ: $P_{2}$, $\tilde{P}_{i} \in \tilde{\Omega}$, $i = 3, 4, 5$, with $\tilde{P}_{3}$, $\tilde{P}_{5}$ locally asymptotically stable and an unstable limit cycle $\Gamma_{\tilde{X}}$ around $\tilde{P}_{5}$. Region Ⅷ: $P_{2}$, $\tilde{P}_{i} \in \tilde{\Omega}$, $i = 3, 4, 5$, with $\tilde{P}_{3}$ globally asymptotically stable and $\tilde{P}_{5}$ unstable.

DownLoad: Full-Size Img PowerPoint

Furthermore, in regions Ⅴ and Ⅷ, we observe the existence of three interior equilibria in the model (3.2), which the chair-node bifurcation formed between regions Ⅱ and Ⅴ shows the collision of $P_{3}$ and $\tilde{P}_{4}$. However, the collision between $\tilde{P}_{4}$ and $\tilde{P}_{5}$ is shown by the curve between regions Ⅶ and Ⅷ. Similarly, the curve separating regions Ⅴ and Ⅲ represents a saddle-node bifurcation.

On the other hand, since the selected fixed parameters were chosen in order to show all possible dynamics of the model (3.2) in the bifurcation plane $(P, E)$ in Figure 6, Figure 7 shows the phase portraits of each bifurcation region, whose dynamics are described as:

Figure 7.  Dynamics of the model (3.2) describing bifurcation regions Ⅰ to Ⅵ shown in 6. The vector fields $\tilde{X}$ and $Y$ are represented by green and orange colors, respectively, and $\tilde{\Sigma}$ is the black dashed line. Blue line: $\varphi_{Z_{c}}$. Blue curve: unstable limit cycle $\Gamma_{\tilde{X}}$. Margenta curve: stable limit cycle $\Gamma_{\tilde{X}}$. Red point: $P_{0}$. Blue point: $P_{2}$. Pink point: $\tilde{P}_{1}$. Green point: $\tilde{P}_{3}$. Orange point: $\tilde{P}_{4}$. Gray point: $\tilde{P}_{5}$. Black point: $\tilde{P}_{*}$.

DownLoad: Full-Size Img PowerPoint

● Regions Ⅰ-Ⅲ are equivalent to that shown in Figure 2, so there is a single equilibrium $\tilde{P}^{*}$ and its local stability is observed in Figures 7(a, b, c).

● For the case where $\tilde{P}_{3}$ and $\tilde{P}_{5}$ are unstable equilibria, we observe the formation of a limit cycle $\Gamma_{\tilde{X}}$ surrounding the three interior equilibria in region Ⅳ as observed in Figure 7(d), so that the trajectories $\varphi_{Z_{c}}$ converge to $\Gamma_{\tilde{X}}$.

● If $\tilde{P}_{3}$ is unstable and $\tilde{P}_{5}$ is stable, in region Ⅴ the trajectories $\varphi_{Z_{c}}$ converge to $\tilde{P}_{3}$ as observed in Figure 7(e).

● For the case where $\tilde{P}_{3}$ and $\tilde{P}_{5}$ are stable, we have either the formation of a locally unstable limit cycle surrounding $\tilde{P}_{5}$, or a homoclinic curve joining $\tilde{P}_{4}$ and surrounding $\tilde{P}_{5}$, or the breaking of such a curve without altering the stability of $\tilde{P}_{5}$. In this case, in region Ⅶ we have a locally unstable limit cycle as observed in Figure 7(g), which vanishes as shown in region Ⅵ of Figure 7(f).

● In region Ⅷ we have that $\varphi_{Z_{c}}$ converges to $\tilde{P}_{3}$ as observed in Figure 7(h).

4.   Conclusions

Since the predator-prey Leslie-Gower models proposed by several researchers which consider constant prey refuge are valid only when the initial population size of the prey is above their refuge ^[1,29], that is, when $P < x(0)$, in this paper, we make a general modification to these models in order to show the dynamics of both species without restrictions on the initial condition of the proposed Leslie-Gower models.

On the other hand, given the hypotheses on the dynamics of predators in the presence or absence of prey for their food, the proposed model (1.3) considers that the predator's diet depends only on prey; in the absence of prey, they are forced to obtain other types of food for their subsistence. However, the model (1.5) considers that the predator's diet is given by prey and alternative food provided by the environment, and in the absence of prey, predators are not forced to change their diet.

In both proposed models (1.3) and (1.5), all trajectories with initial condition in the vector field $Y$ cross the line $x = P$, and remain in the vector field $X$ or $\tilde{X}$, respectively, because the system of differential equations describing the vector field $Y$ are independent and the only stable interior equilibrium does not belong to the biological sense region, though the trajectories try to make a convergence to this point. Similarly, the stable or unstable varieties of each interior equilibrium in both models do not cross the straight line $x = P$, unlike the unstable variety of $P_{0}$, which connects to $P_{1}$ or $\tilde{P}_{1}$, and $P_{2}$. Biologically, the growth of the prey described by both models (1.3) and (1.5) cannot be controlled below $P > 0$, so prey growth over time is higher than $P$.

The model (1.3) has a unique interior equilibrium point if the intrinsic birth rates of the predators are greater than the product of the catchability coefficient and hasvesting effort of the predators, that is $s-qE > 0$; otherwise, the prey converge to their carrying capacity $K > 0$ and the predators become extinct. Moreover, the model (1.3) has a stable limit cycle as long as the interior equilibrium is unstable, so both species stabilize over time if the interior equilibrium is stable. In this model we have two types of bifurcations: Transcritical and Hopf. These are what separates the extinction or not of predators and the long-term stabilization of the population size of both species, respectively.

The bifurcation cases of model (1.3) can be transferred to model (1.5). However, model (1.5) has two additional bifurcations: saddle-node and homoclinic. In this case, the model can have between one or three interior equilibria as long as $s-qE > 0$; otherwise, analogous to model (1.3), predators die out and prey converge to the carrying capacity.

For the case where model (1.5) has three interior equilibria, there exists a stable limit cycle surrounding all its interior equilibria if $\tilde{P}_{3}$ and $\tilde{P}_{5}$ are unstable. Similarly, the model (1.5) could have a locally unstable limit cycle if $\tilde{P}_{3}$ and $\tilde{P}_{5}$ are stable, and the homoclinic bifurcation is given by the formation of a homoclinic curve $\tilde{\gamma}$ given by $\tilde{P}_{4}$.

In addition to the homoclinic bifurcation, the model (1.5) could have a heteroclinic bifurcation, whose bifurcation curve coincides with the homoclinic, from which a variety connection transition is made between $\tilde{P}_{1}$ and $\tilde{P}_{3}$ or $\tilde{P}_{1}$ and $\tilde{P}_{5}$, where the locally unstable limit cycle is given by the formation of the heteroclinic curve connecting $\tilde{P}_{1}$ and $\tilde{P}_{5}$.

On the other hand, if the model (1.5) has a unique interior equilibrium, equivalently to the model (1.3), we have that the growth of both species stabilizes with time, or on the contrary has oscillatory solutions, when the equilibrium is stable or unstable, respectively. However, if the model (1.5) has three interior equilibria, the growth of the prey converges to its minimum quantity, regardless of the initial condition, if $\tilde{P}_{3}$ is stable and $\tilde{P}_{5}$ is unstable. However, if $\tilde{P}_{3}$ and $\tilde{P}_{5}$ are stable, and there is no unstable cycle surrounding $\tilde{P}_{5}$, then prey converges to its minimum quantity if the initial condition of the species is above the stable variety of $\tilde{P}_{4}$, otherwise prey converges to its maximum quantity. Similarly, if the model (1.5) has three equilibria where $\tilde{P}_{3}$ and $\tilde{P}_{5}$ are stable, with an unstable cycle surrounding $\tilde{P}_{5}$, prey converges to its minimum quantity if the initial condition is above the limit cycle.

Similarly, if the model (1.5) has three interior equilibria, where $\tilde{P}_{3}$ is unstable and $\tilde{P}_{5}$ is stable, the prey converge to their maximum quantity regardless of their initial condition. However, if $\tilde{P}_{3}$ and $\tilde{P}_{5}$ are unstable, both species have oscillatory solutions. In particular, alternative food for predators could lead to a possible maximum or minimum choice in the amount of prey over time as long as the model (1.5) has three interior equilibria, with appropriate stability for each of these, and an appropriate selection in the initial condition of prey and predators.

To conclude, if in both proposed models (1.3) and (1.5) no external factors inhibiting the intrinsic growth of prey are considered, that is, if $\dot{x} = rx\left(1-\frac{x}{K}\right)$ in the absence of predators, the bifurcation cases for the two models (1.3) and (1.5) hold, and the interior equilibria still need to be computed computationally, so the alteration of $\dot{x} = rx\left(1-\frac{x}{K}\right)$ or $\dot{x} = \frac{rx}{x+b}\left(1-\frac{x}{K}\right)$ in the absence of predators does not affect the overall mathematical results in both proposed models.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

This research was supported by MCIN/AEI/10.13039/501100011033 through grant BADS, no. PID2019-109320GB-100, and the associated FPI contract PRE2019-088899 to Christian Cortés García. The Spanish MICINN has also funded the "Severo Ochoa" Centers of Excellence to CNB, SEV 2017-0712.

Conflict of interest

The author declare there is no conflict of interest.