Electronic Journal of Differential Equations, Vol. 2020 (2020), No. 55, pp. 1–19.

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu

GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM WITH

A DARBOUX INVARIANT

REGILENE OLIVEIRA, CLAUDIA VALLS

Abstract. We study the global dynamics of the classic May-Leonard model

in R3. Such model depends on two real parameters and its global dynamics

is known when the system is completely integrable. Using the Poincaré com-
pactification on R3 we obtain the global dynamics of the classical May-Leonard

differential system in R3 when β = −1 − α. In this case, the system is non-

integrable and it admits a Darboux invariant. We provide the global phase
portrait in each octant and in the Poincaré ball, that is, the compactification

of R3 in the sphere S2 at infinity. We also describe the ω-limit and α-limit of
each of the orbits. For some values of the parameter α we find a separatrix

cycle F formed by orbits connecting the finite singular points on the boundary

of the first octant and every orbit on this octant has F as the ω-limit. The
same holds for the sixth and eighth octants.

1. Introduction and statement of main results

The Lotka-Volterra systems in R3

ẋi = xi

( 3∑
j=1

aijxj + bi

)
(i = 1, . . . , 3).

were introduce by Lotka and Volterra in the 1920’s for describing the interaction
among species, see [8, 12, 17]. Since their appearance, these systems and related
ones have been intensively investigated because they are also important for de-
scribing different phenomena in biology, ecology and chemistry; see for instance
[3, 5, 6, 7, 9, 11, 13, 15] and the references therein.

Because of its simplicity, special attention has attracted the so-called classical
May-Leonard system introduced by May and Leonard in [14] which again describes
a competition between three species but is a simple model with a rich dynamical
behavior. This system can be written as

ẋ = x(1− x− αy − βz),
ẏ = y(1− βx− y − αz),
ż = z(1− αx− βy − z),

(1.1)

where α and β are real parameters.

2010 Mathematics Subject Classification. 37C15, 37C10.
Key words and phrases. Lotka-Volterra systems; May-Leonard systems;

Darboux invariant; phase portraits; limit sets; Poincaré compactification.
c©2020 Texas State University.

Submitted March 1, 2019. Published June 3, 2020.

1



2 R. OLIVEIRA, C. VALLS EJDE-2020/55

It was showed in [14] that, whenever α+ β 6= −1, system (1.1) has four singular
points in R3

+ = {(x, y, z) ∈ R3, x, y, z ≥ 0}. Three of them are on the boundary
of R3

+,

E1 = (1, 0, 0), E2 = (0, 1, 0), E3 = (0, 0, 1),

and the fourth one is the interior point

C =
(
(1 + α+ β)−1, (1 + α+ β)−1, (1 + α+ β)−1

)
.

If in addition, α+β > 2 and either α < 1, or β > 1, then there is a separatrix cycle
F formed by orbits connecting E1, E2 and E3 on the boundary of R3

+, and every
orbit in R3

+, except of the equilibrium point C has F as the ω-limit. Further, it was
shown in [2] that in the degenerate case α+ β = 2, the cycle F becomes a triangle
on the invariant plane x + y + z = 1 and all orbits inside the triangle are closed
and every orbit in the interior of R3

+ has one of these closed orbits as its ω-limit.
Moreover, in [2] it was studied completely the dynamics of the May-Leonard system
whenever α + β = 2, or α = β. Other dynamical aspects or other values of the
parameters α, β were also considered in [1, 16] and references therein.

In this paper we shall restrict to the cases where either α+ β 6= 2, or α 6= β and
study the May-Leonard systems (1.1) that admit a Darboux invariant of the form
estf(x, y, z), where f(x, y, z) is given by the product of invariant planes (see more
details below). The existence of such a Darboux invariant provides information
about the ω- and α-limit of all orbits of the system. The case where the May-
Leonard system (1.1) admits a Darboux invariant as above is whenever

α+ β + 1 = 0.

Here we consider this case and we describe the global dynamics in the compactifi-
cation of R3 in function of α. In particular, we provide the ω- and α-limits of all
orbits of the system in the positive octant, and so determining the initial and final
evolution of the three species considered in the May-Leonard model. We recall that
the global description of the flow of a differential system in R3 is generally very
difficult. In that paper, using the Poincaré compactification and the fact that we
have an invariant we are able to do it. The case α + β = −1 is usually forgotten
in the literature because in this case the point C does not exist and the dynamics
become more complex. As far as the authors know this is the first paper in which
the May Leonard problem under the above assumptions is completely studied from
the dynamical point of view. We recall that loosely speaking the Poincaré ball is
obtained identifying R3 with the interior of the 3-dimensional ball of radius one
centered at the origin, and extending analytically the flow of system (1.1) to the
boundary of S2 on that ball (see the Appendix for details).

2. Statements of main results

An invariant of system (1.1) on an open subset U of R3 is a nonconstant C1

function I in the variables x, y, z and t such that I(x(t), y(t), z(t), t) is constant on
all solution curves (x(t), y(t), z(t)) of system (1.1) contained in U , i.e.

x(1−x−αy−βz)∂I
∂x

+y(1−βx−y−αz)∂I
∂y

+z(1−αx−βy−z)∂I
∂z

+
∂I

∂t
= 0, (2.1)

for all (x, y, z) ∈ U .



EJDE-2020/55 GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM 3

On the other hand given f ∈ C[x, y, z] we say that the surface f(x, y, z) = 0 is an
invariant algebraic surface of system (1.1) if there exists K ∈ C[x, y, z] such that

x(1−x−αy−βz)∂f
∂x

+y(1−βx−y−αz)∂f
∂y

+z(1−αx−βy−z)∂f
∂z

= Kf. (2.2)

The polynomial K is called the cofactor of the invariant algebraic surface f = 0.
For more details see [4, Chapter 8].

An invariant I is called a Darboux invariant if it can be written in the form

I(x, y, z, t) = fλ1
1 · · · fλp

p es t,

where, for i = 1, . . . p, fi = 0 are invariant algebraic surfaces of system (1.1), λi ∈ C,
and s ∈ R \ {0}.

