Final evolutions of a class of May-Leonard Lotka-Volterra systems

We study a particular class of Lotka-Volterra 3-dimensional systems 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 non-negative octant of ℝ3 including its infinity. This can be done because this differential system possesses a Darboux invariant.


Introduction
Polynomial ordinary differential systems are often used in various branches of applied mathematics, physics, chemist, engineering, etc. Models studying the interaction 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 [8] 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 differential system in R 3 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 [10] and [2] and on Lotka-Volterra systems [9], and the references quoted there.
The Lotka-Volterra systems in R 3 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.1) is reduced tȯ x = x− (1+ a+ b)x 2 , because the other equations do not provide any further information. Therefore, the bisectrix of the non-negative octant is an invariant straight line for this differential system.
In this paper we describe the global dynamics of system (1.1) in function of the parameters a and b when a + b = −1. The system (1.1) is defined in R 3 . In order to study the dynamics of its orbits at infinity we extend analytically its flow by using the Poincaré compactification of R 3 . In the appendix we give precise definitions for this compactification. The region of interest in our study is the non-negative octant of R 3 , i.e. where x ≥ 0, y ≥ 0, z ≥ 0. So we shall study the flow of the Poincaré compactification in the region We remark that the global dynamics of the May-Leonard system (1.1) with a + b = −1 can be studied because this differential system has a Darboux invariant.
The differential system (1.1) has been extensively studied in order to understand the interaction between species and try to predict possible extinction or overpopulation for example. However our interest is purely mathematical, we want to illustrate how the Darboux invariant can be used to describe the global dynamics of a differential system. Note that we are interested in the study of system (1.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 R 3 .

Statement of the main results
We denote by X the polynomial vector field associated to the differential system (1.1), and by p(X ) the Poincaré compactification of X , see the appendix in section 5. The flow of system (1.1) in the region R is described in the next two theorems. For a formal definition of topologically equivalent phase portraits see [5].
Theorem 2.1. For the May-Leonard differential system (1.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.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.1) on R ∞ = ∂ R ∩ {x 2 + y 2 + z 2 = 1} (i.e. the phase portrait at the infinity of the non-negative octant of R 3 ) is topologically equivalent to the one described in Fig. 1 Fig. 2 Let p(γ) denote the orbit γ of the vector field X associated to system (1.1) in the Poincaré compactification p(X ).  (i) The α-limit set of p(γ) is either the origin of R 3 , or the heteroclinic loop connecting the singular points p 1 , p 2 and p 3 , or the heteroclinic loop connecting the singular points p ∞ x , p ∞ y and p ∞ z (see Fig. 1(a)). (ii) The ω-limit set of p(γ) is either p ∞ b (the positive endpoint of the bisectrix x = y = z), or the heteroclinic loop connecting the singular points p ∞ x , p ∞ y and p ∞ z (see Figure 1(b)).
An immediate consequence of Theorem 2.2 is the following result.
Corollary 2.1. All orbit γ of system (1.1) with a+ b = −1 such that p(γ) is contained in the interior of R has their α-limit in xyz = 0 and their ω-limit in R ∞ . PSfrag replacements

Proof of Theorem 1
The following two lemmas will be useful to the proof of Theorem 1. PSfrag replacements Proof. The finite singular points of system (1.1) with a + b = −1 are the solutions of the system Since A > 0 for a ∈ R and the region of interest is R, we have: (i) If a ≤ −2 or a ≥ 1 system (1.1) has only four finite equilibrium points: p 0 , p 1 , p 2 and p 3 .
All these finite equilibrium points are hyperbolic if a = −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 [4].
We note that when a ∈ (−2, 1) and a → 1 we have that p 4 → p 3 , p 5 → p 1 and p 6 → p 2 ; while if a → −2 we have that p 4 → p 2 , p 5 → p 3 and p 6 → p 1 . This behavior of these equilibria allows to determine by continuity the local phase portraits on the boundary of R of the non-hyperbolic equilibrium points p 1 , p 2 and p 3 when a = −2 and a = 1 from the global phase portraits of the boundary of R when a ∈ (−2, 1).
The linear part of system (1.1) at the equilibrium p 0 is the identity matrix. Therefore it is a repelling equilibrium.
The eigenvalues of linear part at equilibrium points p 1 , p 2 and p 3 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 p 4 , p 5 and p 6 are 3, −1 and (−2 + a + a 2 )/A. Since −2 + a + a 2 < 0 for −2 < a < 1 then these equilibria have a 2-dimensional stable manifold and an 1-dimensional unstable one. Moreover, p 4 (respectively p 5 and p 6 ) is an attractor restricted to the invariant boundary x = 0 (respectively y = 0 and z = 0). Now we explain in few words, what we mean by saying that the local phase portraits, for a = −2 and a = 1, can be determined by continuity. For example, if a ∈ (−2, 1) then, on the plane x = 0, we have that p 4 is a node and p 2 is a saddle. If a tends to −2 then p 4 tends to p 2 and we have a saddle-node bifurcation. For a = −2 we have that p 4 = p 2 is a saddle-node singularity with two hyperbolic sectors in the half-plane {x = 0, y < 0} and a parabolic sector in the half-plane {x = 0, y > 0}. So, in a neighborhood of p 2 contained in the half-plane {x = 0, y > 0}, the local phase portraits are the same for a = −2 or a < −2 in the non-negative octant. The same analysis can be done for the other points p 1 and p 3 , and for the case when a tends to 1. Proof. 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.1). See appendix for details. Thus the differential system (1.1) in the local chart U 1 becomeṡ (3.1) So system (1.1), with a + b = −1 and satisfying a ≤ −2 or a ≥ 1, has two equilibrium points at infinity: (0, 0, 0) and (1, 1, 0). We call the first one p ∞ x the positive endpoint of the x-axis, and the second one p ∞ b the positive endpoint of the bisectrix x = y = z. The linear part at the equilibrium p ∞ x has the eigenvalues 1 − a and 2 + a at infinity and eigenvalue 1 in its finite direction. Therefore, on the infinity p ∞ x is a saddle such that its stable separatrix is contained in the z 1 -axis when a ≤ −2 (respectively z 2 -axis when a ≥ 1).
The eigenvalues of linear part at equilibrium p ∞ b are (−3 ± i √ 3(1 + 2a))/2 and 0. Therefore, on the infinity p ∞ b is a stable focus turning clockwise if a ≤ −2 (respectively counterclockwise if a ≥ 1). Now system (1.1) in the local chart U 2 writeṡ

