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

mathematics Article Dynamics of a Two Prey and One Predator System with Indirect Effect Renato Colucci 1, Érika Diz-Pita 2and M. Victoria Otero-Espinar 2,*   Citation: Colucci, R.; Diz-Pita, É.; Otero-Espinar, M.V. Dynamics of a Two Prey and One Predator System with Indirect Effect. Mathematics 2021, 9, 436. https://doi.org/10.3390/ math9040436 Academic Editor: Luigi Fortuna Received: 31 January 2021 Accepted: 19 February 2021 Published: 22 February 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Dipartimento di Ingegneria Industriale e Scienze Matematiche, Università Politecnica delle Marche, Via Brecce Bianche 1, 60131 Ancona, Italy; r[email protected] 2 Departamento de Estatística, Análise Matemática e Optimización, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain; [email protected] *Correspondence: mvictoria.oter[email protected] Abstract: We study a population model with two preys and one predator, considering a Holling type II functional response for the interaction between first 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. Keywords: predator–prey; Hopf bifurcation; indirect effects of predation 1. 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: ˙ x=rx1−x k−cxz a+αηy+x, ˙ y=y(β−δz), ˙ z=bxz a+αηy+x+γyz −mz, (1) 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 – 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. Mathematics 2021,9, 436. https://doi.org/10.3390/math9040436 https://www.mdpi.com/journal/mathematics Mathematics 2021,9, 436 2 of 22 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 become ˙ x=rx1−x k−cxz a+αηy+x+Lyz, ˙ y=y(β−δz)−Lyz, ˙ z=bxz a+αηy+x+γyz −mz. (2) 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=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 – 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 3we 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 6contains some conclusive remarks. Mathematics 2021,9, 436 3 of 22 2. 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)∈R3: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 ˙ x=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 ttends to infinity. See Figure 1. x z E0 y E2 E Figure 1. The dynamics on the axes. Now, we present some considerations about the dynamics on the invariant planes. On the invariant plane, z=0 the system is ˙ x=rx1−x k, ˙ y=βy,(3) 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 is ˙ x=rx1−x k−cxz a+x, ˙ z=bxz a+x−mz, (4) and the equilibria are (0, 0),(k, 0). If kb >m(k+a) there exists a further equilibrium point with coordinates (m/(b−m),¯ z)where ¯ z=r cba b−m1−ma k(b−m). Mathematics 2021,9, 436 4 of 22 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 , then ˙ 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 E01111 x z EE Figure 2. Phase portrait on the plane y=0 with bk/(a+k) = m. 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 E2are λ1,2 =A1±√A2 A3, (5) where A1=−cmr(k(m−b) + a(m+b)), A2=c2mr4bk(b−m)2(m(a+k)−bk) + mr(k(m−b) + a(m+b))2, A3=2bck(b−m). (6) x z E0E1111 E22 x-nullclines z-nullclines Figure 3. Phase portrait on the plane y=0 with bk/(a+k)>m. Note that A3> 0 by the existence conditions of the equilibrium point. These eigenvalues are complex if A2<0, and in that case, they have a positive real part if m<b(k−a) k+a, (7) in which case, E2is unstable. If m>b(k−a) k+a, (8) Mathematics 2021,9, 436 5 of 22 the real part of the eigenvalues is negative, so E2 is asymptotically stable. In the case with A2> 0, as the determinant of the Jacobian matrix is positive, it is not possible that the eigenvalues have different sign. Then, if A1 is positive, both eigenvalues are positive and E2 is unstable. In other case, from the conditions A1+√A2< 0 and A1−√A2< 0 we would get that A1<−√A2< 0, which is a contradiction. In the same way, we obtain that if A1<0, then both eigenvalues are negative and E2is asymptotically stable. We give here some results about the possible existence of periodic orbits surrounding the equilibrium point E2in the plane y=0. The equilibrium point E2 is a Hopf equilibrium if and only if A1= 0 and A2< 0, it is, when m=b(k−a)/(k+a) . Note that this occurs only for a<k . In general, when a differential system ˙ x=F(x , µ) in Rn has an equilibrium x0 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(x0) of the differential system at the equilibrium. When `1(x0)< 0, the equilibrium x0 is a weak focus of the differential system restricted to the central surface of x0 , associated to the pair of complex eigenvalues, which cross the imaginary axis, and the limit cycle that emerges from x0is stable. In this case, we say that the Hopf bifurcation is supercritical. Theorem 2. The equilibrium E2 of system (4) undergoes a supercritical Hopf bifurcation at m0=b(k−a)/(a+k)>0 . For m<m0 , the system has a unique and stable limit cycle bifurcating from the equilibrium point E2. Proof. We use the results presented on Chapter 3 of [ 29 ] for computing the first Lyapunov coefficient `1 at the equilibrium E2 . At first, to simplify calculation, we introduce in system (4) a new time variable τby dt = (a+x)dτ, obtaining the polynomial system: ˙ x=r kx(k−x)(a+x)−cxz, ˙ z=bxz −m(a+x)z.(9) This system has the positive equilibrium E2=am b−m,−abr(m(a+k)−bk) ck(b−m)2, which is the same as (am/(b−m) , z) with the notation introduced in Section 2. The Jacobian matrix at this equilibrium is A(m) =      −amr(k(m−b) + a(b+m)) k(b−m)2−acm b−m −abr(m(a+k)−bk) ck(b−m)0      and it has eigenvalues µ(m)±ω(m)i, where µ(m) = −amr(k(m−b) + a(b+m)) 2k(b−m)2and ω(m) = s−a2bmr(m(a+k)−bk) k(b−m)2. (10) We get µ(m0) = 0 for m0=b(k−a) a+k, (11) Mathematics 2021,9, 436 6 of 22 which is positive because as we have said before, a neccesary condition for Hopf bifurcation is a<k. Moreover, ω2(m0) = −abr(a−k)(a+k) 4k>0. (12) Therefore, at m=m0 , the equilibrium point E2 has a pair of pure imaginary eigenvalues ±iω(m) and the system has a Hopf bifurcation. The equilibrium is stable for m>m0 and unstable for m<m0 . 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 µ0(m0) = r(a−k)(a+k)2 8abk <0, (13) where 0 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=m0, and then, the equilibrium E2has the expression E2=k−a 2,r(a+k)2 4ck . (14) We translate E2to the origin of coordinates obtainig the system ˙ ε1=−r kε3 1−r(k−a) 2kε2 1−cε1ε2−ck(k−a) 2kε2, ˙ ε2=2ab a+kε1ε2+abr(a+k)2 2ck(a+k)ε1, (15) which can be represented as ˙ ε=Aε+1 2B(ε,ε) + 1 6C(ε,ε,ε), (16) where A=A(m0)and the multilinear functions Band Care given by B(ε,η) =     −r(k−a) kε1η1−c(ε1η2+ε2η1) 2ab a+k(ε1η2+ε2η1)     , C(ε,η,ζ) =    −6r kε1η1ζ1 0  . We need to find two eigenvectors p,qof the matrix Averifying Aq =iωq,ATp=−iωp, and <p,q>=1, as for example q=1 c(a−k)ω  c(a−k) 2 ωi and p=  ω c(a−k) 2i . (17) Mathematics 2021,9, 436 7 of 22 Now, we compute g20 =hp,B(q,q)i=4abk +r(a+k)(a−k) 4kω(a+k)+1 k−ai, g11 =hp,B(q,q)i=r(a−k) 4kω,g21 =hp,C(q,q,q)i=−3r 4kω2, and the first Lyapynov coefficient `1=1 2ω2Re(ig20g11 +ωg21) = −1 4rωk3(a+k)2 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 E2 through a Hopf Bifurcation for m<m0 . Now, we include some numerical experiments. We fix the parameters as follows: r=k=1, a=0.9, b=1.5, c=1. In this case, m0=0.1667 and the eigenvalues of J(E2)are λ1,2 =±0.253229 i. In Figure 4, we represent the case m= 0.18 >m0 , in which the equilibrium is locally asymptotically stable. In Figure 5, we represent the case m= 0.14 <m0 , in which the equilibrium loses stability and a limit cycle arises due to Hopf bifurcation. 0 100 200 300 400 500 600 700 800 900 1000 Time t 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x(t) z(t) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x(t) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 z(t) Figure 4. The time histories and the solution for m= 0.18 >m0 . In this case, the equilibrium is locally asymptotically stable and nearby solutions converge to it. 0 100 200 300 400 500 600 700 800 900 1000 Time t 0 0.5 1 1.5 2 2.5 x(t) z(t) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 x(t) 0 0.5 1 1.5 2 2.5 z(t) Figure 5. The time histories and the solution for m= 0.14 <m0 . In this case, the equilibrium is no more locally asymptotically stable and a stable limit cycle attracts nearby solutions. Mathematics 2021,9, 436 8 of 22 We conclude this section by considering a case in which there are no periodic orbits in {y=0}∩R3 +: Theorem 3. If r1−a k+bk −m(k+a) k<0, (18) then system (4)does not admit periodic orbits in the set {(x,z)∈R2:x,z≥0}. Proof. Let f(x,z) = rx1−x k−cxz a+xand g(x,z) = bxz a+x−mz. 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 ∆(x,z) = ∂(ϕf) ∂x+∂(ϕg) ∂z 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: ∆(x,z) = r−a kr+b−m−2r kx−am x. We observe that, since ˙ x<0 for x≥k, there are no periodic orbits in the set {(x,z)∈R2 +:x≥k}, 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≤kfor which we obtain ∆(x,z)<r−a kr+b−m−am k. Then, ∆(x , z)< 0 in {(x , z)∈R2 +:x≤k} if r−a kr+b−m−am k< 0, and we conclude that there are no periodic orbits in the whole set {(x,z)∈R2 +}. Remark 1. We observe that, for a<k , the condition of Theorem 3on the nonexistence of periodic orbits r1−a k+b<ma+k k, that is m>rk−a k+a+bk k+a, implies m>bk−a k+a, and then is compatible with the results of bifurcation analysis. 3. 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. Theorem 4. System (2)has the following boundary equilibria on ∂R3 +: •E0= (0,0,0)for any values of the parameters, •E1= (k,0, 0)for any values of the parameters, Mathematics 2021,9, 436 9 of 22 •If kb >m(a+k), the equilibrium E2= (ma/(b−m),0, ¯ z)with ¯ z=r cba b−m1−ma k(b−m). 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 βcm −b(δ+L)rx∗1−x∗ k>0. (19) where x∗is a solution of the equation C4x4+C3x3+C2x2+C1x+C0=0. (20) with the coefficients Ci,i=0,...,4 defined below. Proof. By direct calculation, we obtain the equilibrium E∗= (x∗,y∗,z∗), with z∗=β δ+L, y∗=mc γc+bL −(δ+L)b β(γc+bL)rx∗1−x∗ k, and x∗is the solution of the equation C4x4+C3x3+C2x2+C1x+C0=0. (21) We must require that x∗verifies βcm −b(δ+L)rx∗1−x∗ k>0 so that y∗is positive. This is always verified if kbr δ+L β<4cm. (22) Otherwise, this is satisfied when x∗∈ 0, k 2−sk2 4−kβcm br(δ+L)!∪ k 2+sk2 4−kβcm br(δ+L),∞!. (23) The expressions of the coefficients in Equation (21) are given by Mathematics 2021,9, 436 16 of 22 Proof. By hypothesis the only invariant sets of ∂R3 + are equilibria. The analysis of the previous sections is sufficient to exclude the existence of cycle between equilibria. In Figure 6below we represent the case in which there are only two boundary equilibria (E0 and E1 ) while in Figure 7we represent the case in which we have also E2 . In the latter situation we distinguish two cases: both eigenvalues with positive and both with negative real part respectively. x z E0 y E Figure 6. Possible connections in the cases in which there exist only two boundary equilibria. x z E0 y E2 E0 E0 E1 x z E0 y E2 E Figure 7. Possible connections in the cases in which there exists three boundary equilibria. In the first case the eigenvalues λ2,3 of J(E2) both have positive real part while in the second they have negative real part. In conclusion, we are not able to prove uniform persistence, however the above results guarantee a sort of weak persistence of the three species. 5. 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 mas bifurcation parameter. In details, we set: a=b=c=α=η=r=1, k=β=1 2,δ=4 10,γ=L=1 10. Mathematics 2021,9, 436 17 of 22 With the above choice of the parameters we have β=δ+L, and as a consequence z∗=1. In this case the equilibrium is: E∗= (x∗,y∗,1) where y∗=5hm−x∗+2(x∗)2i and x∗is the solution of the equation D4x4+D3x3+D2x2+D1x+D0=0, with D4=10 >0, D3=−9<0, D2=3>0, D1=−1 2(m−1),D0=−1 2m(1+5m). Moreover the functional Jacobian at E∗is J(E∗) =    f1xf1yf1z 0 0 −1 2y∗ f3xf3y0   where f1x=1−4x∗−(1+y∗) (1+x∗+y∗)2,f1y=x∗ (1+x∗+y∗)2+1 10,f1z=−x∗ 1+x∗+y∗+1 10y∗, f3x=(1+y∗) (1+x∗+y∗)2,f3y=−x∗ (1+x∗+y∗)2+1 10. The characteristic polynomial is p(λ) = λ3+s1λ2+s2λ+s1 where s1=−f1x=−Tr(J(E∗)),s2=−f1zf3x+y∗(δ+L)f3y,s3−y∗(δ+L)hf1xf3y−f1yf3xi. The Hurwitz matrix of the characteristic polynomial is given by H(p) =     1s20 s1s30 s1s2−s3 s10 0     . If s1> 0, we always have a negative (real) eigenvalue, while if p0(λ)> 0, that is s2 1− 3 s2< 0, we ensure that the other two eigenvalues are complex conjugate. If s1s2−s3> 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 s1s2−s3=f3xf1xf1z−1 2y∗f1y=0, Mathematics 2021,9, 436 18 of 22 that is f1xf1z−1 2y∗f1y=0. We have numerically solved the previous equation by using the software Matlab and we obtained the critical value mH=0.2617. For this value we obtain s1= 0.215694, s2=− 0.0410984, s3=− 0.0306826 and s2 1−3s2=0.169819. For m=mHthe eigenvalues of the Jacobian matrix J(E∗)are λ1=−0.9395, λ2,3 =±0.2027i. We will see below that as m passes trough the value m=mH 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 s1s2−s3 suggests that the cycle is stable and the equilibrium looses stability. In Figure 8below we represents the solutions for m=mH , as we expect, they converge to a periodic solution of period T=2π 0.2027 =30.9975. 615 620 625 630 635 640 645 650 Time t 0.3 0.35 0.4 0.45 0.5 x(t) 680 685 690 695 700 705 710 715 Time t 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 y(t) 660 665 670 675 680 685 690 Time t 0.7 0.8 0.9 1 1.1 1.2 1.3 z(t) Figure 8. A limit cycle arises for m=mH . We have represented the solution and a zoom showing the period ( T= 30.9975) of its three components. 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. x(0) = y(0) = z(0) = 1 10. Mathematics 2021,9, 436 19 of 22 In a first numerical experiment we fix m= 2 / 10 <mH and as we expect solutions converges to the equilibrium (see Figure 9below) E∗= (0.2668,0.3778, 1). In this case the eigenvalues of J(E∗)are λ1=−0.5252, λ2,3 =−0.0257 ±0.1898i, and have negative real parts. 0 100 200 300 400 500 600 700 800 Time t 0 0.2 0.4 0.6 0.8 1 1.2 x(t) 0 100 200 300 400 500 600 700 800 Time t 0 1 2 3 4 5 6 7 8 9 10 y(t) 0 100 200 300 400 500 600 700 800 Time t 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 z(t) Figure 9. Graphic and time series of the solution for m= 2 / 10 <mH . The positive equilibrium E∗ is locally stable and nearby solutions converge to it. In a second numerical experiment, we set m= 4 / 10 >mH . The positive equilibrium point E∗= (0.595,2.57, 1) is unstable, the eigenvalues of J(E∗)are λ1=−1.6157, λ2,3 =0.0136 ±0.3237i. We observe that the eigenvalues λ2,3 now have positive real part and a Hopf bifurcation occurs at mH=0.2617. Solutions converge to a limit cycle as shown in Figure 10. Mathematics 2021,9, 436 20 of 22 0 100 200 300 400 500 600 700 800 Time t 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x(t) 0 100 200 300 400 500 600 700 800 Time t 0 5 10 15 20 25 y(t) 0 100 200 300 400 500 600 700 800 Time t 0 1 2 3 4 5 6 z(t) Figure 10. Graphic and time series of the solution for m= 4 / 10 >mH . The positive equilibrium E∗ is unstable and solutions converge to a stable limit cycle. 6. 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 E2 . 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 Mathematics 2021,9, 436 21 of 22 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. Funding: The first author has been partially supported by G.N.A.M.P.A.—INdAM (Italy) and MIUR (Italy). The second and third authors are partially supported by the Ministerio de Economía, Industria y Competitividad, Agencia Estatal de Investigación (Spain), grant MTM2016-79661-P (European FEDER support included, UE) and the Consellería de Educación, Universidade e Formación Profesional (Xunta de Galicia), grant ED431C 2019/10 with FEDER funds. The second author is also supported by the Ministerio de Educacion, Cultura y Deporte de España, contract FPU17/02125. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Conflicts of Interest: The authors declare no conflict of interest. References 1. Sharma, S.; Samanta, G.P. Dynamical Behaviour of a Two Prey and One Predator System. Differ. Eq. Dyn. Syst. 2014 ,22, 125–145. [CrossRef] 2. Bolker, B.; Holyoak, M.; Krivan, V.; Rowe, L.; Schmitz, O. Connecting theoretical and empirical studies of trait-mediated interactions. Ecology 2003,84, 1101–1114. [CrossRef] 3. Cariveau, D.; Irwin, R.E.; Brody, A.K.; Garcia-Mayeya, L.S.; von der Ohe, A. Direct and indirect effects of pollinators and seed predators to selection on plant and floral traits. Oikos 2004,104, 15–26. [CrossRef] 4. Gomez, J.; Zamora, R. Top-down effects in a tritrophoc system: Parasitoids enhance plant fitness. Ecology 1994 ,75, 1023–1030. [CrossRef] 5. Hessen, D.O.; Andersen, T.; Brettum, P.; Faafeng, B.A. Phytoplankton contribution to sestonic mass and elemental ratios in lakes: Implications for zooplankton nutrition. Limnol. Oceanogr. Meth. 2003,48, 1289–1296. [CrossRef] 6. Lundgren, V.; Granéli, E. Grazer-induced defense in Phaeocystis globosa (Prymnesiophyceae): influence of different nutrient conditions. Limnol. Oceanogr. Meth. 2010,55, 1965–1976. [CrossRef] 7. Margalef, R. Life forms of Phytoplanktos as survival alternative in an unstable environment. Oceanol. Acta 1978,134, 493–509. 8. Menge, B.A. Indirect effects in marine rocky intertidal interaction webs: Patterns and importance. Ecol. Monogr. 1995 ,65, 21–74. [CrossRef] 9. Sarnelle, O. Daphnia as keystone predators: Effects on phytoplankton diversity and grazing resistance. J. Plankton Res. 2005 ,27, 1229–1238. [CrossRef] 10. Snyder, W.E.; Ives, A.R. Generalist predators disrupt biological control by a specialist parasitoid. Ecology 2001 ,82, 705–716. [CrossRef] 11. Walsh, M.R.; Reznick, D.N. Interactions between the direct and indirect effects of predators determine life history evolution in a killifish. Proc. Natl. Acad. Sci. USA 2008,105, 594–599. [CrossRef] 12. Wootton, J.T. Indirect effects, prey susceptibility, and habitat selection: Impacts of birds on limpets and algae. Ecology 1992 ,73, 981–991. [CrossRef] 13. Estes, J.; Crooks, K.; Holt, R. Ecological Role of Predators. Enciclopedia Biodivers. 2001,4, 857–878. Mathematics 2021,9, 436 22 of 22 14. Carusela, M.F.; Momo, F.R.; Romanelli, L. Competition, predation and coexistence in a three trophic system. Ecol. Model. 2009 ,220, 2349–2352. [CrossRef] 15. Colucci, R. Coexistence in a one-predator, two-prey system with indirect effects. J. Appl. Math. 2013,2013, 625391. [CrossRef] 16. Smith, H.L.; Thieme, H.R. Dynamical Systems and Population Persistence; American Mathematical Society: Providence, RI, USA, 2011. 17. Butler, G.; Freedman, H.I.; Waltman, P. Uniformly persistent systems. Proc. Am. Math. Soc. 1968,96, 425–430. [CrossRef] 18. Colucci, R.; Nuñez, D. Periodic orbits for a three-dimensional biological differential systems. Abstr. Appl. Anal. 2013 , 2013, 465183. [CrossRef] 19. Caraballo, T.; Colucci, R.; Han, X. Non-autonomous dynamics of a semi-Kolmogorov population model with periodic forcing. Nonlinear Anal. Real World Appl. 2016,31, 661–680. [CrossRef] 20. Caraballo, T.; Colucci, R.; Han, X. Semi-Kolmogorov models for predation with indirect effects in random environments. Discret. Contin. Dyn. Syst. 2016,21, 2129–2143. [CrossRef] 21. Caraballo, T.; Colucci, R.; Han, X. Predation with indirect effects in fluctuating environments. Nonlinear Dyn. 2016 ,84, 115–126. [CrossRef] 22. Poincaré, H. Mémoire sur les courbes définies par une équation différentielle. J. Math. Pures Appl. 1881,7, 375–422. 23. Hilbert, D. Mathematische probleme. Lecture, Second Internat Congr Math. Paris, 1900. Bull. Am. Math. Soc. 1902 ,8, 437–479. [CrossRef] 24. Ilyashenko, Y. Centennial history of Hilbert’s 16th problem. Bull. Am. Math. Soc. 2002,39, 301–354. [CrossRef] 25. Li, J. Hilbert’s 16th problem and bifurcations of planar polynomial vector fields. Int. J. Bifurc. Chaos 2003,13, 47–106. [CrossRef] 26. van der Pol, B. The London, Edinburgh and Dublin Philosophical Magazine; Taylor and Francis: London, UK, 1926; pp. 978–992. 27. Diz-Pita, E.; Llibre, J.; Otero-Espinar, M.V.; Valls, C. The zero-Hopf bifurcations in the Kolmogorov systems of degree 3 in R3 . Commun. Nonlinear Sci. Numer. Simul. 2021,95C. [CrossRef] 28. Han, M.; Llibre, J.; Tian, Y. On the Zero-Hopf Bifurcation of the Lotka-Volterra Systems in R3 .Mathematics 2020 ,8, 1137. [CrossRef] 29. Kuznetsov, Y. Elements of Applied Bifurcation Theory, 2nd ed.; Springer: New York, NY, USA, 1998. 30. Freedman, H.I.; Waltman, P. Mathematical analysis of some three-species food chain models. Math. Biosci. 1977 ,33, 257–276. [CrossRef] 31. Schuster, P.; Sigmund, K.; Wolff, R. Dynamical systems under constant organization III. Cooperative and competitive behaviour in hypercycles. J. Differ. Eq. 1979,32, 357–368. [CrossRef] 32. Hofbauer, J. A general cooperation theorem for hypercycles. Monatshefte Math. 1981,91, 233–240. [CrossRef] 33. May, R.M.; Leonard, W.J. Nonlinear Aspects of Competition Between Three Species. SIAM J. Appl. Math. 1975 ,29, 243–253. [CrossRef] 34. Buscarino, A.; Fortuna, L.; Frasca, M. Essentials of Nonlinear Circuit Dynamics with MATLAB ® and Laboratory Experiments; CRC Press: Boca Raton, FL, USA. [CrossRef]