Theorem 2.1. The following statements hold for system (1.1) with either α+β 6= 2
or α 6= β.

(a) It has a Darboux invariant of the form I(x, y, z, t) = estf(x, y, z), where f
is given by the product of invariant planes if and only if α+ β + 1 = 0.

(b) It is invariant under the symmetry (x, y, z)→ (y, z, x).
(c) The ω-limit of any of its orbits in R3 is contained in Ω union with its

boundary at infinity in the Poincaré compactification in R3, where

Ω = {(x, y, z) ∈ R3 : x = 0} ∪ {(x, y, z) ∈ R3 : y = 0} ∪ {(x, y, z) ∈ R3 : z = 0}.

From now on, we consider system (1.1) restricted to β = −1 − α with α ∈ R,
that is

ẋ = x(1− x− αy + (1 + α)z),

ẏ = y(1 + (1 + α)x− y − αz),
ż = z(1− αx+ (1 + α)y − z).

(2.3)

To describe de global dynamics of system (2.3) we define the ith-octants, Oi, for
i = {1, 2, . . . , 8}, as

O1 = {(x, y, z) ∈ R3 : x ≥ 0, y ≥ 0, z ≥ 0},
O2 = {(x, y, z) ∈ R3 : x ≤ 0, y ≥ 0, z ≥ 0},
O3 = {(x, y, z) ∈ R3 : x ≤ 0, y ≤ 0, z ≥ 0},
O4 = {(x, y, z) ∈ R3 : x ≥ 0, y ≤ 0, z ≥ 0},
O5 = {(x, y, z) ∈ R3 : x ≥ 0, y ≥ 0, z ≤ 0},
O6 = {(x, y, z) ∈ R3 : x ≤ 0, y ≥ 0, z ≤ 0},
O7 = {(x, y, z) ∈ R3 : x ≤ 0, y ≤ 0, z ≤ 0},
O8 = {(x, y, z) ∈ R3 : x ≥ 0, y ≤ 0, z ≤ 0}.

The global dynamics of system (2.3) in O1 (the positive octant) is given in the
following theorem. In that theorem, we denote by O+

1 the interior of O1, i.e.,

O+
1 = {(x, y, z) : x > 0, y > 0, z > 0}.

Theorem 2.2. The following statements hold for system (2.3) restricted to O1 for
α ∈ (−∞, 1):

(a) The phase portraits in the Poincaré disc of system (2.3) restricted to the
invariant planes are topologically equivalent to one of the phase portraits of
Figure 1.



4 R. OLIVEIRA, C. VALLS EJDE-2020/55

(b) The phase portraits of system (2.3) at the infinity of O1 are topologically
equivalent to one of the phase portraits given in Figure 2. More precisely,

(b.1) for α ≤ −2 the boundary of the infinity of O1 is a heteroclinic cycle
formed by three equilibrium points coming from the ones located at
the end of the three positive half-axes of coordinates, and three orbits
connecting these equilibria, each one coming from the orbit at the end
of every plane of coordinates. In the interior of the infinity of O1 there
is an attractor whose orbits fill completely this interior.

(b.2) For α ∈ (−2, 1) the boundary of the infinity of O1 is a graph formed
by six equilibrium points all of them located at the positive half-axes of
coordinates. Three of them are at the end of the axes and the other
three are between them. Moreover, there are six orbits connecting these
equilibria. Each of them coming from the orbit at the end of every plane
of coordinates. In the interior of the infinity of O1 there is an attractor
of all the orbits coming from the equilibria at the end of every plane of
coordinates.

(c) The phase portrait of system (2.3) on O1 is giving in Figure 3(i) when
α ≤ 2 and Figure 3(ii) when α ∈ (−2, 1). Namely,

(c.1) for α ≤ −2, there exists a separatrix cycle F formed by orbits connect-
ing the finite singular points on the boundary of O1 and every orbit on
O1, except the origin, has F as its ω-limit. The α-limit set of the orbits
on O+

1 \F is formed by four equilibrium points, the origin and the three
equilibria located at the end of the positive half-axes of coordinates.

(c.2) for α ∈ (−2, 1), there exists a graph G formed by orbits connecting the
finite singular points on the boundary of O1, except the origin. The ω-
limit set of every orbit on O1, except the origin, is one of the vertices
of G. The α-limit set of the orbits on O+

1 \ G is formed by seven
equilibrium points, the origin and the equilibria located at the end of
the invariant planes.

Since Theorem 2.2 provides the ω-limit and the α-limit of all orbits inside O1

(which is the octant where system (2.3) has biological meaning), in that theorem
we are determining all the initial and final evolution of the three species considered
by system (2.3) according to the values of the parameter α. We recall that the
statements in Theorem 2.2 also hold for α ≥ 1. Indeed, the dynamics for α > 1 is
the same as the one for α < −2 reversing the time and the dynamics for α = 1 is
the same of the one for α = −2 also reversing the time. It follows from Theorem 2.1
(b) that the dynamics in O6 and O8 are the same as the one in O1.

Theorem 2.3. The phase portrait of system (2.3) on O2 is given in Figure 4(i)
when α < −2, or α > 1; in Figure 4 (ii) when α = −2, or α = 1; and in
Figure 4(iii) when α ∈ (−2, 1). More precisely, all orbits in O2 are heteroclinic and

(a) for α < −2, all orbits starting on O2 \A, where A = {z = 0}∪{x = y = 0},
go in forward time to infinity to an equilibrium point P0 located outside the
end of the invariant planes (see Figure 4 (i)). The same holds for α > 1
in O2 \B, where B = {x = 0} ∪ {y = z = 0}.

(b) for α ∈ [−2, 1], all orbits starting on O2 \ {y = 0} go in forward time to
infinity to the endpoint P0 of the negative half-axes of coordinates (y-axes).
See Figures 4 (ii) and (iii).



EJDE-2020/55 GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM 5

Figure 1. Phase portrait of system (2.3) on the Poincaré disc
restricted to the invariant planes x = 0, y = 0 or z = 0 when: (i)
α < −2, or α > 1; (ii) α = −2, or α = 1; (iii) −2 < α < 1.

