This is a preprint of: “Final evolutions of a class of May-Leonard Lotka-Volterra systems”, Claudio A. Buzzi, Jaume Llibre, Robson A. T. Santos, J. Nonlinear Math. Phys., vol. 27(2), 267–278, 2020. DOI: [10.1080/14029251.2020.1700635] GLOBAL DYNAMICS OF A CLASS OF MAY-LEONARD LOTKA-VOLTERRA SYSTEMS CLAUDIO A. BUZZI1 , JAUME LLIBRE2 , ROBSON A. T. SANTOS3 Abstract. We study a particular class of Lotka-Volterra 3-dimensional sys- tems called May-Leonard systems, which depend on two real parameters a and b, when a + b = −1. For these values of the parameters we shall describe its global dynamics in the compactification of the positive octant of R3 includ- ing its infinity. This can be done because this differential system possesses a Darboux invariant. 1. Introduction Polynomial ordinary differential systems are often used in various branches of applied mathematics, physics, chemist, engineering, etc. Models studying the in- teraction between species of predator-prey type have been extensively analyzed as the classical Lotka-Volterra systems. For more information on the Lotka-Volterra systems see for instance [6] and the references quoted there. In particular, one of these competition models between three species inside the class of 3-dimensional Lotka-Volterra systems is the May-Leonard model given by the polynomial differ- ential system in R3 ẋ = x(1 − x − ay − bz), (1) ẏ = y(1 − bx − y − az), ż = z(1 − ax − by − z), where a and b are real parameters and the dot denotes derivative with respect to the time t. See for more details on the May-Leonard system the papers [7] and [1], and the references quoted there. The Lotka-Volterra systems in R3 have the property that the three coordinate planes are invariant by the flow of these systems. Moreover, at points of straight line x = y = z, system (1) writes ẋ = x, ẏ = y, ż = z. Therefore, the bisectrix of the positive octant is an invariant straight line for this differential system. In this paper we describe the global dynamics of system (1) in function of the parameters a and b when a + b = −1. The system (1) is defined in R3 . In order to study the dynamics of its orbits at infinity we extend analytically its flow by using the Poincaré compactification of R3 . In the appendix we give precise definitions for this compactification. The region of interest in our study is the positive octant 2010 Mathematics Subject Classification. Primary 34D45, 34D05, 37N25, 92D25, 34C12. Key words and phrases. May-Leonard system, Lotka-Volterra system, invariant, global dynamics. 1 http://www.gsd.uab.cat 2 C.A. BUZZI, J. LLIBRE AND R.A.T. SANTOS of R3 , i.e. where x ≥ 0, y ≥ 0, z ≥ 0. So we shall study the flow of the Poincaré compactification in the region R = {(x, y, z) ∈ R3 : x2 + y 2 + z 2 ≤ 1, x ≥ 0, y ≥ 0, z ≥ 0} of the Poincaré ball. We remark that the global dynamics of the May-Leonard system (1) with a+b = −1 can be study because this differential system has a Darboux invariant. Here with this model we illustrate how can be used a Darboux invariant for obtaining the global dynamics of a differential system. The differential system (1) has been extensively studied in order to understand the interaction between species and try to predict possible extinction or overpopula- tion for example. However our interest is purely mathematical, we want to illustrate how the use of Darboux invariants can be used for describing the global dynamics of a differential system. Note that we are interested in the study of system (1) for all real values of the parameters a and b satisfying a + b = −1, and not only for their positive values. Consequently our analysis has no biological meaning. This study could be made in a similar way in the others octants of R3 . 2. Statement of the main results We denote by X the polynomial vector field associated to the differential system (1), and by p(X) the Poincaré compactification of X, see the appendix in section 5. The flow of system (1) in the region R is described in the next two theorems. For a formal definition of topologically equivalent phase portraits see [4]. Theorem 1. For the May-Leonard differential system (1) in the octant R the following statements hold when a + b = −1. (a) The phase portrait of the Poincaré compactification p(X) of system (1) on the boundaries x = 0, y = 0 and z = 0 of R is topologically equivalent to the one described in Fig. 1(a) if a ≤ −2 or a ≥ 1, and in Fig. 2(a) if −2 < a < 1. (b) The phase portrait of the Poincaré compactification p(X) of system (1) on R∞ = ∂R ∩ {x2 + y 2 + z 2 = 1} (i.e. the phase portrait at the infinity of the positive octant of R3 ) is topologically equivalent to the one described in Fig. 1(b) if a ≤ −2 or a ≥ 1, Fig. 2(b) if a = −1/2, Fig. 2(c) if −2 < a < −1/2, and Fig. 2(d) if −1/2 < a < 1. (c) When a = −1/2 the planes x = y, x = z and y = z are invariant by the flow of system (1), and the phase portrait of the Poincaré compactification p(X) of system (1) on R ∩ {x = y}, R ∩ {x = z} and R ∩ {y = z} are topologically equivalent to the ones described in (a), (b) and (c) of Fig. 3 respectively. Let p(γ) denote the orbit γ of the vector field X associated to system (1) in the Poincaré compactification p(X). Theorem 2. Let γ be an orbit of system (1) with a + b = −1 such that p(γ) is contained in the interior of R. Then for any value a ∈ R the following statements hold. MAY-LEONARD LOTKA-VOLTERRA SYSTEMS 3 (a) The α-limit set of p(γ) is the origin of R3 . (b) The ω-limit set of p(γ) is the infinite singular point {( √13 , √13 , √13 )} ∈ R∞ . z z y y x x (a) (b) Figure 1. The global dynamics on the boundary of R for a + b = −1 and a ≤ −2. (a) The dynamics on xyz = 0. (b) The dynamics on R∞ . Reversing the sense of all the orbits we have the global dynamics on the boundary of R for a + b = −1 and a ≥ 1. 3. Proof of Theorem 1 The finite singular points of system (1) with a + b = −1 are the solutions of the system P1 (x, y, z) = x(1 − x + z + a(−y + z)) = 0, P2 (x, y, z) = y(1 + x − y + a(−z + x)) = 0, P3 (x, y, z) = z(1 + y − z + a(−x + y)) = 0, namely p0 = (0, 0, 0), p1 = (1, 0, 0), p2 = (0, 1, 0), p3 = (0, 0, 1),       1−a 2+a 2+a 1−a 1−a 2+a p4 = 0, , , p5 = , 0, , p6 = , ,0 , A A A A A A 2 where A = 1 + a + a . Since A > 0 for a ∈ R and the region of interest is R, we have: (i) If a ≤ −2 or a ≥ 1 system (1) has only four finite equilibrium points: p0 , p1 , p2 and p3 . (ii) If −2 < a < 1 system (1) has the seven finite equilibrium points pj for j = 0, 1 . . . , 6. All these finite equilibrium points are hyperbolic if a 6= −2, 1, and consequently its local phase portrait is topologically equivalent to the phase portrait of its linear part by the Hartman–Grobman Theorem, see for instance [3]. We note that when a ∈ (−2, 1) and a → 1 we have that p4 → p3 , p5 → p1 and p6 → p2 ; while if a → −2 we have that p4 → p2 , p5 → p3 and p6 → p1 . This behavior of these equilibria allows to determine by continuity the local phase 4 C.A. BUZZI, J. LLIBRE AND R.A.T. SANTOS z z p∞ yz p∞ yz p∞ xz p∞ xz p∞ b y y p∞ p∞ xy xy x x (a) (b) z z p∞ yz p∞ yz p∞ b p∞ xz p∞ xz p∞ b y y p∞ xy p∞ xy x x (c) (d) Figure 2. The global dynamics on the boundary of R for a + b = −1 and −2 < a < 1. (a) The dynamics on xyz = 0. The dynamics on R∞ for a = −1/2 in (b), for a ∈ (−2, −1/2) in (c), and for a ∈ (−1/2, 1) in (d). portraits on the boundary of R of the non-hyperbolic equilibrium points p1 , p2 and p3 when a = −2 and a = 1 from the global phase portraits of the boundary of R when a ∈ (−2, 1). We also observe that when a ∈ (−2, 1) and a → 1 we have that p∞ xz → the positive endpoint of the x-axis, p∞yz → the positive endpoint of the z-axis, and p∞ xy → the positive endpoint of the y-axis; while if a → −2 we have that pxz → ∞ ∞ the positive endpoint of the z-axis, pyz → the positive endpoint of the y-axis, and p∞ xy → the positive endpoint of the x-axis. So the behavior of these equilibria allows to determine by continuity the local phase portraits on the boundary of R of the non-hyperbolic equilibrium which are at the positive endpoints of the axes of coordinates when a = −2 and a = 1 from the global phase portraits of the boundary of R when a ∈ (−2, 1). Hence in what follows we only prove Theorem 1 for the values a ∈ R \ {−2, 1}. The linear part of system (1) at the equilibrium p0 is the identity matrix. There- fore it is a repelling equilibrium. MAY-LEONARD LOTKA-VOLTERRA SYSTEMS 5 z z y y x x (a) (b) z y x (c) Figure 3. The global dynamics on R ∩ {x = y}, R ∩ {x = z} and R ∩ {y = z} respectively, when a = b = −1/2. The eigenvalues of linear part at equilibrium points p1 , p2 and p3 are −1, 1 − a, 2 + a. Therefore when a < −2 or a > 1 these equilibria have a 2-dimensional stable manifold and an 1-dimensional unstable one; and for −2 < a < 1 these equilibria have a 2-dimensional unstable manifold and an 1-dimensional stable one. When −2 < a < 1 the eigenvalues of linear part at equilibrium points p4 , p5 and p6 are 3, −1 and (−2 + a + a2 )/A. Since −2 + a + a2 < 0 for −2 < a < 1 then these equilibria have a 2-dimensional stable manifold and an 1-dimensional unstable one. Moreover, p4 (respectively p5 and p6 ) is an attractor restricted to the invariant boundary x = 0 (respectively y = 0 and z = 0). Now we shall study the infinite equilibrium points. For study the dynamics on the infinity R∞ of R we shall use the Poincaré compactification of the differential system (1). See Appendix for details. Thus the differential system (1) in the local chart U1 becomes z˙1 = 2z1 + az1 − z12 + az12 − z1 z2 − 2az1 z2 , (2) z˙2 = z2 − az2 + z1 z2 + 2az1 z2 − 2z22 − az22 , z˙3 = z3 + az1 z3 − z2 z3 − az2 z3 − z32 . So system (1) has two equilibrium points at infinity: (0, 0, 0) and (1, 1, 0). The linear part at the equilibrium (0, 0, 0) has the eigenvalues 1 − a and 2 + a at infinity 6 C.A. BUZZI, J. LLIBRE AND R.A.T. SANTOS and eigenvalue 1 in its finite direction. Therefore, on the infinity (0, 0, 0) is a saddle such that its stable separatrix is contained in the z1 -axis when a < −2 (respectively z2 -axis when a > 1). √ The eigenvalues of linear part at equilibrium (1, 1, 0) are (−3 ± i 3(1 + 2a))/2 and −1. Therefore, on the infinity (1, 1, 0) is a stable focus turning clockwise if a < −2 (respectively counterclockwise if a > 1). Now system (1) in the local chart U2 writes z˙1 = z1 − az1 − 2z12 − az12 + z1 z2 + 2az1 z2 , (3) z˙2 = 2z2 + az2 − z1 z2 − 2az1 z2 − z22 + az22 , z˙3 = z3 − z1 z3 − az1 z3 + az2 z3 − z32 . Since the local chart U2 covers the end part of the plane x = 0 at infinity of the positive octant of R3 , we are interested only in the equilibrium points which are on z1 = 0 and z3 = 0. In this case, there is one equilibrium point at infinity: (0, 0, 0). The eigenvalues of linear part at equilibrium (0, 0, 0) are 1 − a and 2 + a at infinity and eigenvalue 1 in its finite direction. Therefore, on the infinity (0, 0, 0) is a saddle such that its stable separatrix is contained in the z2 -axis when a < −2 (respectively z1 -axis when a > 1). Now for a ≤ −2 or a ≥ 1 we only need to study the equilibrium point at the endpoint of positive z-half-axis, i.e. the equilibrium point at the origin of the local chart U3 . In the local chart U3 system (1) becomes z˙1 = 2z1 + az1 − z12 + az12 − z1 z2 − 2az1 z2 , z˙2 = z2 − az2 + z1 z2 + 2az1 z2 − 2z22 − az22 , z˙3 = z3 + az1 z3 − z2 z3 − az2 z3 − z32 . The linear part at the equilibrium point (0, 0, 0) has eigenvalues 1 − a and 2 + a at infinity and eigenvalue 1 in its finite direction. Therefore the origin of the local chart U3 is a saddle such that its stable separatrix is contained in the z1 -axis when a < −2 (respectively z2 -axis when a > 1). It remains to study the infinite equilibrium points of system (1) in case −2 < a < 1. In the local chart U1 the system (2) has four equilibrium points at infinity: (0, 0, 0), ((2 + a)/(1 − a), 0, 0), (0, (1 − a)/(2 + a), 0) and (1, 1, 0). The eigenvalues of linear part at equilibrium (0, 0, 0) are 1 − a and 2 + a at infinity and eigenvalue 1 in its finite direction. Then this equilibrium is an unstable node. The linear part at the equilibrium ((2 + a)/(1 − a), 0, 0) (respectively (0, (1 − a)/(2 + a), 0)) has the eigenvalues −(2 + a), A/(1 − a) and 3A/(a − 1) (respectively (a − 1), A/(2 + a) and 3A/(2 + a)). Therefore the equilibria ((2 + a)/(1 − a), 0, 0) and (0, (1 − a)/(2 + a), 0) are saddles such that its stable separatrix is contained in the z1 -axis and z2 -axis respectively. √ The eigenvalues of the linear part at the equilibrium (1, 1, 0) are (−3 ± i 3(1 + 2a))/2 and −1. Therefore, on the infinity (1, 1, 0) is a stable focus turning clockwise if −2 < a < −1/2 (respectively counterclockwise if −1/2 < a < 1). When a = −1/2 on the infinity (1, 1, 0) is stable node. Now since we are interested only in the equilibrium points which are on z1 = 0 and z3 = 0, in the local chart U2 the system (3) has two equilibrium points: (0, 0, 0) and (0, (2 + a)/(1 − a), 0). The eigenvalues of linear part at equilibrium (0, 0, 0) are 1 − a and 2 + a at infinity and eigenvalue 1 in its finite direction. Then, on the infinity this equilibrium is an unstable node. The linear part at the equilibrium MAY-LEONARD LOTKA-VOLTERRA SYSTEMS 7 (0, (2 + a)/(1 − a), 0) has the eigenvalues −(2 + a), A/(1 − a) and 3A/(1 − a) at infinity. Therefore, on the infinity (0, (2 + a)/(1 − a), 0) is a saddle such that its stable separatrix is contained in the z2 -axis. In the local chart U3 we only need to study the equilibrium point (0, 0, 0). The eigenvalues of linear part at this equilibrium are 1 − a and 2 + a at infinity and eigenvalue 1 in its finite direction. Then, on the infinity this equilibrium is an unstable node. So, the proof of statements (a) and (b) of Theorem 1 is complete. Now since the bisectrix x = y = z is an invariant straight line for the system, it is easy to check for a = −1/2 that the global phase portrait on the invariant planes R ∩ {x = y}, R ∩ {x = z} and R ∩ {y = z} are topologically equivalents to the ones described in (a), (b) and (c) of Fig. 3 respectively. This completes the proof of Theorem 1. 4. Proof of Theorem 2 We say that a C 1 function I(x, y, z, t) is an invariant of the polynomial differen- tial system (1) if I(x(t), y(t), z(t), t) is constant, for all the values of t for which the solution (x(t), y(t), z(t)) of (1) is defined. When the function I is independent of the time, then it is called a first integral of differential system (1). Also if an invariant I(x, y, z, t) is of the form f (x, y, z)est , then it is called a Darboux invariant. Proposition 1. System (1) has the Darboux invariant I = I(t, x, y, z) = xyze−3t . Proof. It is immediate to check that dI ∂I ∂I ∂I ∂I = ẋ + ẏ + ż + = 0, dt ∂x ∂y ∂z ∂t where ẋ, ẏ and ż are given in (1). Therefore I is a Darboux invariant of system (1).  For knowing how to obtain the Darboux invariant given in Proposition 1 see statement (vi) of Theorem 8.7 of [4], there the theory is described for polynomial vector fields in R2 , but the results and the proofs extend to R3 . Proposition 2. Let I(x, y, z, t) = f (x, y, z)est be a Darboux invariant of system (1). Let p ∈ R3 and ϕp (t) the solution of system (1) such that ϕp (0) = p. Then α(p), ω(p) ⊂ {f (x, y, z) = 0} ∪ S2 . Here α(p) and ω(p) denote the α-limit and ω-limit sets of p respectively, and S2 denotes the boundary of the Poincaré ball, i.e the infinity of R3 . For a proof of Proposition 2 see [5]. Proof of Theorem 2. Let p(γ) = {ϕp (t) = (x(t), y(t), z(t)) : t ∈ R} be the orbit of the Poincaré compactification of system (1) with a + b = −1 such that ϕp (0) = p with p in the interior of R. We recall that all the orbits of a differential system defined on a compact set are defined for all t ∈ R. By Propositions 1 and 2 the α- and ω-limit set of p(γ) is contained in boundary of R, i.e. in {(x, y, z) ∈ R : 8 C.A. BUZZI, J. LLIBRE AND R.A.T. SANTOS xyz = 0} ∪ S2 . Furthermore, by Proposition 1, I(t, x(t), y(t), z(t)) = k constant with k > 0. So (4) x(t)y(t)z(t) = ke3t , for all t. Taking limit in (4) when t → −∞, we obtain lim x(t)y(t)z(t) = 0. t→−∞ Assume that lim x(t) = 0. If lim y(t) = 0 or lim z(t) = 0 the proof is t→−∞ t→−∞ t→−∞ completely similar. Then the points of α-limit set of p(γ) have coordinate x = 0. We distinguish two cases. Case 1. Suppose that a ≤ −2 or a ≥ 1. Therefore, by Theorem 1(a) and Fig.1(a) we have that α-limit set of p(γ) can be only the equilibrium points p0 , p2 , p3 , p∞ ∞ y (positive endpoint of the y-axis) or pz (the positive endpoint of the z-axis). Except the origin p0 , the remaining equilibrium points are of saddle type with their unstable separatrix of each of them contained in one of the coordinate invariant planes. Thus, as p(γ) is on interior of R, then the α-limit set of p(γ) is the origin of R3 . Case 2. Assume that −2 < a < 1. By Theorem 1(a) and Fig.2(a) we have that the α-limit set of p(γ) can be only one of the equilibrium points p0 , p2 , p3 , p4 , p∞ ∞ y , pz ∞ or pyz . Again p0 is the only of these equilibrium points that is not of saddle type. In summary the α-limit set of p(γ) is the origin of R3 for any value a ∈ R. So statement (a) is proved. Now we study the ω-limit set of p(γ). In a similar way taking limit in (4) when t → +∞ we obtain lim x(t)y(t)z(t) = +∞. t→+∞ Suppose that lim x(t) = +∞. The cases lim y(t) = +∞ or lim z(t) = +∞ t→+∞ t→+∞ t→+∞ are proved in a similar way. So by Proposition 2 we conclude that ω-limit set of p(γ) is contained in R∞ . Again, we distinguish two cases. Case 1. Assume that a ≤ −2 or a ≥ 1. By Theorem 1(b) and Fig.1(b) the α-limit set of p(γ) can be only the positive endpoint of x-axis p∞ x , or the end of invariant bisectrix p∞ ∞ b . However the infinite equilibrium point px is saddle type with stable separatrix contained in R∞ ∩ {z = 0} and the orbit p(γ) is in the interior of R. Therefore the ω-limit set of p(γ) is the infinite singular point p∞ √1 √1 √1 b = {( 3 , 3 , 3 )}. Case 2. Suppose that −2 < a < 1. By Theorem 1(b) and Figs.2(b)(c)(d) we have that ω-limit set of P (γ) can be only the equilibrium points p∞ ∞ ∞ ∞ x , pxy , pxz or pb . In ∞ ∞ ∞ this case the infinite singular point px is a repeller, and pxy , pxz are of saddle type with stable separatrix contained in R∞ ∩ {z = 0} and R∞ ∩ {y = 0} respectively. So the ω-limit set of p(γ) is p∞b because the orbit p(γ) is in the interior of R. In short the ω-limit set of p(γ) is the infinite singular point {( √13 , √13 , √13 )} ∈ R∞ for any value a ∈ R. Hence statement (b) of Theorem 2 is proved.  MAY-LEONARD LOTKA-VOLTERRA SYSTEMS 9 5. Appendix: The Poincaré compactification in R3 For more details on the Poincaré compactification in R3 see [2]. In R3 we consider the polynomial differential system ẋ = P1 (x, y, z), ẏ = P2 (x, y, z), ż = P3 (x, y, z), or equivalently its associated polynomial vector field X = (P1 , P2 , P3 ). The degree n of X is defined as n = max{deg(Pi ) : i = 1, 2, 3}. Let S3 = {y = (y1 , y2 , y3 , y4 ) ∈ R4 : kyk = 1} be the unit sphere in R4 , and H+ = {y ∈ S3 : y4 > 0} and H− = {y ∈ S3 : y4 < 0} be the northern and southern hemispheres, respectively. The tangent space to S3 at the point y is denoted by Ty S3 . Then we identify the tangent hyperplane T(0,0,0,1) S3 = {(x1 , x2 , x3 , 1) ∈ R4 : (x1 , x2 , x3 ) ∈ R3 } with R3 . We consider the central projections f+ : R3 = T(0,0,0,1) S3 → H+ and f− : R3 = T(0,0,0,1) S3 → H− , defined by 1 1 f+ (x) = (x1 , x2 , x3 , 1) and f− (x) = − (x1 , x2 , x3 , 1), ∆x ∆x P3 where ∆x = (1 + i=1 x2i )1/2 . Through these central projection R3 is identified with the northern and the southern hemispheres, respectively. The equator of the sphere S3 is S2 = {y ∈ S3 : y4 = 0}. Clearly S2 can be identified with the infinity of R3 . The diffeomorphisms f+ and f− define two copies of X, one Df+ ◦ X in the northern hemisphere and the other Df− ◦ X in the southern one. Denote by X̄ the vector field on S3 \ S2 = H+ ∪ H− such that restricted to H+ coincides with Df+ ◦ X and restricted to H− coincides with Df− ◦ X. We extend analytically the polynomial vector field X̄ to the equator of S3 , i.e. to the infinity of R3 , in such a way that the flow on the boundary is invariant. This is done defining the vector field p(X)(y) = y4n−1 X̄(y), for all y ∈ S3 . This extended vector field p(X) is called the Poincaré compactifica- tion of X on the Poincaré sphere S3 . In what follows we shall work with the orthogonal projection of the closed north- ern hemisphere to y4 = 0. Note that this projection is a closed ball B of radius one, whose interior is diffeomorphic to R3 and whose boundary S2 corresponds to the in- finity of R3 . The projected vector field on B is called the Poincaré compactification on the Poincaré ball B. As S3 is a differentiable manifold, to compute the expression for p(X) we can consider the eight local charts (Ui , Fi ), (Vi , Gi ) where Ui = {y ∈ S3 : yi > 0} and Vi = {y ∈ S3 : yi < 0} for i = 1, 2, 3, 4; the diffeomorphisms Fi : Ui → R3 and Gi : Vi → R3 for i = 1, 2, 3, 4, are the inverses of the central projections from the 10 C.A. BUZZI, J. LLIBRE AND R.A.T. SANTOS origin to the tangent planes at the points (±1, 0, 0, 0), (0, ±1, 0, 0), (0, 0, ±1, 0) and (0, 0, 0, ±1), respectively. The expression of p(X) on the local chart U1 is (5) z3n (−z1 P1 + P2 , −z2 P1 + P3 , −z3 P1 ), where Pi = Pi (1/z3 , z1 /z3 , z2 /z3 ), and the expressions of p(X) in U2 is (6) z3n (−z1 P2 + P1 , −z2 P2 + P3 , −z3 P2 ), where Pi = Pi (z1 /z3 , 1/z3 , z2 /z3 ) in U2 , and in U3 is z3n (7) (−z1 P3 + P1 , −z2 P3 + P2 , −z3 P3 ), (∆z)n−1 where Pi = Pi (z1 /z3 , z2 /z3 , 1/z3 ) in U3 . The expression for p(X) in U4 is z3n+1 (P1 , P2 , P3 ) where the component Pi = Pi (z1 , z2 , z3 ). The expression for p(X) in the local chart Vi is the same as in Ui multiplied by (−1)n−1 . We remark that all the points on the sphere at infinity in the coordinates of any local chart have z3 = 0. Acknowledgments The first author is supported by FAPESP grant 2013/2454-1, CAPES grant 88881.068462/2014-01 and EU Marie-Curie IRSES Brazilian-European partner- ship in Dynamical Systems (FPT-PEOPLE-2012-IRSES 31899 BREUDS). The second author is partially supported by the MINECO grants MTM2016-77278-P and MTM2013-40998-P, an AGAUR grant number 2014SGR-568, the grant FP7- PEOPLE-2012-IRSES 318999, and a CAPES grant number 88881.030454/2013-01 from the program CSF-PVE. The third author is supported by CAPES grant num- ber 99999.006888/2015-01 from the program CAPES-PDSE. References [1] G. Blé, V. Castellanos, J. Llibre, I. Quilantán, Integrability and global dynamics of the May-Leonard model, Nonlinear Anal. 14 (2013), 280–293. [2] A. Cima and J. Llibre, Bounded polynomial vector fields, Trans. Amer. Math. Soc. 318 (1990), 557–579. [3] P. Hartman, A lemma in the theory of structural stability of differential equations, Proc. Amer. Math. Soc. 11 610–620. [4] F. Dumortier, J. Llibre and J.C. Artés, Qualitative Theory of Planar Differential Sys- tems, Universitext (Springer, New York, 2006). [5] J. Llibre and R.D.S. Oliveira, Quadratic systems with invariant straight lines of total multiplicity two having Darboux invariants, Communications in Contemporary Mathematics 17 (2015), 145001, pp 17. [6] J. Llibre and C. Valls, Polynomial, rational and analytic first integrals for a family of 3-dimensional Lotka-Volterra systems, Z. Angew. Math. Phys. 62 (2011), 761–777. [7] R.M. May and W.J. Leonard, Nonlinear aspects of competition between three species, SIAM J. Appl. Math. 29 (1975) 243–253. 1 IBILCE–UNESP, CEP 15054–000, S. J. Rio Preto, São Paulo, Brazil E-mail address: buzzi@ibilce.unesp.br, robson.trevizan@yahoo.com.br MAY-LEONARD LOTKA-VOLTERRA SYSTEMS 11 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bel- laterra, Barcelona, Catalonia, Spain E-mail address: jllibre@mat.uab.cat