(3.2)
Since the local chart U 2 covers the end part of the plane x = 0 at infinity of the non-negative octant of R 3 , we are interested only in the equilibrium points which are on z 1 = 0 and z 3 = 0. In this case, with a + b = −1 and satisfying a ≤ −2 or a ≥ 1, there is one equilibrium point at infinity: (0, 0, 0). We call this equilibrium p ∞ y the positive endpoint of the y-axis. The eigenvalues of linear part at equilibrium p ∞ y are 1 − a and 2 + a at infinity and eigenvalue 1 in its finite direction. Therefore, on the infinity p ∞ y is a saddle such that its stable separatrix is contained in the z 2 -axis when a ≤ −2 (respectively z 1 -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 U 3 . We call this equilibrium point p ∞ z . In the local chart U 3 system (1.1) becomeṡ The linear part at the equilibrium point p ∞ z has eigenvalues 1 − a and 2 + a at infinity and eigenvalue 1 in its finite direction. Therefore the origin of the local chart U 3 is a saddle such that its stable separatrix is contained in the z 1 -axis when a ≤ −2 (respectively z 2 -axis when a ≥ 1).
We have proved that the phase portrait of system (1.1) on R ∞ , for a + b = −1 and a ≤ −2, is the one presented in Figure 1(b). In the same figure, reversing the sense of all the orbits, we have the global dynamics on R ∞ , for a + b = −1 and a ≥ 1.
It remains to study the infinite equilibrium points of system (1.1) in case −2 < a < 1. In the local chart U 1 (1, 1, 0). The eigenvalues of linear part at equilibrium p ∞ x 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 p ∞ xy (respectively p ∞ xz ) has the eigenvalues −(2 + a), A/(1 − a) and 3A/(1 − a) (respectively (a − 1), A/(2 + a) and 3A/(2 + a)). Therefore the equilibria p ∞ xy and p ∞ xz are saddles such that their stable separatrix is contained in the z 1 -axis and z 2 -axis respectively. The eigenvalues of the linear part at the equilibrium p ∞ b are (−3 ± i √ 3(1 + 2a))/2 and 0. Therefore, on the infinity p ∞ b 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 p ∞ b is stable node. Now since we are interested only in the equilibrium points which are on z 1 = 0 and z 3 = 0, in the local chart U 2 the system (3.2) has two equilibrium points: p ∞ y = (0, 0, 0) and p ∞ yz = (0, (2 + a) /(1 − a), 0). The eigenvalues of linear part at equilibrium p ∞ y 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 (0, (2 + a)/ (1 − a), 0) has the eigenvalues −(2 + a), A/(1 − a) and 3A/(1 − a) at infinity. Therefore, on the infinity p ∞ yz = (0, (2 + a)/ (1 − a), 0) is a saddle such that its stable separatrix is contained in the z 2 -axis.
In the local chart U 3 we only need to study the equilibrium point p ∞ z = (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.
We also observe that when a ∈ (−2, 1) and a → 1 we have that 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).
We will show now that does not exist a periodic orbit on R ∞ . For this we will need Bautin's Theorem, which is proved in [1].
Theorem 3.1 (Bautin's Theorem). A quadratic polynomial differential system of the forṁ has no limit cycles. Proof. The differential system (1.1) in the local chart U 1 is given by (3.1). Making z 3 = 0 to determine the dynamics on R ∞ we have the systeṁ

Proof of Theorem 2
We say that a C 1 function I(x, y, z,t) is an invariant of the polynomial differential system (1.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.1) is defined. When the function I is independent of the time, then it is called a first integral of differential system (1.1). Also if an invariant I(x, y, z,t) is of the form f (x, y, z)e st , then it is called a Darboux invariant. Proof. It is immediate to check that whereẋ,ẏ andż are given in (1.1). Therefore I is a Darboux invariant of system (1.1).
For knowing how to obtain the Darboux invariant given in Proposition 4.1 see statement (vi) of Theorem 8.7 of [5], there the theory is described for polynomial vector fields in R 2 , but the results and the proofs extend to R 3 . Here α(p) and ω(p) denote the α-limit and ω-limit sets of p respectively, and S 2 denotes the boundary of the Poincaré ball, i.e. the infinity of R 3 .
For a proof of Proposition 4.2 see [7]. Proof. Let q ∈ ω(p). Then there exists a sequence (t n ) with t n → +∞, such that ϕ p (t n ) = (x(t n ), y(t n ), z(t n )) → q when n → +∞. By Proposition 4.2 we know that q ∈ {(x, y, z) ∈ R : xyz = 0} ∪ R ∞ . Suppose that q / ∈ R ∞ and take ε > 0. Then there exist M > 0 and n 0 ∈ N such that xyz ≤ M for all x, y, z ∈ B(q, ε) and ϕ p (t n ) ∈ B(q, ε) for all n ≥ n 0 . Therefore x(t n )y(t n )z(t n ) ≤ M for all n ≥ n 0 . On the other hand, as lim t→+∞ x(t)y(t)z(t) = +∞ then exists t 0 ∈ R such that x(t)y(t)z(t) ≥ M for all t ≥ t 0 , which is a contradiction. Therefore q ∈ R ∞ .
Proof. [Proof of Theorem 2.2] Let p(γ) = {ϕ p (t) = (x(t), y(t), z(t)) : t ∈ R} be the orbit of the Poincaré compactification of system (1.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 4.1 and 4.2 the αand ω-limit set of p(γ) is contained in boundary of R, i.e. in {(x, y, z) ∈ R : xyz = 0} ∪ R ∞ . Furthermore, by Proposition 4.1, I(t, x(t), y(t), z(t)) = k constant with k > 0. So This implies that the orbit p(γ) tends to the set xyz = 0 when t → −∞. Looking at the dynamics of the flow of the compactified vector field on xyz = 0 described in Fig. 1(a) if a ≤ −2 or a ≥ 1, and in Fig. 2(a) if −2 < a < 1, we consider the following two cases.