Figure 2. Phase portrait of system (2.3) on the Poincaré sphere
when: (i) α < −2 or α > 1; (ii) α = −2 or α = 1; (iii) −2 < α < 1.

Theorem 2.4. The phase portrait of system (2.3) on O3 is given in Figure 5(i)
when α < −2, or α > 1; in Figure 5(ii) when α = −2, or α = 1; and in Figure 5(iii)
when α ∈ (−2, 1). Namely, all orbits in O3 are heteroclinic and

(a) for α < −2, all orbits starting on O3\{y = 0} go in forward time to infinity
to an equilibrium point Q0 located outside the end of the invariant planes
(see Figure 5(i)). The same is true for α > 1 in O3 \ {x = 0}.

(b) for α = −2, all orbits starting on O3\{y = 0} go in forward time to infinity
to the endpoint Q0 of the negative half-axes of coordinates (x-axes). For
α = 1, all orbits starting on O3 \ {x = 0} go in forward time to infinity
to the endpoint Q0 of the negative half-axes of coordinates (y-axes). See
Figure 5(ii).

(c) for α ∈ (−2, 1), all orbits starting on O3 \ {x = 0} go in forward time to
infinity to an equilibrium point Q0 located outside the end of the invariant
planes (see Figure 5(iii)).

Theorem 2.5. The phase portrait of system (2.3) on O5 is given in Figure 6(i)
when α < −2, or α > 1; in Figure 6(ii) when α = −2, or α = 1; and in Figure 6(iii)
when α ∈ (−2, 1). Namely, all orbits in O5 are heteroclinic and

(a) for α < −2, all orbits starting on O5\C, where C = {x = 0}∪{y = z = 0},
go in forward time to infinity to an equilibrium point R0 located outside the



6 R. OLIVEIRA, C. VALLS EJDE-2020/55

Figure 3. Phase portrait of system (2.3) in the first octant O1:
(i) for α ≤ −2 or α ≥ 1 and (ii) for −2 < α < 1. On the left-hand
side we draw the phase portrait of system (2.3) restricted to the
positive invariant planes and on the right-hand side we draw the
phase portrait of system (2.3) on the Poincaré sphere in O1.

end of the invariant planes (see Figure 6(i)). The same holds for α > 1 in
O5 \D, where D = {y = 0} ∪ {x = z = 0}.

(b) for α ∈ [−2, 1], all orbits starting on O5 \ {z = 0} go in forward time to
infinity to the endpoint R0 of the positive half-axes of coordinates (z-axes).
See Figures 6 (ii) and (iii).

Theorems 2.2–2.5 are proved in section 6. We remark that for Figures 4-6 not
be too crowed, which makes difficult to comprehend, we draw only the separatrices
of each equilibrium. The structure of this article is as follows. In sections 3 and
4 we describe the dynamics on the Poincaré sphere and on the invariant planes,
respectively. In section 5 we provide the results that will be used in the proof of
Theorem 2.1 regarding the α- and ω-limits. We have also included an Appendix
with the Poincaré compactification in both R2 and R3.

3. Dynamics on the Poincaré sphere S2

In this section we present the analysis of the flow of system (2.3) at infinity
using the Poincaré compactification of the system in R3, described in the local
charts Ui, Vi for i = 1, 2, 3.



EJDE-2020/55 GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM 7

Figure 4. Phase portrait of system (2.3) in the second octant
O2: (i) for α < −2 or α > 1, (ii) for α = −2 or α = 1 and
(iii) for −2 < α < 1. On the left-hand side we draw the phase
portrait of system (2.3) restricted to the invariant planes and on
the right-hand side we draw the phase portrait of system (2.3) on
the Poincaré sphere in O2.



8 R. OLIVEIRA, C. VALLS EJDE-2020/55

Figure 5. Phase portrait of system (2.3) in the third octant O3:
(i) for α < −2 or α > 1, (ii) for α = −2 or α = 1 and (iii) for −2 <
α < 1. On the left-hand side we draw the phase portrait of system
(2.3) restricted to the invariant planes and on the right-hand side
we draw the phase portrait of system (2.3) on the Poincaré sphere
in O3.

From Appendix 7.2 the expression of the Poincaré compactification p(X) of
system (2.3) in the local chart U1 is given by

ż1 = z1(2− z1 − z2 + α+ z1α− 2z2α),

ż2 = −z2(−1− z1 + 2z2 + α− 2z1α+ z2α),

ż3 = −z3(−1 + z2 + z3 − z1α+ z2α).

(3.1)

For z3 = 0 (which correspond to the points on the sphere S2 at infinity) system
(3.1) becomes

ż1 = z1(2− z1 − z2 + α+ z1α− 2z2α),

ż2 = −z2(−1− z1 + 2z2 + α− 2z1α+ z2α).
(3.2)

System (3.2) has the foloowing equilibrium points

p1 = (0, 0), p2 = (1, 1), p3 =
(

0 ,
1− α
2 + α

)
, p4 =

(−2− α
−1 + α

, 0
)
.

Note that p3 exists whenever α 6= −2 and p4, whenever α 6= 1. The eigenvalues of
the Jacobian matrix evaluated in each of the equilibria are 1− α and 2 + α for p1,
(−3±

√
3i|(2α−1)|)/2 for p2, −1+α and 3(1 + α+ α2)/(2 + α) for p3, and −(2+α)

and −3(1+α+α2)/(−1+α) for p4. So, p2 is a stable focus and, except when α = 1



EJDE-2020/55 GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM 9

Figure 6. Phase portrait of system (2.3) in the fifth octant O5:
(i) for α < −2 or α > 1, (ii) for α = −2 or α = 1 and (iii) for −2 <
α < 1. On the left-hand side we draw the phase portrait of system
(2.3) restricted to the invariant planes and on the right-hand side
we draw the phase portrait of system (2.3) on the Poincaré sphere
in O5.

or α = −2, these points are hyperbolic, whose topological type depend on the sign
of α− 1 and α+ 2.

The flow in the local chart V1 is the same as the flow in the local chart U1 because
the compacted vector field p(X) in V1 coincides with the vector field p(X) in U1

