Dynamics of a Two Prey and One Predator System with Indirect Effect

: We study a population model with two preys and one predator, considering a Holling type II functional response for the interaction between ﬁrst prey and predator and taking into account indirect effect of predation. We perform the stability analysis of equilibria and study the possibility of Hopf bifurcation. We also include a detailed discussion on the problem of persistence. Several numerical simulations are provided in order to illustrate the theoretical results of the paper.


Introduction
Population dynamics has been extensively studied by researchers in biomathematics, particularly the predator-prey models.
In [1], the authors consider the following two prey one predator model: where x, y, z represent the population densities of the two preys and of the predator, respectively. In the previous model, for the interaction between the first prey and the predator, they considered a Holling type II functional response where the handling time of predator for the second prey is also involved, whereas for the interaction between the second prey and the predator, they considered a Lotka-Volterra functional response. It is also assumed that there is no intraspecific interaction in the second prey population and its growth is exponential; as a consequence, there is a huge availability of second prey in the absence of a predator, and there is no searching time for the second prey population. They found necessary and sufficient conditions for existence and stability of the nontrivial equilibrium E * (see [1]).
In order to recover a more complex behavior, we consider a modification of the model that takes the indirect effects of predations into account.
The role played by indirect effects in population dynamics has been investigated in the last several decades (see [2][3][4][5][6][7][8][9][10][11][12]). In the case of predation, it has been pointed out (see [13]) that predator can alter the morphology (see [6]) or the behavior of the preys. In particular, the preys, in order to avoid contacts with predators, may reduce their normal activity or may stay hidden most of the time. Many kinds of indirect effects have been described in the literature (see, for example, [8] for a detailed discussion); an interesting example (see [14]) is the case of refuge indirect effect.
A model was proposed in order to take into account indirect interactions in a plankton community (see [15] and the references quoted therein). They analyzed the effects of predator Daphnia over two groups of phytoplankton of different morphology (see [7,9]), having phosphorous as a resource (see [5,9]). In this case, the predator prefers to predate the smaller size prey group and the other one take advantages of it.
The model was analytically studied in [15] using persistence theory (see [16,17]) and in [18] using bifurcation theory. Both studies suggest the importance of indirect effects of predation in order to describe cases of coexistence in real life. In [19], seasonal indirect effects were considered to show the possibility of chaotic motion, whereas in [20,21], the authors considered the stochastic version of the model. Regarding the model (1), since there is a higher availability of the second prey, it is natural to suppose that the predator prefer to predate the second prey and the first one take advantages of it. The easiest way to model this situation consists of adding the indirect effect term −Lyz in the second equation and the term Lyz in the first one, with the parameter L > 0 describing the intensity of indirect effects. The systems becomė We will consider initial conditions x(0) ≥ 0, y(0) ≥ 0, z(0) ≥ 0 and we assume all the parameters are positive and with the following meaning: r and k are the intrinsic growth rate and carrying capacity of the first prey, respectively; β and δ are the intrinsic growth rate and predation rate of the second prey, respectively; a is the half saturation value of the predator; b is the maximum growth rate of the predator; c is the maximum rate of predation for first prey item; m is the death rate of the predator in the absence of prey; α is the quotient of the handling time of the predator per second prey item and the handling time of the predator per first prey item; η is the quotient of the capture rate of the second prey and the capture rate of the first prey; γ is the efficiency with which the second prey consumed by the predator gets converted into predator biomass (see [1] for more details).
We note that the system is not of Kolmogorov type; indeed, the first equation cannot be written in the formẋ = x f (x, y, z).
Such systems can be regarded as semi-Kolmogorov systems using a terminology introduced in [19].
For system (2), we perform the stability analysis of equilibria and we analyze the existence of limit cycles by Hopf bifurcation, as they play an important role in the qualitative theory of differential systems. The study of limit cycles was initiated by Poincaré [22] and motivated by the famous 16th Hilbert problem [23][24][25] and by the fact that the behavior of many natural phenomena has been modelized by limit cycles, as, for instance, the famous limit cycle of van der Pol [26].
In the last several decades, the existence of limit cycles for systems with biological meaning, such as the Lotka-Volterra or Kolmogorov system, have been studied through Hopf and zero-Hopf bifurcation (see, for example, [27,28] for recent results).
The rest of the paper is organised as follows: in Section 2, we present a preliminary analysis of the features of the model. In Section 3 we provide a study of existence and stability of equilibria. In Section 4, we discuss the problem of persistence of the three species. In Section 5, we present a case of Hopf bifurcation. Finally, Section 6 contains some conclusive remarks.

