This is a preprint of: “Hopf periodic orbits for a ratio-dependent predator-prey model with stage structure”, Jaume Llibre, Claudio Vidal, Discrete Contin. Dyn. Syst. Ser. B, vol. 21(6), 1859– 1867, 2016. DOI: [10.3934/dcdsb.2016026] DISCRETE AND CONTINUOUS doi:10.3934/dcdsb.2016017 DYNAMICAL SYSTEMS SERIES B Volume 21, Number 6, August 2016 pp. X–XX HOPF PERIODIC ORBITS FOR A RATIO–DEPENDENT PREDATOR–PREY MODEL WITH STAGE STRUCTURE Jaume Llibre Departament de Matem`atiques, Universitat Aut` onoma de Barcelona 08193 Bellaterra, Barcelona, Catalonia, Spain Claudio Vidal Grupo de Investigaci´ on en Sistemas Din´amicos y Aplicaciones-GISDA Departamento de Matem´ atica, Universidad del Bio Bio Concepci´on, Avda. Collao 1202, Chile (Communicated by Shigui Ruan) Abstract. A ratio–dependent predator-prey model with stage structure for prey was investigated in [8]. There the authors mentioned that they were unable to show if such a model admits limit cycles when the unique equilibrium point E ∗ at the positive octant is unstable. Here we characterize the existence of Hopf bifurcations for the systems. In particular we provide a positive answer to the above question showing for such models the existence of small–amplitude Hopf limit cycles being the equilibrium point E ∗ unstable. 1. Introduction and statement of the main results. Frequently it is assumed in the ratio–dependent predator-prey models that each individual prey admits the same risk to be attacked by a predator. This hypothesis is not always realistic for many species. Thus there are many species whose individuals along his life pass through two stages, immature and mature. Here the prey individuals are classified in immature or mature, and we assume that the immature ones cannot be attacked by the predators. This assumption is reasonable for many mammals, because the immature preys concealed in a mountain cave, are raised by their parents and they do not necessarily go out for seeking food, so the possibility of being attacked by the predators is negligible. Stage structured models have been studied with attention in these last years. Thus a stage–structured model of single species growth with of immature and ma- ture individuals was stated and analyzed in [1]. Later on in [2] it was also assumed that the time from immaturity to maturity is itself state dependent. More recently, in the articles [4, 5, 6, 7, 9] the authors considered predator–prey models with stage structure for prey or predator in order to analyze the influence of a stage structure for the prey or the predator. Xu, Chaplain and Davidson [8] studied the effect of 2010 Mathematics Subject Classification. Primary: 34D23, 92D25. Key words and phrases. Ratio–dependence, predator-prey model, Hopf bifurcation, averaging theory. 1 http://www.gsd.uab.cat 2 JAUME LLIBRE AND CLAUDIO VIDAL stage structure for prey on the dynamics of ratio–dependent predator-prey system by considering the differential system x˙ = −(b + r1 )x + ay, a1 yz y˙ = bx − b1 y 2 − , y + mz (1)   a2 y z˙ = z −r + , y + mz where x = x(t) represents the density of immature individual preys at time t, and y = y(t) denotes the density of mature individual preys at time t, z = z(t) represents the density of the predator at time t. As usual the eight parameters a, a1 , a2 , b, b1 , m, r and r1 of the system are positive, see [8] for their meaning. When y ∗ > 0 and a2 > r the differential system (1) has a unique equilibrium point   ay ∗ (a2 − r)y ∗ aa2 bm − a1 (a2 − r)(b + r1 ) E∗ = , y∗ , with y ∗ = , b + r1 mr a2 b1 m(b + r1 ) in the positive octant, i.e., in {(x, y, z) ∈ R3 : x > 0, y > 0, z > 0}. In the conclusions of the paper [8] the authors said: We would like to mention here that we are unable to show system (1) admits limit cycles when E ∗ exists and is unstable. We leave this for future work. As far as we know this question remained open. The main difficulty for solving it was that the differential system (1) depends on 8 parameters. Here we shall characterize the existence of Hopf bifurcations for systems (1). Moreover using the Hopf bifurcation for these systems we provide a positive answer to the above question restricted to the small–amplitude Hopf limit cycles. The following result characterizes when the equilibrium point E ∗ has eigenvalues of the form ρ, ε ± ωi. Proposition 1. The equilibrium point E ∗ of the differential system (1) has eigen- values of the form ρ, ε ± ωi if and only if β(b + r1 )(a22 r2 − 2a2 r3 + r4 + 2a22 rε − 2a2 r2 ε + a22 ε2 + a2 rε2 + a22 ω 2 + a2 rω 2 ) a= , αb βa2 m(a22 r2 − 2a2 r3 + r4 + 2a22 rε − 2a2 r2 ε + a22 ε2 + a22 ω 2 ) a1 = , α(a2 − r) where α= −a22 br2 + a2 b2 r2 + 2a2 br3 − b2 r3 − br4 − a22 r2 r1 + 2a2 br2 r1 + 2a2 r3 r1 −2br3 r1 − r4 r1 + a2 r2 r12 − r3 r12 − 2a22 brε + 4a2 br2 ε − 2br3 ε − 2a22 rr1 ε +4a2 r2 r1 ε − 2r3 r1 ε − a22 bε2 − a2 brε2 + 2a2 r2 ε2 − 2r3 ε2 − a22 r1 ε2 −a2 rr1 ε2 − a22 bω 2 − a2 brω 2 + 2a2 r2 ω 2 − 2r3 ω 2 − a22 r1 ω 2 − a2 rr1 ω 2 , β= b2 + 2br1 + r12 + 2bε + 2r1 ε + 2ε2 + 2ω 2 . Proposition 1 is proved in section 2. Hypothesis. In the rest of the paper we assume that the parameters a and a1 are the ones given in the statement of Proposition 1. Since a = a(b, m, r, r1 , ε, ω) and a1 = a1 (a2 , b, m, r, r1 , ε, ω), the parameters a2 , b, m, r, r1 , ε and ω must be chosen in such a way that a > 0 and a1 > 0, later on we shall show that such election of the parameters exist. Note that in fact we have changed the parameters a and a1 by the parameters ε and ω. HOPF PERIODIC ORBITS FOR A PREDATOR–PREY MODEL 3 The linearization of the differential system (1) at E ∗ when ε = 0 has the pair of conjugate purely imaginary eigenvalues ±ωi, and if the other one eigenvalues ρ is non–zero this is the setting for a Hopf bifurcation. We can expect to see a small–amplitude limit cycle bifurcating from the equilibrium point E ∗ . It remains to compute the first Liapunov coefficient `1 (E ∗ |ε=0 ) of (1) at E ∗ when ε = 0 and see that it is non–zero. When `1 (E ∗ |ε=0 ) < 0 the point E ∗ is a weak focus on the corresponding two–dimensional central manifold of the system (1) and the limit cycle that borns from E ∗ is stable restricted to this central manifold, consequently the weak focus becomes unstable restricted to the central surface. Due to the fact that the eigenvalues of E ∗ are ρ and ε ± ωi where the small–amplitude limit cycle which appears in the Hopf bifurcation is stable, it follows that ε > 0. In this case we say that the Hopf bifurcation is supercritical. When `1 (E ∗ |ε=0 ) > 0 the point E ∗ is also a weak focus of system (1) but the limit cycle that borns from E ∗ is unstable, and the weak focus becomes stable, and consequently ε < 0. In this second case we say that the Hopf bifurcation is subcritical. Here we use the results on the Hopf bifurcation stated in Chapters 3 and 5 of the book of Kuznetsov [3] for computing `1 (E ∗ |ε=0 ). We have computed in function of all parameters of the system the Liapunov coefficient `1 (E ∗ |ε=0 ), but his expression is so huge that we shall need tens of pages for writing it. All the steps of the computation of `1 (E ∗ |ε=0 ) are described in section 2. These steps as we shall see can be used for computing effectively `1 (E ∗ |ε=0 ) for a given set of parameters. Using them we have proved the following result. Theorem 2. The following system ratio–dependent predator–prey model with stage structure   7y ε2 + 14ε + 50 6ε2 + 4ε + 7 x˙ = −7x + , −40ε2 − 14ε + 2   2 2z 2ε2 + 14ε + 51 4ε2 + 4ε + 5 y y˙ = x − y + , (2) (y + z) (40ε2 + 14ε − 2)   2y z˙ = z −1 + , y+z has a Hopf bifurcation at the equilibrium point E ∗ = (x∗ , y ∗ , z ∗ ) where     ∗ −2ε2 + 24ε + 95 ε2 + 1 ε2 + 14ε + 50 6ε2 + 4ε + 7 16625 x = 2 = + O(ε), 2 4 (20ε + 7ε − 1) 2   ∗ ∗ −2ε2 + 24ε + 95 ε2 + 1 95 y =z = 2 = + O(ε), 2 (−20ε − 7ε + 1) 2 when ε = 0. The eigenvalues at E ∗ are  7 2ε2 − 24ε − 95 665 − 2 =− + O(ε) and ε ± i. 4 (20ε + 7ε − 1) 4 Since `1 (E ∗ |ε=0 ) = −5796656446586/592156593152863031275 < 0 the Hopf bifur- cation is supercritical and an stable Hopf limit cycle on the central manifold through E ∗ there exist for ε > 0 sufficiently small. So locally the equilibrium point E ∗ is unstable on the central manifold when it exhibits the small–amplitude Hopf limit cycle. Theorem 2 is proved in Section 4. We also have the following result. 4 JAUME LLIBRE AND CLAUDIO VIDAL Proposition 3. Assume that a Hopf limit cycle borns at the equilibrium point E ∗ of the differential system (1), then the real eigenvalue at this equilibrium is negative. Proposition 3 is also proved in Section 4. 2. Proof of Proposition 1. The characteristic polynomial of the linear part of the differential system (1) at the equilibrium point E ∗ is p(λ) = λ3 + d2 λ2 + d1 λ + d0 , where 1 d2 = − a1 a22 b + 2aa22 bm + a22 b2 m + a22 bmr + a1 br2 − a2 bmr2 a22 m(b + r1 )  −a1 a22 r1 + 2a22 bmr1 + a22 mrr1 + a1 r2 r1 − a2 mr2 r1 + a22 mr12 , 1 d1 = − a1 a22 b2 + aa22 b2 m − a1 a22 br + 2aa22 bmr + a22 b2 mr a22 m(b + r1 ) +2a1 a2 br2 + a1 b2 r2 − 2aa2 bmr2 − a2 b2 mr2 − a1 br3 − 2a1 a22 br1 +aa22 bmr1 − a1 a22 rr1 + 2a22 bmrr1 + 2a1 a2 r2 r1 + 2a 2 2  1 br r1 − 2a2 bmr r1 3 2 2 2 2 2 2 2 2 −a1 r r1 − a1 a2 r1 + a2 mrr1 + a1 r r1 − a2 mr r1 , r(−a2 + r)(−a1 a2 b + aa2 bm + a1 br − a1 a2 r1 + a1 rr1 ) d0 = − . a22 m Forcing that the characteristic polynomial p(λ) be equal to (λ − ρ)(λ − ε − ωi)(λ − ε + ωi), we get the expressions for a and a1 given in the statement of Proposition 1. Moreover (a2 − r)r(b + r1 ) ρ= − a2 b2 r + a22 r2 − 2a2 r3 + r4 − 2a2 brr1 − a2 rr12 αa2  +2a22 rε − 2a2 brε − 2a2 r2 ε − 2a2 rr1 ε + a22 ε2 − a2 rε2 + a22 ω 2 − a2 rω 2 . Hence, the proposition follows. 3. The computation of the Liapunov coefficient. The following result comes from [3]. Theorem 4. Let x˙ = F (x) be a differential system in Rn having E ∗ as an equilib- rium point. Consider the third order Taylor approximation of F around E ∗ given 1 1 by F (x) = Ax + B(x, x) + C(x, x, x) + O(|x|4 ). Assume that A has a pair of 2! 3! purely imaginary eigenvalues ±ωi. Let q be the eigenvector of A corresponding to the eigenvalue ωi, normalized so that q.q = 1, where q is the conjugate vector of q. Let p be the adjoint eigenvector such that AT p = −ωip and p.q = 1. If I denotes the n × n identity matrix, then `1 (E ∗ |ε=0 ) is equal to 1 Re(p · C(q, q, q) − 2p · B(q, A−1 B(q, q)) + p · B(q, (2ωiI − A)−1 B(q, q))). (3) 2ω The linear part of system (1) at the equilibrium E ∗ is   −b − r1 a 0  a (a2 − r 2 )(b + r1 ) − 2aa22 bm a1 r2   b 1 2 − 2  A=  2 a2 m(b + r1 ) a2 ,   (a2 − r)2 r(r − a2 )  0 a2 m a2 HOPF PERIODIC ORBITS FOR A PREDATOR–PREY MODEL 5 and its eigenvalues are ε ± ωi and ρ = −δ(a2 − r)r(b + r1 )/(σa2 ), where δ = r4 − 2a2 r3 + a22 r2 − 2a2 εr2 − a2 b2 r − a2 r12 r − a2 ε2 r − a2 ω 2 r − 2a2 br1 r+ 2a22 εr − 2a2 bεr − 2a2 r1 εr + a22 ε2 + a22 ω 2 , σ = br4 + r1 r4 + b2 r3 + r12 r3 + 2ε2 r3 + 2ω 2 r3 − 2a2 br3 − 2a2 r1 r3 + 2br1 r3 +2bεr3 + 2r1 εr3 − a2 b2 r2 − a2 r12 r2 − 2a2 ε2 r2 − 2a2 ω 2 r2 + a22 br2 + a22 r1 r2 −2a2 br1 r2 − 4a2 bεr2 − 4a2 r1 εr2 + a2 bε2 r + a2 r1 ε2 r + a2 bω 2 r + a2 r1 ω 2 r +2a22 bεr + 2a22 r1 εr + a22 bε2 + a22 r1 ε2 + a22 bω 2 + a22 r1 ω 2 . Now we compute the bi– and tri–linear functions B and C of Theorem 4, and we obtain  B(x1 , y1 , z1 , x2 , y2 , z2 ) = B1 , B2 , B3 , where Bi = Bi (x1 , y1 , z1 , x2 , y2 , z2 ) are B1 = 0, b1 B2 = − a1 a32 by1 y2 − aa32 bmy1 y2 − 2a1 a2 br2 y1 y2 + a1 br3 y1 y2 + a1 a32 r1 y1 y2 ηa22 −2a1 a2 r2 r1 y1 y2 + a1 r3 r1 y1 y2 − 2a1 a2 bmr2 y2 z1 + 2a1 bmr3 y2 z1 −2a1 a2 mr2 r1 y2 z1 + 2a1 mr3 r1 y2 z1 − 2a1 a2 bmr2 y1 z2 + 2a1 bmr3 y1 z2 −2a1 a2 mr2 r1 y1 z2 + 2a1 mr3 r1 y1 z2 + a1 bm2 r3 z1 z2 + a1 m2 r3 r1 z1 z2 , b1 r(b + r1 ) 2 B3 = a2 y1 y2 − 2a2 ry1 y2 + r2 y1 y2 − 2a2 mry2 z1 + 2mr2 y2 z1 ηa2  −2a2 mry1 z2 + 2mr2 y1 z2 + m2 r2 z1 z2 , η = a1 a2 b − aa2 bm − a1 br + a1 a2 r1 − a1 rr1 ; and  C(x1 , y1 , z1 , x2 , y2 , z2 , x3 , y3 , z3 ) = C1 , C2 , C3 , where Ci = Ci (x1 , y1 , z1 , x2 , y2 , z2 , x3 , y3 , z3 ) are C1 = 0, a1 b21 m(a2 − r)r2 (b + r1 )2 C2 = − a2 y1 y2 y3 − ry1 y2 y3 + a2 my2 y3 z1 η 2 a22  −3mry2 y3 z1 + a2 my1 y3 z2 − 3mry1 y3 z2 + a2 my1 y2 z3 − 3mry1 y2 z3 , a2 C3 = C2 . a1 Computing the normalized eigenvector q of A when ε = 0, associated to the eigenvalue ωi, we obtain q = (q1 , q2 , q3 ) where im(b + r1 )(ib + ir1 + ω)(a2 r − r2 + ia2 ω)(a22 r2 − 2a2 r3 + r4 + a22 ω 2 + a2 rω 2 ) q1 = , −α0 µb(a2 − r)2 m(a2 r − r2 + ia2 ω) q2 = , µ(a2 − r)2 a q3 = , µ where α0 = α|ε=0 and h m2 (a22 r2 − 2a2 r3 + r4 + a22 ω 2 ) m2 (b + r1 )2 2 µ= 1+ 4 + 2 2 (b + 2br1 + r12 + ω 2 ) (a2 − r) α0 b (a2 − r)4 i1/2 (a22 r2 − 2a2 r3 + r4 + a22 ω 2 )(a22 r2 − 2a2 r3 + r4 + a22 ω 2 + a2 rω 2 )2 . 6 JAUME LLIBRE AND CLAUDIO VIDAL The adjoint eigenvector p = (p1 , p2 , p3 ) is bν p1 = , b + r1 − iω p2 = ν, (mr2 (a2 r − r2 + ia2 ω)(b2 + 2br1 + r12 + 2ω 2 )ν p3 = , α0 (a2 − r) where α0 µ(a2 − r)2 (b + r1 + iω) ν= , 2mων1 with ν1 = −ia22 b3 r2 + ia32 br3 + ia2 b3 r3 − 3ia22 br4 + 3ia2 br5 − ibr6 − 3ia22 b2 r2 r1 +ia32 r3 r1 + 3ia2 b2 r3 r1 − 3ia22 r4 r1 + 3ia2 r5 r1 − ir6 r1 − 3ia22 br2 r12 +3ia2 br3 r12 − ia22 r2 r13 + ia2 r3 r13 − a32 br2 ω + a22 b2 r2 ω + 2a22 br3 ω −a2 b2 r3 ω − a2 br4 ω − a32 r2 r1 ω + 2a22 br2 r1 ω + 2a22 r3 r1 ω − 2a2 br3 r1 ω −a2 r4 r1 ω + a22 r2 r12 ω − a2 r3 r12 ω + ia32 brω 2 − 2ia22 br2 ω 2 + ia2 br3 ω 2 +ia32 rr1 ω 2 − 2ia22 r2 r1 ω 2 + ia2 r3 r1 ω 2 − a32 bω 3 − a22 brω 3 + 2a22 r2 ω 3 −2a2 r3 ω 3 − a32 r1 ω 3 − a22 rr1 ω 3 . In short, all the elements appearing in Theorem 4 are given explicitly for the differential system (1), only remains to compute the Liapunov coefficient `1 (E ∗ |ε=0 ) using formula (3). As we mentioned in the introduction we have computed explicitly `1 (E ∗ |ε=0 ) in function of the coefficients of the system (1), using the algebraic manipulator mathematica, but this expression occupies tens of pages. 4. Proofs of Theorem 2 and Proposition 3. Proof of Theorem 2. We shall use the notations and definitions introduced in sec- tion 3. Then, for the differential system (2) we have the matrix     7 ε2 + 14ε + 50 6ε2 + 4ε + 7  −7 − 2 (20ε2 + 7ε − 1) 0         160ε3 + 642ε2 + 370ε + 635 2ε2 + 14ε + 51 4ε2 + 4ε + 5  A= 1 ,  4 (20ε2 + 7ε − 1) 4 (20ε2 + 7ε − 1)      1 1 0 − 2 2 having the eigenvalues given in the statement of Theorem 2. The normalized eigenvector q of A when ε = 0 associated to the eigenvalue ωi is r r ! 441 + 637i 2 2 q= − √ ,− (1 + 2i), − , 600274 300137 300137 The adjoint eigenvector p = (p1 , p2 , p3 ) is r 4559 − 15988i 300137 p1 = , 11056025 2 r 637 − 4659i 300137 p2 = , 442241 2 r −507705 + 172635i 300137 p3 = . 884482 2 HOPF PERIODIC ORBITS FOR A PREDATOR–PREY MODEL 7 The bi–linear function B(x1 , y1 , z1 , x2 , y2 , z2 ) = (B1 , B2 , B3 ) with Bi = Bi (x1 , y1 , z1 , x2 , y2 , z2 ) is given by B1 = 0, 1 B2 = (−875y1 y2 − 3570y2 z1 − 3570y1 z2 + 1785z1 z2 ), 2660 1 B3 = (−y1 y2 + 2y2 z1 + 2y1 z2 − z1 z2 ); 190 and the tri–linear function C(x1 , y1 , z1 , x2 , y2 , z2 , x3 , y3 , z3 ) = (C1 , C2 , C3 ) with Ci = Ci (x1 , y1 , z1 , x2 , y2 , z2 , x3 , y3 , z3 ) is given by C1 = 0, 51 C2 = − (y1 y2 y3 − y2 y3 z1 − y1 y3 z2 − y1 y2 z3 ), 7220 1 C3 = (y1 y2 y3 − y2 y3 z1 − y1 y3 z2 − y1 y2 z3 ). 18050 According to Theorem 4 for computing `1 (E ∗ |ε=0 ), we need to compute first A−1 and a (2iI − A)−1 . Note that now ω = 1. We have that   89 70 8925 −  133 −  19 19    2 2 51  A−1 =   − 665 − 95 ,  19    2 2 13  − − 665 95 19 and   7 + 2i −1225 0  635 255    (2iI − A)−1 =  −1 + 2i .  4 4    1 1 0 − + 2i 2 2 The first, second and third terms of `1 (E ∗ |ε=0 ) given in Theorem 4 are 950436 + 129948i p · C(q, q, q) = , 239582861065685 −152367568 + 104825176i −2p · B(q, A−1 B(q, q)) = , 4552074360248015 928154260416 + 1952607677788i p · B(q, (2iI − A)−1 B(q, q)) = . 93498409445188899675 Consequently we get 5796656446586 `1 (E ∗ |ε=0 ) = − . 592156593152863031275 This completes the proof of the theorem. Proof of Proposition 3. From the expression of ρ, given in section 2, and for |ε| sufficiently small it follows that ρ > 0 if and only if α0 (−a2 b2 r + a22 r2 − 2a2 r3 + r4 − 2a2 brr1 − a2 rr12 + a22 ω 2 − a2 rω 2 ) > 0. (4) Here we have used that a2 > r, recall that this is due to the fact that the equilibrium point E ∗ must be in the positive octant. 8 JAUME LLIBRE AND CLAUDIO VIDAL By Proposition 1 and for |ε| sufficiently small we get that a > 0 if and only if α0 (a22 r2 − 2a2 r3 + r4 + a22 ω 2 + a2 rω 2 ) > 0. (5) Note that β|ε=0 > 0, see its definition in the statement of Proposition 1. On the other hand, the coordinate x of the equilibrium point E ∗ is positive if and only −(−a2 b2 r + a22 r2 − 2a2 r3 + r4 − 2a2 brr1 − a2 rr12 + a22 ω 2 − a2 rω 2 )· (6) (a22 r2 − 2a2 r3 + r4 + a22 ω 2 + a2 rω 2 ) > 0. Now, since the x–coordinate of E ∗ and the parameter a are positive, the inequal- ities (5) and (6) hold. Then, from (4) we obtain that ρ < 0. That is, the real eigenvalue of the equilibrium point E ∗ is negative for |ε| sufficiently small. Hence the proposition is proved. 5. Conclusions. Xu, Chaplain and Davidson [8] studied the effect of stage struc- ture for prey on the dynamics of ratio–dependent predator-prey system by consid- ering the differential system x˙ = −(b + r1 )x + ay, a1 yz y˙ = bx − b1 y 2 − , y + mz (7)   a2 y z˙ = z −r + , y + mz where x = x(t) represents the density of immature individual preys at time t, and y = y(t) denotes the density of mature individual preys at time t, z = z(t) represents the density of the predator at time t. As usual the eight parameters a, a1 , a2 , b, b1 , m, r and r1 of the system are positive. It is verified that when y ∗ > 0 and a2 > r the differential system (7) has a unique equilibrium point   ∗ ay ∗ ∗ (a2 − r)y ∗ aa2 bm − a1 (a2 − r)(b + r1 ) E = ,y , with y ∗ = , b + r1 mr a2 b1 m(b + r1 ) in the positive octant, i.e., in {(x, y, z) ∈ R3 : x > 0, y > 0, z > 0}. In the conclusions of the paper [8] the authors said: We would like to mention here that we are unable to show system (1) admits limit cycles when E ∗ exists and is unstable. We leave this for future work. As far as we know this question remained open. In this work we answer affirmatively this question, see Theorem 2. The main difficulty for solving it was that the differential system (7) depends on 8 parameters, our main result show the existence of a Hopf bifurcation at the equilibrium point E ∗ for appropriate values of the parameters. Acknowledgments. The first author is partially supported by a MINECO grant MTM2013-40998-P, an AGAUR grant number 2013SGR-568, the grants FP7-PEO- PLE-2012-IRSES 318999 and 316338. The second author is partially supported by Direcci´on de Investigaci´on DIUBB 1204084/R. REFERENCES [1] W. G. Aiello and H. I. Freedman, A time delay model of single–species growth with stage structure, Math. Biosci., 101 (1990), 139–153. [2] W. G. Aiello, H. I. Freedman and J. Wu, Analysis of a model representing stage-structured population growth with state–dependent time delay, SIAM J. Appl. Math., 52 (1992), 855– 869. HOPF PERIODIC ORBITS FOR A PREDATOR–PREY MODEL 9 [3] Y. Kuznetsov, Elements of Applied Bifurcation Theory, Applied Mathematical Sciences, Vol. 112, Springer–Verlag, New York, 2004. [4] Z. Li, M. Han and F. Chen, Global stability of a predator–prey system with stage structure and mutual interference, Discrete Contin. Dyn. Syst. Ser. B, 19 (2014), 173–187. [5] K. G. Magnusson, Destabilizing effect of cannibalism on a structured predator–prey system, Math. Biosci., 155 (1999), 61–75. [6] W. Wang and L. Chen, A predator–prey system with stage structure for predator, Comput. Math. Appl., 33 (1997), 83–91. [7] R. Xu, Global convergence of a predator–prey model with stage structure and spatio–temporal delay, Discrete Contin. Dyn. Syst. Ser. B, 15 (2011), 273–291. [8] R. Xu, M. A. J. Chaplain and F. A. Davidson, Persistence and global stability of a ratio– dependent predator–prey model with stage structure, Appl. Math. Comput., 158 (2004), 729–744. [9] X. Zhang and L. Chen, The stage–structured predator–prey model and optimal harvesting policy, Math. Biosci., 168 (2000), 201–210. Received October 2014; revised March 2016. E-mail address: jllibre@mat.uab.cat E-mail address: clvidal@ubiobio.cl