multiplied by −1. Hence the phase portrait on the chart V1 is the same as the one
in U1 reserving in an appropriate way the direction of the time.

The expression of the Poincaré compactification p(X) of system (2.3) in the local
chart U2 is

ż1 = −z1(−1 + 2z1 − z2 + α+ z1α− 2z2α)

ż2 = z2(2− z1 − z2 + α− 2z1α+ z2α),

ż3 = −z3(−1 + z1 + z3 + z1α− z2α).

(3.3)

System (3.3) restricted to z3 = 0 becomes

ż1 = −z1(−1 + 2z1 − z2 + α+ z1α− 2z2α)

ż2 = z2(2− z1 − z2 + α− 2z1α+ z2α),
(3.4)



10 R. OLIVEIRA, C. VALLS EJDE-2020/55

which has the following equilibria

q1 = (0, 0), q2 = (1, 1), q3 =
(

0,−2 + α

α− 1

)
, q4 =

(1− α
2 + α

, 0
)
.

Note that the equilibria q2 = (1, 1) and q4 = (−(−1 + α)/(2 + α), 0) were already
studied in the local chart U1. The eigenvalues of the equilibria q1 are 1 − α and
2+α and of the equilibria q3 are −(2+α) and −3(1+α+α2)/(2+α) for q3. Again
the topological type of these two singular points depend on the sign of 2 + α and
1− α and q3 does not exist when α = 1.

The flow in the local V2 is the same as the flow in the local chart U2 reserving
in an appropriate way the direction of the time.

The expression of the Poincaré compactification p(X) of system (2.3) in the local
chart U3 is exactly the same as in in the local chart U1, because of the symmetry.
Consequently, the flow at infinity in the local chart V3 is the same as the flow in
the local chart U3 and U1 reversing appropriately the time.

From the study presented in section 3 the topological type of each equilibrium
point in the local chart Ui, i = 1, 2, 3 depends on the sign of α − 1 and α + 2 so
we split the study of the dynamics on the Poincaré sphere in five cases: α < −2,
α > 1, α ∈ (−2, 1), α = −2 and α = 1. In Table 1 we provide the description of the
topological type of each equilibria (pi and qj , i = 1, 2, 3, 4 and j = 1, 3) according
to the value of the parameter α in each the local charts U1 and U2, respectively (we
recall that the dynamics in U3 coincide with the one in U1).

Table 1. Topological type of each equilibria in the local chart U1

and U2 according to the value of α, where pi ∈ U1, i = 1, 2, 3, 4
and qj ∈ U2, j = 1, 3.

α < −2 α = −2 α ∈ (−2, 1) α = 1 α > 1

p1 saddle saddle-node node saddle-node saddle

p2 unstable focus unstable focus
unstable focus or node

when α = −1/2.
unstable focus unstable focus

p3 stable node – saddle – unstable node

p4 unstable node – saddle – stable node

q1 saddle saddle-node node saddle-node saddle

q3 unstable node – saddle – stable node

From the local behavior of the orbits of system (2.3) in each local chart of the
Poincaré sphere together with the existence of the invariant planes x = 0, y = 0
and z = 0 and because the boundary of the invariant planes at infinity of R3 are
invariant, we get five global phase portraits in the Poincaré sphere, according to the
cases studied above. Note that the global phase portraits for α = −2 and α = 1 are
topologically equivalent up to a rescaling of the time, as well as the global phase
portraits for α < −2 and α > 1. In Figure 2 we only draw the non topologically
equivalent phase portraits in the Poincaré sphere for system (2.3).



EJDE-2020/55 GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM 11

4. Dynamics of system (2.3) restricted to the invariant planes

Now we restrict our study to each coordinate plane. On the invariant plane
x = 0 system (2.3) becomes

ẏ =y(1− y − αz),
ż =z(1 + (1 + α)y − z).

(4.1)

This system has the following equilibria

r0 = (0, 0), r1 = (1, 0), r2 = (0, 1), r3 =
( 1− α

1 + α+ α2
,

α+ 2

1 + α+ α2

)
.

The eigenvalues of the Jacobian matrix evaluated at each of the equilibria are 1, 1
for r0; −1, 2+α for r1; −1, 1−α for r2; and −1, (2 + α)(α− 1)/(1 + α+ α2) for r3.
So r0 is a unstable node and the other points are hyperbolic whose topological type
depend on the sign of α− 1 and α+ 2. Note that r3 coalesces with one of the other
equilibria when either α = −2 or α = 1. The topological type of each of these
equilibria is described in Table 2 according to the values of α.

Table 2. Topological type of each equilibria in the invariant plane
x = 0 according to the value of α.

α < −2 α = −2 α ∈ (−2, 1) α = 1 α > 1

r0 unstable node unstable node unstable node unstable node unstable node

r1 stable node saddle–node saddle saddle saddle

r2 saddle saddle saddle saddle–node stable node

r3 saddle – stable node – saddle

From Appendix 7.1 the expression of the Poincaré compactification of system
(4.1) in the local chart U1 is

ż1 =2z1 − z21 + z1α+ z21α,

ż2 =z2 − z22 + z1z2α.

There are two equilibria on z2 = 0, namely s0 = (0, 0) and s1 = ((α+2)/(1−α), 0).
Note that s1 exists whenever α 6= 1 and coalesces with s0 when α = −2. The
eigenvalues of the Jacobian matrix evaluated at these two equilibria are 1, (2 + α)
for s0 and −2 − α, (1 + α + α2)/(1 − α) for s1. Again they are hyperbolic whose
topological type depend on the sign of α− 1 and α+ 2.

The expression of the Poincaré compactification of system (4.1) in the local
chart U2 is

ż1 =z1 − 2z21 − z1α− z21α,
ż2 =z2 − z1z2 − z22 − z1z2α.

So the point s2, the origin of the local chart U2, is an equilibrium point whose
eigenvalues of its Jacobian matrix are 1, 1−α. The topological types of each of the
equilibria in the local charts U1 and U2 are described in Table 3 according to the
values of α.