Analysis of the System (2) on the Invariant Planes
First, we will show that the dynamics of the system, considering positive initial conditions, is contained in the first octant. Theorem 1. The set {(x, y, z) ∈ R 3 : x, y, z ≥ 0} is positively invariant for system (2).
Proof. At first, we must note that the planes z = 0 and y = 0 are invariant. On the plane x = 0, we haveẋ = Lyz > 0; then, solutions do not leave the positive octant, that is, the set {(x, y, z) ∈ R : x, y, z ≥ 0} is positively invariant.
We analyze the dynamics on the boundary of {(x, y, z) ∈ R : x, y, z ≥ 0}. In order to do that, we first study the dynamics on the coordinate axes and then on the planes y = 0 and z = 0.
The three axes are invariant for the dynamics; in particular, any solution with initial conditions on the x-axis tends to the equilibrium (k, 0, 0), any solution with initial conditions on the z-axis tends to the equilibrium (0, 0, 0), and any solution with initial conditions on the y-axis verifies that y(t) tends to infinity when t tends to infinity. See Figure 1. Now, we present some considerations about the dynamics on the invariant planes. On the invariant plane, z = 0 the system iṡ and for this system, solutions are unbounded except that on the positive x-axis. The equilibrium points are (0, 0) and (k, 0). The eigenvalues of DF(0, 0) are r and β, which are both positive, so this equilibrium is an unstable node. The eigenvalues of DF(k, 0) are −r and β, so the equilibrium point is a saddle. On the plane y = 0, the system iṡ and the equilibria are (0, 0), (k, 0). If there exists a further equilibrium point with coordinates (m/(b − m),z) wherē The eigenvalues of DF(0, 0) are r and −m, so it is a saddle. The eigenvalues of DF(k, 0) are −r and bk/(a + k) − m. The first one is negative and the second changes it sign when bk = (a + k)m. When bk/(a + k) − m < 0, it is a stable node, and when bk/(a + k) − m > 0, it loses its stability; it becomes a saddle and the third equilibrium appears.
The x-nullclines are x = 0 and z = (r/(ck))(k − x)(a + x). If bk/(a + k) = m, theṅ z = (ba(x − k))/((a + x)(a + k))z is positive if x > k and negative if x < k. The local phase portrait in this case is given in Figure 2.
x-nullclines z-nullclines When bk/(a + k) > m, the equilibrium (ma/(b − m),z) appears, as is it shown in Figure 3. The eigenvalues of the Jacobian matrix of system (4) at E 2 are where x z E0 E1 1 1 1 E2 2 x-nullclines z-nullclines Note that A 3 > 0 by the existence conditions of the equilibrium point. These eigenvalues are complex if A 2 < 0, and in that case, they have a positive real part if in which case, E 2 is unstable. If the real part of the eigenvalues is negative, so E 2 is asymptotically stable. In the case with A 2 > 0, as the determinant of the Jacobian matrix is positive, it is not possible that the eigenvalues have different sign. Then, if A 1 is positive, both eigenvalues are positive and E 2 is unstable. In other case, from the conditions which is a contradiction. In the same way, we obtain that if A 1 < 0, then both eigenvalues are negative and E 2 is asymptotically stable.
We give here some results about the possible existence of periodic orbits surrounding the equilibrium point E 2 in the plane y = 0.
The equilibrium point E 2 is a Hopf equilibrium if and only if A 1 = 0 and A 2 < 0, it is, when m = b(k − a)/(k + a). Note that this occurs only for a < k. In general, when a differential systemẋ = F(x, µ) in R n has an equilibrium x 0 with eigenvalues ±ωi, it can exhibit a Hopf bifurcation, that is, a local bifurcation in which the equilibrium point loses stability as a pair of complex conjugate eigenvalues of the linearization around the equilibrium point, cross the imaginary axis of the complex plane. To show that this bifurcation takes place, it is necessary to compute the first Lyapunov coefficient 1 (x 0 ) of the differential system at the equilibrium. When 1 (x 0 ) < 0, the equilibrium x 0 is a weak focus of the differential system restricted to the central surface of x 0 , associated to the pair of complex eigenvalues, which cross the imaginary axis, and the limit cycle that emerges from x 0 is stable. In this case, we say that the Hopf bifurcation is supercritical.

