DYNAMICS OF SOME THREE DIMENSIONAL LOTKA–VOLTERRA SYSTEMS

We characterize the dynamics of the following two Lotka– Volterra differential systems ẋ = x(r + ay + bz), ẋ = x(r + ax+ by + cz), ẏ = y(r − ax+ cz), and ẏ = y(r + ax+ dy + ez), ż = z(r − bx− cy), ż = z(r + ax+ dy + fz). We analyze the biological meaning of the dynamics of these Lotka– Volterra systems.


Introduction and statement of the main results
The Lotka-Volterra differential systems appeared at the work of A.J. Lotka in 1910 for modeling autocatalytic chemical reactions, and was extended by himself in 1920 to for modeling the dynamics of a plant species and a herbivorous animal species, where x and y are respectively the numbers of preys and predators, and α, β, γ and δ are positive real parameters describing the interaction of the two species. Volterra developed this last model independently from Lotka to explain the exchange of the fish catches between fish and predatory fish in the Adriatic Sea during the first World War. Kolmogorov [9] in 1936 studied these systems again and extended them to arbitrary dimension, and for this reason these kinds of systems are also called Kolmogorov systems.
This last Lotka-Volterra differential system has been modified in different ways for studying the dynamics of the interaction between the competition of two or more species. Lotka-Volterra systems has also been used to model dynamical phenomena from different subjects, such as hydrodynamics [3], plasma physics [10], chemical reactions [7], and evolution of conflicting species in biology [8,18] and so on. Lotka-Volterra systems have been studied from different points of view, see for instance [1,2,4,11,12,13,14,17].
Here we study two kinds of differential systems of Lotka-Volterra type. The first one is the next system.
x = x(r + ay + bz), where the dot denotes the derivative with respect to the time t, and a, b, c are nonzero constants, and (x, y, z) are located in the positive octant of R 3 , here this means the set {(x, y, z) ∈ R 3 : x ≥ 0, y ≥ 0, z ≥ 0}. Under the change of variables where we still use the variables x, y, z instead of X, Y, Z.
Our first result characterizes the global dynamics of system (3) in the positive octant of R 3 . Theorem 1. For system (3) the following statements hold.
(a) System (3) has the two functionally independent first integrals H 1 = x + y + z and H 2 = x c y −b z a . (b) If abc = 0 then the global phase portraits of system (3) on the invariant set {(x, y, z)| x + y + z = h, x ≥ 0, y ≥ 0, z ≥ 0} are topologically equivalent either to the one of Fig. 1, or to the one of Fig. 2.
The biological meaning of the conclusion in Fig. 1 of statement (b) is that only the species x will survive. Whereas the biological meaning of the conclusion in Fig. 2 of statement (b) is that the three species persist on a periodic solution, or in an equilibrium point.
The proof of Theorem 1 is given in section 2. We remark that if abc = 0, system (3) has less biological meaning and its dynamics can be easily obtained by our methods, and so it is omitted here.
We now study the global dynamics of system (3) in the positive octant of R 3 . Theorem 2. Any orbit of system (1) when r = 0 starting at some point of the positive octant except at the origin, either goes in forward time to the origin and in backward time approaches to infinity, or vice-versa.
The proof of Theorem 2 will be given in section 3. 3 We note that the dynamics of systems (1) and (3) are completely different, but they are related through the transformation (2). Each orbit of system (3) is bounded and it is limited by one of the invariant triangles L + h defined at the beginning of section 2. If r > 0, this orbit under the inverse change of variables (2) goes to the origin when t → −∞ and to infinity when t → ∞. If r < 0, the converse happens.

DYNAMICS OF THREE DIMENSIONAL LOTKA-VOLTERRA SYSTEMS
Next we consider the following three dimensional differential systems of Lotka-Volterra type.ẋ = x(r + ax + by + cz), y = y(r + ax + dy + ez), In a similar way to the study of system (1), under the transformation (2) system (4) becomesẋ = x(ax + by + cz), y = y(ax + dy + ez), where again we still use the variables x, y, z instead of X, Y, Z. In order to avoid degeneracies we assume that For system (5) we have the next result.  Since statement (b) provide the α-and ω-limits sets of all the orbits inside the positive octant of R 3 , we can determine the initial and final evolutions of the three species modeled by the Lotka-Volterra system (5).
The proof of Theorem 3 will be given in section 4. The heteroclinic orbits in (c) will be precisely described in the proof of that statement.
For a definition of Poincaré disc and Poincaré sphere see [6] and [5,15], respectively. When we say the singularities at the endpoints of the axes, we mean the singular points which are on the boundary of the Poincaré disc or of the Poincaré sphere at the endpoints of the coordinate axes.