Combining the above analysis in the finite part and at each local chart at infinity
for each of the values of α we get five global phase portraits of system (4.1) in the
Poincaré disc. We remark that the global phase portraits when α > 1 and α = 1
are topologically equivalent to the cases α < −2 and α = −2, respectively. So,



12 R. OLIVEIRA, C. VALLS EJDE-2020/55

Table 3. Topological type of each equilibria of system (4.1) in the
local chart U1 and U2 according to the value of α.

α < −2 α = −2 α ∈ (−2, 1) α = 1 α > 1

s0 saddle saddle–node unstable node unstable node unstable node

s1 unstable node – saddle — unstable node

s2 unstable node unstable node unstable node saddle–node saddle

in Figure 1 we only draw the three non topologically equivalent phase portraits of
system (4.1).

On the invariant plane y = 0 system (2.3) becomes

ẋ =x(1− x+ (1 + α)z),

ż =z(1− z − αx).
(4.2)

Such system is equivalent to system (4.1) by the change of coordinates (x, z) →
(z, y), in other words, the dynamics of system (4.2) is equivalent to the one of
system (4.1) up to a rotation, see Figure 1.

On the invariant plane z = 0 system (2.3) becomes

ẋ =x(1− x− αy),

ẏ =y(1 + (1 + α)x− y).
(4.3)

Such system is equivalent to system (4.1) by the change of coordinates (x, y) →
(y, z). So the dynamics of system (4.3) is topologically equivalent to the dynamics
of system (4.1) up to a rotation and it is drawn in Figure 1.

5. On the ω-limit of the orbits of system (2.3) on O1

In this section we provide the ω-limit of all orbits of system (2.3) in O1. To do
it, we first introduce the notion of ω-limit and α-limit and we state and proof an
auxiliary result that will be used to prove the existence of the heteroclinic cycle
and of the graph in the positive octant O1.

Let φp(t) be the solution of system (1.1) passing through the point p ∈ R3,
defined on its maximal interval (αp, ωp) such that φp(0) = p. If ωp =∞, we define
the ω-limit set of p as

ω(p) = {q ∈ R3 : ∃{tn} with tn =∞ and φp(tn) = q when n =∞}.

In the same way, if αp = −∞, we define the α-limit set of p as

α(p) = {q ∈ R3 : ∃{tn} with tn = −∞ and φp(tn) = q when n =∞}.

For more details on the ω- and α-limit sets see for instance [4, section 1.4].
The existence of a Darboux invariant of system (1.1) provides information about

the ω- and α-limit sets of all orbits of system (1.1). More precisely, we have the
following result, where the definition of Poincaré compactification and Poincaré
sphere is given in subsection 7.2. Its proof can be found in [10] for the 2-dimensional
case so here we repeat it for the 3-dimensional case.

Proposition 5.1. Let S2 be the infinity of the Poincaré sphere and I(x, y, z, t) =
f(x, y, z)est be a Darboux invariant of system (1.1). Let also p ∈ R3 and φp(t)
be the solution of system (1.1) with maximal interval (αp, ωp) such that φp(0) =



EJDE-2020/55 GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM 13

p. If ωp = ∞ then ω(p) ⊂ {f(x, y, z) = 0} ∪ S2 and if αp = −∞ then α(p) ⊂
{f(x, y, z) = 0} ∪ S2.

Proof. Here we prove the first statement in the proposition since the second one
follows in the same lines. Assume that s > 0 and let φp(t) = (xp(t), yp(t), zp(t)).
Since I(x, y, z, t) is an invariant I(xp(t), yp(t), z(t)t) = a ∈ R for all t ∈ (αp, ωp), it
follows that

a = lim
t→∞

I(xp(t), yp(t), zp(t), t) = lim
t→∞

f(xp(t), yp(t), zp(t))e
st.

As limt→∞ est = ∞, we have that limt→∞ f(xp(t), yp(t), zr(t)) = 0. So, by conti-
nuity and the definition of ω-limit set it follows that ω(p) ⊂ {f(x, y, z) = 0}, and
for the α-limit set α(p) ∈ S2. �

Now we shall provide the ω-limit of all orbits of system (2.3) in O. For this, we
will distinguish between the cases α ≤ −2 and α ∈ (−2, 1). The results for α ≥ 1
are the same as the ones for α ≤ −2 since the dynamics is the same reversing the
time.

Proposition 5.2. Assume α ≤ 2. There exists a separatrix cycle F formed by
orbits connecting the finite equilibrium points on the boundary of O1 and every
orbit on O1 has F as the ω-limit.

Proof. Let V (x, y, z) = x+ y + z and denote by γ(t) = (x(t), y(t), z(t)) an orbit of
system (2.3). We define V (t) = V (γ(t)) and

V̇ =
dV

dt
(γ(t))

=
∂V

∂x
x(1− x− αy + (1 + α)z) +

∂V

∂y
y(1 + (1 + α)x− y − αz)

+
∂V

∂z
z(1− αx+ (1 + α)y − z).

Hence

V̇ = x+ y + z − (x2 + y2 + z2 + xy + xz + yz).

Applying the change of coordinates

x→ 1

3
(3w − 1), y → 1

3
(3w − 3v − 3u+ 2), z → 1

3
(3w − 6u+ 2),

the set {(x, y, z) ∈ R3 : V̇ = 0} is the elliptic paraboloid

E = {(u, v, w) ∈ R3 : 3w − 3u2 − v2 = 0}.

Clearly V̇ = 0 on the paraboloid E that contains the equilibrium points

(1, 0, 0), (0, 1, 0), (0, 0, 1).

On that equilibria the function V is equal to 1. Moreover,

V̇ = x+ y + z − (x2 + y2 + z2 + xy + xz + yz)

≤ x+ y + z − (x2 + y2 + z2 + 2xy + 2xz + 2yz)

= V − V 2 = V (1− V ).

So, V̇ < 0 on V > 1. Taking this into account, that the points of the interior of O1

satisfying V̇ = 0 are the ones on E and that V on the above equilibrium points is



14 R. OLIVEIRA, C. VALLS EJDE-2020/55