Theorem 2.
The equilibrium E 2 of system (4) undergoes a supercritical Hopf bifurcation at m 0 = b(k − a)/(a + k) > 0. For m < m 0 , the system has a unique and stable limit cycle bifurcating from the equilibrium point E 2 .
Proof. We use the results presented on Chapter 3 of [29] for computing the first Lyapunov coefficient 1 at the equilibrium E 2 . At first, to simplify calculation, we introduce in system (4) a new time variable τ by dt = (a + x)dτ, obtaining the polynomial system: This system has the positive equilibrium which is the same as (am/(b − m), z) with the notation introduced in Section 2. The Jacobian matrix at this equilibrium is We get µ(m 0 ) = 0 for which is positive because as we have said before, a neccesary condition for Hopf bifurcation is a < k. Moreover, Therefore, at m = m 0 , the equilibrium point E 2 has a pair of pure imaginary eigenvalues ±iω(m) and the system has a Hopf bifurcation. The equilibrium is stable for m > m 0 and unstable for m < m 0 . In order to analyze this Hopf bifurcation, we will apply Theorem 3.3 in [29], so we must prove if the genericity conditions are satisfied. We check that the transversality condition is satisfied as where represents the derivative with respect to m, and the sign is determined because a < k.
To check the second condition, we must compute the first Lyapunov coefficient. We fix the value m = m 0 , and then, the equilibrium E 2 has the expression We translate E 2 to the origin of coordinates obtainig the systeṁ which can be represented asε where A = A(m 0 ) and the multilinear functions B and C are given by We need to find two eigenvectors p, q of the matrix A verifying Aq = iωq, A T p = −iωp, and < p, q >= 1, as for example Now, we compute , g 21 = p, C(q, q, q) = − 3r 4kω 2 , and the first Lyapynov coefficient which is negative for any values of the parameters, and so the second condition of the theorem we are applying is satisfied, and we can conclude that a unique and stable limit cycle bifurcates from the equilibrium point E 2 through a Hopf Bifurcation for m < m 0 . Now, we include some numerical experiments. We fix the parameters as follows: In this case, m 0 = 0.1667 and the eigenvalues of J(E 2 ) are In Figure 4, we represent the case m = 0.18 > m 0 , in which the equilibrium is locally asymptotically stable. In Figure 5, we represent the case m = 0.14 < m 0 , in which the equilibrium loses stability and a limit cycle arises due to Hopf bifurcation.  z(t) Figure 5. The time histories and the solution for m = 0.14 < m 0 . In this case, the equilibrium is no more locally asymptotically stable and a stable limit cycle attracts nearby solutions.
We conclude this section by considering a case in which there are no periodic orbits in {y = 0} ∩ R 3 + : then system (4) does not admit periodic orbits in the set {(x, z) ∈ R 2 : x, z ≥ 0} .
Proof. Let In order to prove the nonexistence of periodic orbits, we use the Bendixson-Dulac theorem that states that if there exists a function ϕ(x, z) such that the term does not change sign in a simply connected set S, then there are no periodic orbits on S. We consider then function ϕ(y, z) = (a + x)/x; then: We observe that, sinceẋ < 0 for x ≥ k, there are no periodic orbits in the set and for the same reason, there are no periodic orbits crossing the half line {x = k, z ≥ 0}. As a consequence we will restrict to the case x ≤ k for which we obtain and we conclude that there are no periodic orbits in the whole set {(x, z) ∈ R 2 + }. Remark 1. We observe that, for a < k, the condition of Theorem 3 on the nonexistence of periodic orbits and then is compatible with the results of bifurcation analysis.

Existence and Stability Analysis of Equilibria
The first step for studying the dynamics of the system (2) is to find all the equilibrium points and analyzing their stability. • If kb > m(a + k), the equilibrium Proof. From direct calculation.
We also analyze the existence of nontrivial positive equilibria for system (2).
Theorem 5. The system (2) has at least one positive equilibrium E * (x * , y * , z * ) if and only if where x * is a solution of the equation with the coefficients C i , i = 0, ..., 4 defined below.
Proof. By direct calculation, we obtain the equilibrium E * = (x * , y * , z * ), with and x * is the solution of the equation We must require that x * verifies Otherwise, this is satisfied when The expressions of the coefficients in Equation (21) are given by We apply Descartes rule of signs in Equation (21) to study the existence of positive roots. The coefficient of degree zero, C 0 , is always negative and the coefficient of degree four, C 4 , is always positive. This means that there always exist a real positive and a real negative zero of the polynomial. The other two zeroes can be complex or real with the same sign.
The other coefficients of the polynomial can be positive, negative or zero, and combining all the possible signs we obtain that: • If one of the following conditions holds The signs of C 3 , C 1 and −C 2 are equal. then there exists three or one positive roots of Equation (21). • In any other case, there exists one positive root of Equation (21).
So that there exist at least one non-trivial interior equilibrium E * = (x * , y * , z * ) of the system (2) if condition (19) is satisfied.