Proof of Theorem 1
Proof of statement (a). It can be verified by direct calculations.
Proof of statement (b). For studying the global dynamics of system (3) we consider each level surface L h := {(x, y, z) ∈ R 2 | x + y + z = h} and denote by L + h its restriction to the positive octant. So we must have h ≥ 0. For h = 0 the level surface is limited to the origin. For h > 0, L + h is a triangle with the three boundaries denoted by B x h , B y h and B z h , which are invariant and located respectively on the yz, xz and xy planes. This follows from the invariance of the three coordinate planes and of the level surface L h . The three vertices of L + h are denoted by P x h = (h, 0, 0), P y h = (0, h, 0) and P z h = (0, 0, h), which are singularities of system (3) and are located respectively on the x, y and z axes.
According to the above analysis we only need to study global dynamics of system (3) on L + h . Permuting the names of the variables and changing the sign of the time (if necessary), we only need to study two cases: • Case 1: a, b and c are positive, • Case 2: a and c positive and b negative.

Restricting system (3) to the invariant set L +
h we obtain Beside the singularities P x h , P y h and P z h , the necessary condition for system (6) to have a fourth singularity in L + h is a − b + c = 0. If it exists, it should be of the form This shows that system (6) has a fourth singularity in L + h if and only if the case 2 happens.
The Jacobian matrix of system (6) is It is easy to show that the singularities S x h = (h, 0), S y h = (0, h) and S z h = (0, 0) of (6), associated to P x h , P y h and P z h respectively, have the eigenvalues (−ah, −bh), (ah, −ch), (bh, ch),  In case 1 system (6) has only the three singularities S x h , S y h and S z h , which are respectively stable node, saddle and unstable node. This shows that system (3) has the phase portrait on L + h ⊂ {x + y + z = h} given in Fig. 1. In case 2 the three singularities S x h , S y h and S z h are all saddles. The singularity S + h of system (6), associated to P + h in the interior of L + h , has a pair of pure imaginary eigenvalues In order to distinguish if this singularity is either a focus or a center, we compute a first integral of system (6). From statement (a) system (6) has the Darboux first integral For more details on Darboux theory of integrablity, see for example [6,Chapter 8] or [16]. Clearly the Darboux first integral H + h is an analytic function in the interior of the triangle Γ + h , which is projection of L + h onto the xy plane. By the classical Poincaré-Lyapunov theorem which says that an elementary monodromy singularity of a planar analytic differential system is a center if and only if it has an analytic first integral defined in a neighborhood of the singularity, it follows that the singularity S + h is a center. Moreover, we can prove using the first integral H + h that the periodic orbits of system (6) fill up the interior of L + h except the singularity S + h . So system (3) has the phase portrait in L + h given in Fig. 2.

Proof of Theorem 2
Note that r = 0. We can check that system (1) has the three invariant planes x = 0, y = 0 and z = 0, and has always the two functionally independent Darboux invariants Recall that a Darboux invariant of a polynomial vector field Y(w) in R n is a function of the form e −σt f (w) with σ a nonzero constant and f a polynomial satisfying Y(f (w)) = σf (w), i.e. f is a Darboux polynomial of the vector field Y(w) with cofactor σ.
We note that if a − b + c = 0, the function V 2 is a first integral. From these last two invariants V 1 and V 2 it follows that system (1) has always the first integral F (x, y, z) = x c y −b z a (x + y + z) −(a−b+c) .
Some easy calculations show that system (1) has a unique singularity in the positive octant, the origin.
Since the origin is the unique finite singularity of system (1), and the three coordinate planes are invariant, it follows that any orbit starting from the interior of the positive octant will stay in it.
For proving the theorem, let γ + be an orbit with its initial point, say P + , located in the positive octant except at the origin, and let (x, y, z) = (ξ(t), η(t), ζ(t)) be the expression of this orbit. By the Darboux invariant V 1 we have That is ξ(t) + η(t) + ζ(t) = he rt .