equal to one, we conclude that every orbit in the interior of O1 enters, and remains
inside, the set

S = {(x, y, z) ∈ O+
1 : 0 < V ≤ 1}.

Furthermore, in view of Theorem 2.1(c) all orbits in S approach Ω, so all orbits
on O1 approach S ∩ Ω and have their ω-limit in this set. In view of the Poincaré-
Bendixson theorem and that the equilibria (1, 0, 0), (0, 1, 0) and (0, 0, 1) can not be
the ω-limit of any orbit because they are saddles, the unique possibility is to be the
heteroclinic cycle F formed by orbits connecting the equilibrium points (1, 0, 0),
(0, 1, 0) and (0, 0, 1). This proves the existence of the separatrix cycle F as stated
in the proposition. �

Proposition 5.3. Assume α ∈ (−2, 1). There exists a graph G formed by orbits
connecting the finite equilibrium points on O1, except the origin. The ω-limit set
of all orbits on O1 is the set formed by the three equilibria that are not on the
half-positive axes.

Proof. Consider V as in the proof Proposition 5.2. Proceeding in the same manner
as in that proof, we have that the points of the interior of O1 satisfying V̇ = 0 are
the ones on the paraboloid E, that V on the equilibrium points (1, 0, 0), (0, 1, 0),
(0, 0, 1) is equal to 1 and that V on the equilibrium points

P1 =
1

1 + α+ α2
(0, α− 1, α+ 2), P2 =

1

1 + α+ α2
(α− 1, α+ 2, 0),

P3 =
1

1 + α+ α2
(α+ 2, 0, α− 1)

is equal to 3/(1 + α + α2) > 1. So, as in the proof of Proposition 5.2, every orbit
in the interior of O1 enters, and remains inside, the set

S1 = {(x, y, z) ∈ O+
1 : 0 < V ≤ 3/(1 + α+ α2)}.

Furthermore, in view of Theorem 2.1 (c) all orbits in S1 approach Ω, so all orbits
on O+

1 approach S1 ∩Ω and have their ω-limit in this set. In view of the Poincaré-
Bendixson theorem and that the equilibria (1, 0, 0), (0, 1, 0) and (0, 0, 1) can not
be the ω-limit of any orbits because they are saddles, the unique possibility is to
be either P1, or P2, or P3. This guarantees the existence of the graph G as in the
statement of the proposition. �

6. Proofs of Theorems 2.1–2.5

Proof of Theorem 2.1. Taking into account that α + β 6= 2 and α 6= β it follows
easily from (2.2) that the unique irreducible invariant planes of system (1.1) are
f1(x, y, z) = x = 0, f2(x, y, z) = y = 0 and f3(x, y, z) = z = 0. Moreover, the
f1(x, y, z) = 0, f2(x, y, z) = 0 and f3(x, y, z) = 0 have cofactors, 1 − x − αy − βz,
1 − βx − y − αz and 1 − αx − βy − z, respectively. Hence, system (1.1) admits
a Darboux invariant of the form I(x, y, z, t) = estf(x, y, z), with f = fr11 fr22 fr33 , if
and only if, equation (2.1) is satisfied for a real s 6= 0. Doing so we get

s = −3r2, r1 = r3 = r2, α+ β + 1 = 0.

Choosing r2 = 1 we have s = −3 and r1 = r2 = r3 = 1, being the Darboux invariant
I = e−3txyz. This concludes the proof of Theorem 2.1(a).

Statement (b) can be easily checked and statement (c) follows from Proposi-
tion 5.1 and the expression of the Darboux invariant obtained in statement (a). �



EJDE-2020/55 GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM 15

Proof of Theorem 2.2. The proof of statement (a) of Theorem 2.2 follows directly
from the study done in section 3.

To prove statement (b) we note that since the planes of coordinates and the
sphere at infinity are invariant, the intersection if formed by orbits.

If α ≤ −2 the open arc of the infinity of O1 in the local chart U1 corresponding
to the end of the plane z = 0 is formed by an orbit having as α-limit the equilibrium
(0, 0, 0). The open arc of the infinity of O1 corresponding to the end of the plane
y = 0 is formed by an orbit having as ω-limit the equilibrium (0, 0, 0). The orbit
on the x-axes near of (0, 0, 0) has this equilibrium as its ω-limit. Similar studies
can be done for the equilibria located at the origin of the local charts U2 and U3

that is, at the end of the y- and z-axes, respectively. Hence the boundary of the
infinity of O1 is formed by a heteroclinic cycle formed by three equilibria coming
from the ones located at the end of the positive half-axes, and the three orbits
living on the three open arcs connecting these three points and contained in the
boundary of the infinity of O1. From the study of the global dynamics on the
Poincaré sphere in section 3 we see that there is an additional equilibria in the
interior of the heteroclinic cycle that is an stable attractor and it is the ω-limit set
of each point on the interior of the heteroclinic cycle. This completes the proof of
statement (b.1). The orientation of the heteroclinic cycle is reversed for the case
α ≥ 1.

If −2 < α < 1, as shown in section 3 there are three additional equilibria on
the open arcs connecting the end of the invariant planes. In the local chart U1

the middle equilibria in the x-axis is the ω-limit of the equilibrium at the end of
the plane z = 0 and of (0, 0, 0). The same happens on the y-axis. Similar studies
can be done for the y-axis on the local charts U2 and for the equilibrium at the
end of the z-axis in the local chart U3. Thus the boundary of the infinity of O1

is formed by a graph formed by six equilibria coming from the ones located at the
positive half-axes, and the six orbits living on the six open arcs connecting these six
points and contained in the boundary of the infinity of O1. From the study of the
global dynamics on the Poincaré sphere made in section 3 we see that there is an
additional equilibrium point in the interior of the graph. It is an stable attractor
and it is the ω-limit set of each point on the interior of the graph. This completes
the proof of statement (b.2).

The proof concerning the ω-limit part of statement (c.1) follows directly from
Proposition 5.2. Combining the results from sections 3 and 4.1 we get the phase
portraits in Figure 3. From them and the existence of the separatrix cycle F we
conclude the proof of statement (c.1).