Remark 2.
The positive equilibria, if they exist, are all on the plane z = β/(δ + L). • There exist systems for which there is only one positive solution x * of Equation (21), and this solution is such that verifies condition (19). As an example the system with parameters: • There are systems for which there is only one positive solution x * of Equation (21), and this solution does not verify condition (19), so there are no positive equilibria for system (2). As an example the system with parameters: • There are systems for which there are three positive solutions x * of Equation (21), and only for one of this solutions condition (19) holds, so there are one positive equilibrium for system (2).
As an example the system with parameters: • Also there are systems for which there are three positive solutions x * of Equation (21), and none of them verifies condition (19), so system (2) has none positive equilibria. As an example the system with parameters: We have not found, numerically, any conditions for which Equation (21) has three positive solutions x * and the three of them verify condition (19) but we are not able to exclude this case analytically.
We will analyze the local stability of the equilibria. To do this we consider the Jacobian functional of the vector field and where we have set for simplicity ω = αη. Theorem 6. The stability of the boundary equilibria is the following: • The equilibrium point E 0 is always a saddle, the z-axis is the stable manifold and it has an unstable manifold of dimension two.

•
The equilibrium point E 1 is a saddle with a stable manifold of dimension one if bk > m(a + k) and with a stable manifold of dimension two if bk < m(a + k).

•
The equilibrium point E 2 is unstable if β > z(L + δ) or if β < z(L + δ) and one of the following statements holds: -A 2 < 0 and condition (7) holds, -A 2 > 0 and A 1 > 0. The equilibrium point E 2 is asymptotically stable if β < z(L + δ) and one of the following statements holds: -A 2 < 0 and condition (8) holds, -A 2 > 0 and A 1 < 0. where the coefficients A i are the ones given in (6).
Proof. The local stability analysis of the equilibria E 0 and E 1 is the same as in the case without indirect effects.
For the equilibrium E 0 = (0, 0, 0) the Jacobian matrix is so there are two positive eigenvalues and one negative eigenvalue so that E 0 is a saddle with a stable manifold of dimension one, which is the z-axis, and an unstable manifold of dimension two.
For the equilibrium E 1 = (k, 0, 0) we have When the value of bk surpasses the value of m(a + k), a new equilibrium E 2 = (ma/(b − m), 0,z) appears, and the equilibrium E 1 changes from having an unstable manifold of dimension two to an unstable manifold of dimension one. The functional Jacobian of E 2 is with eigenvalues with the coefficients A i given in (6).
The sign of eigenvalues λ 2,3 has been analyzed in Section 2. We note that the expression can be written as z(L + δ), so we conclude that the eigenvalue λ 1 is positive if and it is negative in the other case. Combining the different possibilities for the three eigenvalues we obtain the conditions for the stability of E 2 .

Remark 4.
We recall that in the case β = z(L + δ), the equilibrium point E 2 has two non-zero eigenvalues in the plane y = 0, as stated in Section 2, and the third eigenvalue is zero. The direction on the flow in the y-direction is determined by the z coordinate. Note that in this casė y = y(L + δ)(z − z), soẏ is positive if z < z andẏ is negative if z > z.

Remark 5.
We have seen that if (22) is verified then there exists at least one positive equilibrium E * . We can prove that, at least in this case, the first eigenvalue λ 1 of J(E 2 ) is positive, that is E 2 is unstable. We recall that λ 1 is positive if In this case we observe thatz ≤ kbr 4cm , in fact can be rewritten as which is always verified.
We use Hurwitz criterion to study the stability of the equilibrium point E * , and we conclude that it is asymptotically stable if and only if s 1 , s 2 , s 3 > 0 and also the expression s 1 s 2 − s 3 is positive.

Remark 6.
We observe that if x * ∈ k 2 + k 2 4 − kβcm br(δ+L) , ∞ then f 1 x is negative, that is Trace(J(E * )) < 0. This means that at least one eigenvalue has negative real part.