DYNAMICS OF THREE DIMENSIONAL LOTKA-VOLTERRA SYSTEMS 7
Then the theorem follows from this last expression and the facts that ξ(t), η(t), ζ(t) > 0.

Proof of Theorem 3
We can check that system (5) has the first integral which is of Darboux type. This first integral is difficult to use for studying dynamics of system (5). But we can check that By assumption (H 0 ) we can assume without loss of generality that e > f, otherwise we can interchange the variables y and z. Some easy calculations show that system (5) has only singularities located on the three invariant coordinate planes. In addition, since y = 0 is invariant, we get from (7) that each orbit starting at the point outside y = 0 will transversally pass through all the planes y/z =constant = 0. This implies that any orbit starting at the point outside the three invariant coordinate planes will either approach to the coordinate planes, or approach to infinity. So we focus our next studies on the coordinate planes and the infinity.
Under the assumption (H 0 ), after rescaling the spacial variables

system (5) can be written in an equivalent way aṡ
x = x(x + y + z), y = y(x + dy + ez), For studying the dynamics of system (8) at infinity, we use the Poincaré compactification, see [5]. Using the local chart of the infinity at the x direction, under the change of variables and the time rescaling dt = wdτ , we get where the prime denotes the derivative with respect to τ . By assumption (H 0 ) and the rescaling we have d = 1, e = 1 and f = e. If f = 1, system (9) has a unique singularity on the invariant plane w = 0, i.e. the origin, and has an one dimensional stable manifold (i.e. the w axis), and a two dimensional center manifold (i.e. the plane at infinity). Combining the directions of orbits on the y and z axes and the fact that we can get the local dynamics of system (9) at the origin, it is shown in one of the first four phase portraits of Fig. 3. If f = 1, system (9) has the line {w = 0} ∩ {y = 0} fulfilled with singularities. Since e > f one have e > 1. Distinguishing d > 1 and d < 1 we get the last two phase portraits of Fig.  3.
Using the local chart at infinity in the y direction, under the change of variables x → x w , y → 1 w , z → z w , and the time rescaling dt = wdτ , we get z = (f − e)z 2 .
Since d = 1, e = 1 and e > f, system (10) has the unique singularity at the origin, which is a saddle-node on w = 0. System (10) has the local  phase portraits given in Fig. 4. Using the local chart at infinity in the z direction, under the change of variables x → x w , y → y w , z → 1 w , and the time rescaling dt = wdτ , we get y = (e − f )y.
If f = 1, system (11) has the unique singularity at the origin, which is hyperbolic. System (11) has either the first, or the third, or the fourth local phase portrait given in Fig. 5. If f = 1, on w = 0 the line y = 0 is fulfilled with singularities. Since e > f = 1 we have the second phase portrait of Fig. 5.
Combining the above analysis in each local charts at infinity, we get the global phase portraits of system (8) in the Poincaré sphere. Fig. 6 shows the phase portraits of system (8) on the semi-sphere, which faces us. The phase portraits on the back side of the Poincaré sphere is a projection of those in Fig. 6. We complete the proof of statement (a). Figure 6. Phase portraits of system (8) on the Poincaré sphere. To prove statement (b), we restrict our study to each coordinate plane. On the invariant plane x = 0, system (8) is reduced to (12)ẏ = y(dy + ez),ż = z(dy + f z).
This system has the same form as system (9) when restricted to w = 0. So they have the same phase portraits as those in Fig. 5 without the w axis, here we have used the fact that d, e, f, e − f > 0. Then the phase portraits of system (12) on the Poincaré disc are the same as those in Fig. 6 only with different parameter conditons, see Fig. 7.
It has the same formula than (13) and so has either the first, or the third, or the fourth phase portrait of Fig. 8 because d = 1 by assumption. For a better understanding the global phase portraits of system (8), we draw the phase portraits of system (14) in Fig 9. This proves statement (b).  This shows that all orbits with the initial points in the interior of the positive octant are heteroclinic. We complete the proof of statement (c) and consequently the theorem.