The proof concerning the ω-limit part of statement (c.2) follows directly from
Proposition 5.3. Combining the results from sections 3 and 4.1 we get the phase
portraits in Figure 4. From them and the existence of the graph G we conclude the
proof of statement (c.2). Hence, the proof of Theorem 2.2 is concluded. �

Proof of Theorems 2.3–2.5. Combining the analysis of the dynamics in the Poincaré
sphere and in each invariant plane obtained in sections 3 and 4.1, respectively, we
get the global phase portraits of system (2.3) in each octant given in Figures 4-6.
We note that these are all possible phase portraits up to a rotation (see Theo-
rem 2.1 (b)). The behavior in forward time of all orbits in each octant described
in Theorems 2.3 to 2.5 follow directly from the previous analysis and the Poincaré



16 R. OLIVEIRA, C. VALLS EJDE-2020/55

Bendixson theorem (see again Figures 4-6). The proof of Theorems 2.3 to 2.5 is
complete. �

7. Appendix

7.1. Poincaré compactification in R2. Let X = (P1(x, y), P2(x, y)) be a planar
polynomial vector field of degree n = max{deg(Pi) : i = 1, 2}. The Poincaré
compactified vector field p(X) corresponding to X is an analytic vector field induced
on S2 as follows (for more details, see [4]).

Let
S2 = {y = (y1, y2, y3) ∈ R3 : y21 + y22 + y23 = 1}

and TyS2 be the tangent plane to S2 at point y. We identify R2 with T(0,0,1)S2
and we consider the central projection f : T(0,0,1)S2 = S2. The map f defines two

copies of X on S2, one in the southern hemisphere and the other in the northern
hemisphere. Denote by X̄ the vector field D(f ◦ X) defined on S2 \ S1, where
S1 = {y ∈ S2 : y3 = 0} is identified with the infinity of R2.

For extending X̄ to a vector field on S2, including S1, X must satisfy convenient
conditions. Since the degree of X is n, p(X) is the unique analytic extension of
yn−13 X̄ to S2. On S2 \ S1 there is two symmetric copies of X, and once we know
the behavior of p(X) near S1, we know the behavior of X in a neighborhood of the
infinity. The Poincaré compactification has the property that S1 is invariant under
the flow of p(X). The projection of the closed northern hemisphere of S2 on y3 = 0
under (y1, y2, y3) 7→ (y1, y2) is called the Poincaré disc, and its boundary is S1.

Two polynomial vector fields X and Y on R2 are topologically equivalent if there
exists a homeomorphism on S2, preserving the infinity S1, carrying orbits of the
flow induced by p(X) into orbits of the flow induced by p(Y ) preserving or not the
orientation of all orbits.

As S2 is a differentiable manifold, in order to compute the explicit expression of
p(X), we consider six local charts

Ui = {y ∈ S2 : yi > 0} and Vi = {y ∈ S2 : yi < 0},
where i = 1, 2, 3, and the diffeomorphisms Fi : Ui = R2 and Gi : Vi = R2, for i =
1, 2, 3, which are the inverses of the central projections from the tangent planes at
the points (1, 0, 0), (−1, 0, 0), (0, 1, 0), (0,−1, 0), (0, 0, 1) and (0, 0,−1), respectively.
After some computations and a rescaling of the time, p(X) in the local charts U1

and U2 is given, respectively, by:

zn2 (P2 − z1P1,−z2P1) , where Pi = Pi(1/z2, z1/z2),

zn2 (P1 − z2P2,−z1P2) , where Pi = Pi(z1/z2, 1/z2).

The expression for p(X) in U3 is zn2 (P1, P2) and the expression for p(X) in Vi’s are
the same as that for Ui’s but multiplied by the factor (−1)n−1. In these coordinates
z2 = 0 always denotes the points of the infinity S1.

7.2. Poincaré compactification in R3. Let

X =
(
P1(x, y, z), P2(x, y, z), P3(x, y, z)

)
be a planar polynomial vector field of degree n = max{deg(Pi) : i = 1, 2, 3} and let

S3 = {y = (y1, y2, y3, y4) : ‖y‖ = 1},
S+ = {y ∈ S3 : y4 > 0}, S− = {y ∈ S3 : y4 < 0}



EJDE-2020/55 GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM 17

be, respectively, the unit sphere in R4, the northern hemisphere of S3 and the
southern hemisphere of S3. The tangent space of S3 at the point y will be denoted
by TyS3 and the tangent plane

T(0,0,0,1)S3 = {(x1, x2, x3, 1) ∈ R4 : (x1, x2, x3) ∈ R3}
can be identified with R3.

Consider the central projections f± : R3 = T(0,0,0,1)S3 → S± given by

f±(x) = ± (x1, x2, x3, 1)

∆(x)
with ∆(x) =

(
1 +

3∑
i=1

x2i

)1/2
.

Using these central projections, R3 is identified with S+ and S−. Note that the
equator of S3 is S2 = {y ∈ S3 : y4 = 0}.

The maps f± define two copies of X on S3, one Df+ ◦X in S+, and the other,
Df− ◦X in S−. Denote by X the vector field on S3 \ S2 = S+ ∪ S−, that restricted
to S+ coincides with Df+ ◦X, and restricted to S− coincides with Df− ◦X. We
can extend analytically the vector field X(y) to the whole sphere S3 setting p(X) =
yn−14 X(y). This extended vector field p(X) is called the Poincaré compactification
of X on S3.

Using that 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.

Note that the diffeomorphisms Fi : Ui → R3 and Gi : Vi → R3 for i = 1, 2, 3, 4 are
the inverse of the central projections from the origin to the tangent hyperplane at
the points (±1, 0, 0, 0), (0,±1, 0, 0), (0, 0,±1, 0) and (0, 0, 0,±1), respectively.

Assume that (0, 0, 0, 0), (y1, y2, y3, y4) ∈ S3 and (1, z1, z2, z3) in the tangent hy-
perplane to S3 at (1, 0, 0, 0) are collinear. Then we have 1/y1 = z1/y2 = z2/y3 =
z3/y4 and, so