Some Remarks about Coexistence of the Three Species
The problem of the coexistence of the three species can be reformulated in mathematical terms by finding the conditions for which the positive solutions starting in the interior of R 3 + do not approach the boundary of the set, ∂R 3 + as t → +∞. These ideas can be made rigorous in the context of persistence theory (see [16]). There are several definitions that are used by mathematicians depending on the context. If we consider a nonlinear system of the formẋ then we have persistence (see [30]) if while we have permanence (see [31]) if there exists m, Finally, we have uniform persistence (see [32]) if there exists ε > 0, independent of In many cases, the favorite choice for analysis is uniform persistence, since in real cases, requiring that lim sup t→+∞ x i (t) > 0 is not sufficient. In fact, a small stochastic or non-autonomous perturbations may lead solutions converges to the boundary. For this reason, in general, it is important to require that lim sup t→+∞ x i (t) ≥ ε > 0. We recall the definition introduced in Unfortunately, in the case of system (29), we cannot prove the existence of a compact attractor of bounded set of R 3 + and as a consequence we are not able to obtain uniform persistence of the system. This lack of dissipativity can be guessed by considering the unbounded solutions on the invariant plane z = 0. Moreover, we recall that the divergence of the vector field F associated to the system is related to the evolution of the three dimensional volumes elements under the flow of the system. We observe that it has a complex expression, in particular, there are no values of the parameters for which it has a negative sign (which means contraction), on a region of ∂R 3 + . Moreover, while the first term and the second term are negative for high values of x and z respectively, the third term is positive for high value of y.
Then, we are able only to prove conditions (H2)-(H3) which only ensures that the invariant sets of ∂R 3 + does not attracts positive solutions. For condition (H2) we first recall the definition: + is called is called weakly ρ−repelling if there is no solution x(t) of system (29) starting at x 0 , with ρ(x 0 ) > 0, such that x(t) → M as t → +∞. Theorem 8. We suppose that hypothesis of Theorem 3 is verified. Then the invariant set of ∂R 3 + are weakly ρ−repelling in any of the following cases 1.
bk < k(a + m); 2. bk > k(a + m) and β > ab(bk − m(a + k))r(L + δ) Proof. In order to prove the theorem, it is sufficient to show that the stable manifolds of the invariant sets of ∂R 3 + are contained all contained in ∂R 3 + . The equilibria E 0 and E 1 have their stable manifolds on ∂R 3 + for any value of the parameters. Then if b − k(m + a) < 0, there are no further equilibria. Otherwise, if E 2 exists, a sufficient condition that ensures that its stable manifold is in ∂R 3 + consists in requiring that that is, its first eigenvalues is positive. To conclude the proof we have to exclude the existence of other invariant set contained in ∂R 3 + . We use Theorem 3 that guarantees the non existence of periodic orbits in the plane y = 0, that is the only invariant sets of ∂R 3 + are the equilibria. Now we pass to check condition (H3).

Definition 3. Let A, B ⊂ ∂R 3
+ . Then A is chained to B in ∂R 3 + , and we write A ⇒ if there exists a total trajectory x(t) with x(0) / ∈ A ∪ B such that Thanks to this notion we are able to exclude the case in which orbits are attracted by an heteroclinic cycle that connects the equilibria such as in the well known case of [33]. Proof. By hypothesis the only invariant sets of ∂R 3 + are equilibria. The analysis of the previous sections is sufficient to exclude the existence of cycle between equilibria. In Figure 6 below we represent the case in which there are only two boundary equilibria (E 0 and E 1 ) while in Figure 7 we represent the case in which we have also E 2 . In the latter situation we distinguish two cases: both eigenvalues with positive and both with negative real part respectively. In conclusion, we are not able to prove uniform persistence, however the above results guarantee a sort of weak persistence of the three species.