F1(y) = (y2/y1, y3/y1, y4/y1) = (z1, z2, z3)

defines the coordinates on U1. As

DF1(y) =

−y2/y21 1/y1 0 0
−y3/y21 0 1/y1 0
−y4/y21 0 1/y1 0


and yn−14 = (z3/∆(z)n−1), the analytical vector field p(X) in the local chart U1

becomes, after a rescaling of the time variable,

zn3
(
− z1P1 + P2,−z2P1 + P3, z3P1

)
, where Pi = Pi(1/z3, z1/z3, z2/z3).

Similarly, the expressions of p(X) in U2 and U3, after a rescaling of the time variable,
are, respectively,

zn3
(
− z1P2 + P1,−z2P2 + P3, z3P2

)
, where Pi = Pi(z1/z3, 1/z3, z2/z3),

zn3
(
− z1P3 + P1,−z2P3 + P2, z3P3

)
, where Pi = Pi(z1/z3, z2/z3, 1/z3).

The expression for p(X) in U4 is zn+1
3 (P1, P2, P3) and the expression for p(X) in Vi

is the same as in Ui multiplied by (−1)n−1, for all i = 1, 2, 3, 4.
From now on we will consider only the orthogonal projection of p(X) from S+

to y4 = 0 and we will denote it again by p(X). Observe that the projection of the
closed S+ is a closed ball of radius one, denoted by B, whose interior is diffeomorphic
to R3. Its boundary, S2, corresponds to the infinity of R3. Moreover, p(X) is defined



18 R. OLIVEIRA, C. VALLS EJDE-2020/55

in the whole closed ball B in such way that the flow on the boundary, given by
z3 = 0 is invariant. The vector field induced by p(X) on B is called the Poincaré
compactification of X and B is called the Poincaré sphere.

We recall that two polynomial vector fields X and Y on R3 are topologically
equivalent if there exists a homeomorphism on S3, preserving the infinity S2, car-
rying orbits of the flow induced by p(X) into orbits of the flow induced by p(Y )
preserving or not the orientation of the orbits.

Acknowledgements. R. Oliveira was partially supported by FAPESP grants “Pro-
jeto Temático” 2019/21181-0 and project number 2017/20854-5. C. Valls was par-
tially supported by FCT/Portugal through UID/MAT/04459/2019.

References

[1] A. Battauz, F. Zanolin; Coexistence states for periodic competitive Kolmogorov system, J.
Math. Anal. Appl. 219 (1998), 179–199.

[2] G. Blé, V. Castellanos, J. Llibre, I. Quilantán; Integrability and global dynamics of the May-

Leonard model, Nonlinear Anal. Real World Appl. 14 (2013), 280–293.
[3] F. H. Busse; Transition to turbulence via the statistical limit cycle rout, Syneretics, Springer-

Verlag, 1978.

[4] F. Dumortier, J. Llibre, J. C. Artés; Qualitative Theory of Planar Differential Systems,
UniversiText, Springer-verlag, New York, 2006.

[5] P. Glansdorff, I. Prigogine; Thermodynamic theory of structure, stability and fluctuations,

John Wiley & Sons Ltd, London 1971.
[6] R. Hannesson; Optimal harvesting of ecologically interdependent fish species, J. Environ.

Econ. Manag. 10 (1983), 329–345.
[7] J. Jimenez, J. Llibre, J. C. Medrado; Crossing limit cycles for a class of piecewise linear

differential centers separated by a conic, Electron. J. Differential Equations 2020 no. 41

(2020), 1–36.
[8] A. Kolmogorov; Sulla teoria di Volterra della lotta per l’esistenza, G. Ist. Ital. Degli Attuari

7 (1936), 74–80.

[9] G. Laval, R. Pellat; Plasma Physics. Proceedings of Summer School of Theoretical Physics,
Gordon and Breach, New York, 1975.

[10] J. Llibre, R. Oliveira; Quadratic systems with invariant straight lines of total multiplicity

two having Darboux invariants, Commun. Contemp. Math. 17 (2015), 1450018, 17 pp.
[11] J. Llibre, X. Zhang; Dynamics of some three-dimensional Lotka-Volterra systems, Mediter-

ranean J. Math. 14 (2017), 126–139.

[12] A. J. Lotka; Analytical note on certain rhythmic relations in organic systems, Proc. Natl.
Acad. Sci. USA, 6 (1920), 410–415.

[13] R. M. May; Stability and Complexity in Model Ecosystems, Princeton, New Jersey, 1974.
[14] R. M. May, W. J. Leonard; Nonlinear aspects of competition between three species, SIAM J.

Appl. Math., 29 (1975), 243–253.

[15] L. P. Peng, Z. Feng; Limit cycles from a cubic reversible system via the third-order averaging
method, Electron. J. Differential Equations, 2015 no. 111 (2015), 1–27.

[16] P. Schuster, K. Sigmund, R. Wolf; On ω-limits for competition between three species, SIAM

J. Math., 37 (1979), 49–54.
[17] V. Volterra; Lecons sur la Théorie Mathématique de la Lutte pour la vie, Gauthier-Villars,

Paris, 1931.

Regilene Oliveira

Departamento de Matemática, ICMC-Universidade de São Paulo, Avenida Trabalhador

São-carlense, 400 - 13566-590, São Carlos, SP, Brazil
Email address: regilene@icmc.usp.br



EJDE-2020/55 GLOBAL DYNAMICS OF THE MAY-LEONARD SYSTEM 19

Claudia Valls

Departamento de Matemática, Instituto Superior Técnico, Universidade de Lisboa, 1049-

001 Lisboa, Portugal
Email address: cvalls@math.tecnico.ulisboa.pt


	1. Introduction and statement of main results
	2. Statements of main results
	3. Dynamics on the Poincaré sphere S2
	4. Dynamics of system (2.3) restricted to the invariant planes
	5. On the -limit of the orbits of system (2.3) on O1
	6. Proofs of Theorems 2.1–2.5 
	7. Appendix
	7.1. Poincaré compactification in R2
	7.2. Poincaré compactification in R3
	Acknowledgements

	References