Hopf Bifurcation
In this section, we analyze the possible existence of a limit cycle by Hopf bifurcation for the positive equilibrium E * . We recall that Hopf bifurcation occurs when a pair of complex conjugate eigenvalues of the Jacobian matrix of an equilibrium crosses the imaginary axis. In this case a limit cycle arises and its stability character can be obtained by the analysis of the first Lyapunov coefficient. If it is negative the cycle is stable and the bifurcation is called supercritical; otherwise, it is unstable and the bifurcation is called subcritical.
Because of the complexity of the system and the high number of parameters we are not able to perform a general bifurcation analysis. In this section, we simplify this task by fixing the value of parameters and using m as bifurcation parameter.
In details, we set: With the above choice of the parameters we have and as a consequence z * = 1. In this case the equilibrium is: and x * is the solution of the equation Moreover the functional Jacobian at E * is y * , The characteristic polynomial is The Hurwitz matrix of the characteristic polynomial is given by If s 1 > 0, we always have a negative (real) eigenvalue, while if p (λ) > 0, that is s 2 1 − 3s 2 < 0, we ensure that the other two eigenvalues are complex conjugate. If s 1 s 2 − s 3 > 0 (resp. < 0) then we have at least two eigenvalues with negative (resp. positive) real part.
From Hurwitz-Routh criterion we obtain that a necessary condition for Hopf Bifurcation in this case becomes We have numerically solved the previous equation by using the software Matlab and we obtained the critical value m H = 0.2617.
For this value we obtain s 1 = 0.215694, s 2 = −0.0410984, s 3 = −0.0306826 and s 2 1 − 3s 2 = 0.169819. For m = m H the eigenvalues of the Jacobian matrix J(E * ) are We will see below that as m passes trough the value m = m H the real part of the eigenvalues λ 2,3 change sign from negative to positive.
Then a cycle appears due to Hopf Bifurcation of the equilibrium E * . We are not able to exactly compute the first Lyapunov exponent of the system, however simulations and the sign of the term s 1 s 2 − s 3 suggests that the cycle is stable and the equilibrium looses stability. In Figure 8   In order to illustrate the results of this section we present several numerical simulations. We fix the parameters as above and initial data as follows.
In a first numerical experiment we fix m = 2/10 < m H and as we expect solutions converges to the equilibrium (see Figure 9 below) E * = (0.2668, 0.3778, 1).
In this case the eigenvalues of J(E * ) are and have negative real parts. is unstable, the eigenvalues of J(E * ) are We observe that the eigenvalues λ 2,3 now have positive real part and a Hopf bifurcation occurs at m H = 0.2617. Solutions converge to a limit cycle as shown in Figure 10.

Conclusions
In this paper, we have considered a model describing the dynamics of an ecological system with two prey species and a predator species, which is a modification of the model proposed in [1]. In particular, due to the high availability of one of the two prey populations, we supposed that the predator prefer to predate the more available prey population whereas the other one take advantages of it. In previous articles in the literature (see Introduction for details), it has been pointed out the importance of indirect effects in order to describe real cases of coexistence.
We have performed the stability analysis of equilibria and we have made a detailed analysis of the system on the invariant planes, including the study of the existence of Hopf bifurcation at the equilibrium point E 2 . We have proved that through this bifurcation a stable limit cycle appears.
Regarding the existence of positive equilibria, the expression (21) and the verification of one of the conditions (22) or (23) give us all the positive equilibrium points. The expressions of the equilibria as a function of the parameters are too complicated and not easy to handle. Because of this, we have not been able to determine, in general, for which conditions appear none, one or three equilibria. We have obtained sufficient conditions for the existence of at least one positive equilibrium (see Corollary 1). Also, for fixed values of the parameters it is easy to compute the positive equilibria and maybe it would be possible for certain subfamilies on which less parameters are considered.We have found values of the parameters for which there exist one equilibrium point, and others for which there are not any positive equilibria, but numerically, we have not found values for which three positive equilibria exist, and therefore we think that probably this situation is not feasible, although we have not been able to prove it.
We have shown that Hopf bifurcation can occur also at the positive equilibrium E * and as a consequence, coexistence of the three species via the existence of an attracting limit cycle is possible, by taking into account indirect effects of predation. It is worth mentioning that in [1] Hopf bifurcation is obtained only by considering a version of the model with time delay. Furthermore, we have included a detailed discussion about the problem of persistence of the system.
Throughout the paper, several numerical simulations are provided in order to illustrate the theoretical results.
Due to the complexity of the model, we have not been able to perform a complete bifurcation analysis, then this point remains as an open problem and it would be interesting to study it in the future. A further analysis of predator prey models incorporating indirect effects can be done, for example, considering time delay or non autonomous (seasonal) terms.
A final comment for further research: the proposed model in this paper could be easily implemented by using analog devices ( [34]). Therefore, it could be used to perform experiments in this way.
Author Contributions: All the authors have participate equally in all the aspects of this paper: conceptualization, methodology, investigation, formal analysis, writing-original draft preparation, writing-review and editing. All authors have read and agreed to the published version of the manuscript.