Manuscript submitted to DCDS Supplement Volume 2007

Website: www.aimSciences.org pp. 1–10

PERIODIC ORBITS OF TRITROPHIC SLOW–FAST SYSTEM AND DOUBLE HOMOCLINIC BIFURCATIONS

Alexandre VIDAL Universit´ e Pierre et Marie Curie-Paris6, Paris, F-75005 France Laboratoire J.L. Lions, UMR 7598 CNRS 175 rue du Chevaleret 75013 Paris, France Abstract. The biological models - particularly the ecological ones - must be understood through the bifurcations they undergo as the parameters vary. However, the transition between two dynamical behaviours of a same system for diverse values of parameters may be sometimes quite involved. For instance, the analysis of the non generic motions near the transition states is the first step to understand fully the bifurcations occurring in complex dynamics. In this article, we address the question to describe and explain a double bursting behaviour occuring for a tritrophic slow–fast system. We focus therefore on the appearance of a double homoclinic bifurcation of the fast subsystem as the predator death rate parameter evolves. The first part of this article introduces the slow–fast system which extends Lotka– Volterra dynamics by adding a superpredator. The second part displays the analysis of singular points and bifurcations undergone by fast dynamics. The third part is devoted to the flow analysis near the homoclinic points. Finally, the fourth part is concerned with the main results about the existence of periodic orbits of different periods as the two homoclinic orbits are close enough to each other.

1. Introduction. Let us consider the following system, which is built from a natural plane dynamics:     X P1 Y dX = X R 1− − dT K S1 + X   dY P1 X P2 Z = Y E1 − D1 − (1) dT S1 + X S1 + Y   dZ P2 Y = εZ E2 − D2 dT S2 + Y where X, Y and Z represent the membership of three populations. The variable X is the prey, Y its predator and Z a superpredator for Y . The threshold constant K and the intrinsic growth rate of the prey R characterize the logistic evolution of X (see [6]).

2000 Mathematics Subject Classification. Prim: 34C05, 34C25, 34C26, 34C37; Second: 92D25. Key words and phrases. Slow–Fast System, Tritrophic System, Homoclinic Bifurcation, Periodic Orbit.

1

2

ALEXANDRE VIDAL

The predator–prey interactions are described by two Holling type II matching defined by the following positive parameters: Pj : the maximum predation rates Sj : the half-saturation constants Dj : the death rates Ej : the efficiencies of predation j = 1 concerning the predator Y , j = 2 concerning the superpredator Z. Finally, it is natural to consider that the evolution of the superpredator is slower than that of the preys. Then we introduce different time scales by means of the constant ε << 1. Thus, in the following, we adopt the point of view and methods of singular perturbations theory applied to slow–fast systems. In order to obtain a simpler and more useful analytic form, Klebanoff and Hastings proposed in [4] the following rescalings: x=

X , K

y=

Y , KE1

z=

Z , KE1 E2

t = RT

(2)

which yields (Sε ):   a1 y = f (x, y, z) x˙ = x 1 − x − 1 + b1 x   a1 x a2 z y˙ = y − d1 − = g(x, y, z) 1 + b1 x 1 + b2 y   a2 y − d2 = h(x, y, z) z˙ = εz 1 + b2 y

(3)

All axes and faces of the positive octant R3+ are invariant sets of the system (Sε ). Thus, we limit the phase space to this positive octant. We use also the “reduced system” (S0 ), obtained by posing ε = 0. It describes, for each value of z considered then as a parameter, the behaviour of the “fast system” (Fz ) formed by the dynamics on x and y. Moreover, it is useful to consider the “critical system” obtained from (Sε ) after the following rescaling: τ = εt

;

t=

τ ε

;

dτ =ε dt

and, thereafter, by setting ε = 0. It is written as:   a1 y f (x, y, z) = x 1 − x − =0 1 + b1 x   a1 x a2 z g(x, y, z) = y − d1 − =0 1 + b1 x 1 + b2 y   a2 y z˙ = z − d2 1 + b2 y

(4)

(5)

and is defined on the “critical set”: C = {(x, y, z) ∈ R3+ |f (x, y, z) = 0, g(x, y, z) = 0} Note that C is the set of singular points of the reduced system (S0 ).

(6)

PERIODIC ORBITS OF TRITROPHIC SLOW–FAST SYSTEM

3

2. Singular Points. We later assume a1 − d1 (b1 + 1) > 0 so that the steady point: ! d1 a1 − d1 (b1 + 1) M = (xM , yM , 0) = ,0 (7) , a1 − b1 d1 (a1 − b1 d1 )2 lays in the phase space R3+ . 2.1. Critical Set Analysis. As z, considered as parameter, varies, the number of critical points of the fast dynamics (Fz ) evolves. To explain its dependency in terms of z, it is more convenient to consider the critical set C defined in (6). One obtains easily that C = ∆ ∪ L where: ∆

= {(1, 0, z)|z ∈ R+ }  (x, yL (x), zL (x)) ∈ R3+ |x ∈ [0, 1]

(8)

L =

(9)

and: 1 (1 − x) (1 + b1 x) a1 (a1 x − d1 (b1 + 1)) (a1 + b2 (1 − x)(1 + b1 x)) zL (x) = a1 a2 (1 + b1 x)

yL (x) =

(10) (11)

Note that ∆ always intersect L at the point T = (1, 0, zT ) where: zT =

(a1 − d1 (b1 + 1)) a1 d1 = − >0 a2 (1 + b1 ) a2 (1 + b1 ) a2

(12)

Moreover, we assume that d1 is small enough such that: 0 ∃xP > 0, zL (xP ) = 0

(13)

Consequently, L is ∩-shaped in R3+ and we denote: L - for each z ∈ ]0, zP [, Rz ± the unique point (x, yL (x), z) ∈ L such as x < xP , - for each z ∈ ]zT , zP [, RzLS the unique point (x, yL (x), z) ∈ L such as x > xP . ∆ Furthermore, we name Rz∆S the point (1, 0, z) ∈ ∆ for each 0 < z < zT and Rz − for each z > zT . Then, let us pose (see figure 1): S S ∆ ∆− = ∆ ∩ {z > zT } = {Rz − }, ∆S = ∆ ∩ {z < zT } = {Rz∆S } z>zT 0 xP } = {RzLS }, L± = L ∩ {x < xP } = {Rz ± } zT
0
Thus, P = (xP , yL (xP ), zL (xP )) splits L into LS and L± . Note that the critical set C is the closure of the union of the smooth manifolds ∆S , ∆− , LS and L± .  We introduce the projection π : R3 → R2 ; (x, y, z) → (x, y) . Then, one may easily verify that (Fz ) admits the following singular points: L - if 0 < z < zT : π(Rz ± ) and the saddle π(Rz∆S ) whose attractive direction is x, L ∆ - if zT < z < zP : π(Rz ± ), the attractive node π(Rz − ) and the saddle π(RzLS ), ∆− - if z > zP : the attractive node π(Rz ). Finally, for d1 small enough: ( L 0 < z < zH =⇒ π(Rz ± ) is a repulsive focus ∃zH ∈ ]0, zP [ , (14) L zH < z < zP =⇒ π(Rz ± ) is attractive (node or focus)

4

ALEXANDRE VIDAL

Figure 1. Critical set and its different parts according to the nature of the singular points of the fast dynamics. The heteroclinic orbit γP of (Fz ) in {z = zP } and the orbit γR connecting R and CzR in {z = zR }. Then, H = (xH , yL (xH ), zH ) splits L according to the stability of its points for the fast dynamics: [ [ {RzL− } L+ = L± ∩ {z < zH } = {RzL+ }, L− = L± ∩ {z > zH } = zH
0
2.2. Bifurcation analysis. We can now describe the bifurcations of the fast system according to z and d1 for a1 , b1 , a2 , b2 , d2 fixed. The singular point (0, 1) undergoes a saddle–node transcritical bifurcation for (Fz ) on the line: T :z=

a1 d1 + > 0, a2 (1 + b1 ) a2

d1 > 0

(15)

For z < zL (xP ) great enough, (Fz ) displays the saddle π(RzLS ) and an attractive L node π(Rz − ) which disappear, as z increases, through a saddle–node bifurcation when: P : z = zL (xP ) (16) L

A Hopf bifurcation may transform the repulsive foci π(Rz + ) into the attractive L focus π(Rz − ) and it occurs on the curve: H : hp (z, d1 ) = 0

(17)

that is z = zH which depends on d1 and is defined by 14. It generates limit cycle Cz L of (Fz ), at least for z < zH close to zH , which surrounds the singular point π(Rz + ). Definition 1. When the limit cycles Cz exist for all values z1 < z < z2 , the set: [ Mzz21 = Cz z1
is a smooth two-dimensional manifold.

PERIODIC ORBITS OF TRITROPHIC SLOW–FAST SYSTEM

5

For d1 small enough, the limit cycle may disappear for zH1 < z < zH2 < zH through two homoclinic bifurcations: CzH1 for the saddle π(H1 ) where H1 = RzLHS1 (“direct”) and CzH2 for the saddle π(H2 ) where H2 = RzLHS2 (“reverse”). This shows the decomposition into connected components of the invariant manifold into zH z M0H1 ∪ Czh1 ∪ Czh2 ∪ MzHp1 . The whole occurrence of the bifurcations of (Fz ) according to z and d1 is shown in figure 2.

Figure 2. Bifurcation diagram of the fast dynamics according to z and d1 . For d1 = dhc 1 , a single homoclinic bifurcation appear. For 0 ≤ d1 < dhc , the “direct” and “reverse” homoclinic bifurcations 1 appear for z = zH1 and z = zH2 . We have the following: Proposition 1. For any value of parameter z in ]zT , zP ], the fast system (Fz ) displays a heteroclinic orbit which connects the saddle π(RzLS ) to the attractive node ∆ π(Rz − ) = (1, 0). We denote by γP the particular heteroclinic orbit for z = zP . In the plane {z = zR }, for all 0 ≤ zR ≤ zT , if the focus on L+ is surrounded by a limit cycle CzR , the fast system (Fz ) displays an orbit γR which connects R = (1, 0, zR ) to the attractive limit cycle CzR (see figure 1). Proof. It is a direct consequence of the global attractivity of the node (resp. the limit cycle CzR ) under the flow of (Fz ) for the values ]zT , zP ] (resp. [0, zT ]) of parameter z. 3. Bursting oscillations. We are interested in the case where d1 < dhc 1 such that the two homoclinic bifurcations occur for z = zH1 and z = zH2 close enough (zH2 − zH1 < γ where γ will be introduced in the main result). As z increases, the saddle–node transcritical bifurcation of (0, 1) first appears , then the two homoclinic bifurcations, the Hopf bifurcation and finally the saddle-node bifurcation. 3.1. Qualitative description of phase portrait. The typical motion of an or−1 bit starting near the branch ∆− is as follows. As z˙ ≤ 0 for y ≤ (b2 d2 − a2 ) , z decreases near ∆. The orbit reaches the neighborhood of T and then follows the unstable center manifold of ∆S . It joins the neighborhood of the normally hyperbolic manifold formed by limit cycles persisting from the homoclinic bifurcation at

6

ALEXANDRE VIDAL z

H1 . Then, as z increases, the orbit follows while turning around M0H1 . Near the boundary CzH1 , the orbit follows the branch of the unstable central manifold of LS . Under the influence of the fast dynamics, it approaches ∆− again. But another behaviour near CzH1 may also occur. Let W0s (resp. W0u ) be the union of all the stable (resp. unstable) manifolds of the singular points RzLS for zT < z < zP of (S0 ). These are smooth two dimensional invariant manifolds which perturb, for ε > 0, respectively into the stable center manifold Wεs and the unstable center manifold Wεu . Let us pose Wεc = Wεs ∩ Wεu the center manifold. z The orbit rises turning around the manifold M0H1 until z = zH1 and may lie precisely on Wεs . Then, it tracks for a while near the middle branch Wεc slowly going up near LS as in figure 3. So, when z = zH2 , the orbit spirals again this zH time around the upper cup MzH1P and next follows the branch L− formed by stable singular points of (Fz ). After this, it passes near the “junction point” P and, under the influence of the fast dynamics, approaches ∆− again.

Figure 3. Two types of Bursting orbit. The first one joins the vicinity of ∆− as soon as it passes close to H1 . The second one follows the stable center manifold for a while, spirals again around the upper cup and passes close to P before joining the vicinity of ∆− . To characterize the occurrence of such orbits, the key point of the analysis is to understand how the orbit may follow the stable center manifold Wεs . The next subsection points to characterize the properties of the flow and the choice of initial data which leads to this situation. 3.2. Flow Analysis near H1 . We can choose local coordinates near the homoclinic point H1 for which the system (Sε ) is linearized as follows: x˙ = λx x y˙ = −λy y z˙ = λz ε

(18)

where λx , λy and λz are positive contants. As the homoclinic orbit bifurcates into stable limit cycles, we may assume λy > λx .

PERIODIC ORBITS OF TRITROPHIC SLOW–FAST SYSTEM

7

We analyse the flow in the rectangular neighborhood of H1 : N = {(x, y, z) ∈ R3 | |x| ≤ α, |y| ≤ β, |z| ≤ γ}

(19)

where α, β and γ do not depend on ε. We define the top, front, back and right sides of this rectangle: ΣT = {(x, y, z) ∈ N |z = γ}, ΣF = {(x, y, z) ∈ N |x = α} ΣB = {(x, y, z) ∈ N |x = −α}, ΣR = {(x, y, z) ∈ N |y = β} Then, every orbit which begin on the right side ΣR enter N and thereafter exit through ΣT or ΣF or ΣB . If the orbit leaves through ΣF , it keeps bursting and if it leaves through ΣB , it joins the vicinity of ∆− , as the green orbit in figure 4. Furthermore, a fine choice of α and β allows us to characterize the orbits which leave N through ΣT as the ones which follow Wεc for a sufficient while to burst again around the upper cup (see the red orbit in figure 4).

Figure 4. The neighborhood N of the homoclinic point H1 . The set U ⊂ ΣR is shown in grey. The green orbit does not pass through U. Thus, it joins the vicinity of ∆− as soon as it crosses z = zH1 . The red orbit passes through U. Thus, it escapes N through ΣT . Let us integrate the linearized system (18) for an initial data (x0 , β, z0 ) ∈ ΣR : x(t) = x0 exp(λx t) y(t) = β exp(−λy t) z(t) = z0 + λz εt Then, this orbit leaves N through ΣT if and only if:   1 λx (γ − z0 ) |x0 | < α exp − ε λz

(20)

8

ALEXANDRE VIDAL

This condition defines a set U included in ΣR . This set is shown in grey on the figure 4. 3.3. Construction of the return map near CzH1 . The trajectories which begin on Σ+ R = {(x, β, z) ∈ ΣR |x > 0} exit through ΣF or ΣT . So, let D(φε ) be the set of all points lying on Σ+ R whose trajectories exit through ΣF and φε : D(φε ) → ΣF be the corresponding map defined by the flow. Remark 1. lim D(φε ) = Σ+ R. ε→0

Remark 2. We can explicitly write: φε φyε (x, z)

(x, β, z) → (α, φyε (x, z), φzε (x, z)) λz α x λy = β.( ) λx , φzε (x, z) = z + ε ln( ) α λx x :

(21) (22)

The homoclinic orbit CzH1 intersects ΣF at (α, 0, 0) and ΣR at (0, β, 0). Then, when ε = 0, the flow defines a diffeomorphism from a neighborhood of (α, 0, 0) in ΣF onto a neighborhood of (0, β, 0) in ΣR . We can find a fixed set D(ψ) independent of ε and containing (α, 0, 0) in ΣF such that the map ψε : D(ψ) → ΣR is a diffeomorphism from a subset of ΣF into ΣR for all ε small enough. Then we can define the following return map: πε = ψε ◦ φε : D(φε ) ⊂ ΣR → ΣR

(23) z M0H1

We look the trajectories crossing U backwards as it spiral around until they cross a section Σ = {y = δ} for the first time. We then may prove that: Proposition 2. For δ small independent of ε, the subset of Σ = {y = δ} which is the image of U by the function induced by the backward flow of (Sε ) is formed by strips whose width are exponentially small as ε → 0. 4. Periodic Orbits. In this part, we present results relating to the existence of periodic orbit in the case of Homoclinic Bursting. Limit–periodic sets have been introduced in [1] to which we refer for more details. It is sufficient here to say that they are the limits for the convergence via Hausdorff distance1 of periodic orbits families. 4.1. Existence and Unicity of Bursting Attractive Periodic Orbit. Definition 2. We define (see figure 5): ∆R1 −H1 = ∆ ∩ {zR1 < z ≤ zH1 } z

H1 R1 −H1 1 ΓR ∪ γR1 ∪ MzR ∪ γH1 1 =∆ 1

∆R2 −P = ∆ ∩ {zR2 < z ≤ zP } LSH1 −H2 = LS ∩ {zH1 ≤ z ≤ zH2 } z

zH S R2 −P 1 2 ΓR ∪ γR ∪ MzH R2 ∪ LH1 −H2 ∪ MzR2 ∪ L− ∪ γP 2 =∆

1 The Hausdorff distance between two compacts K and K 0 of a metric space (E, d) is the smallest r > 0 such that ∀x ∈ K, d(x, K 0 ) < r and ∀y ∈ K 0 , d(y, K) < r.

PERIODIC ORBITS OF TRITROPHIC SLOW–FAST SYSTEM

9

Figure 5. Two various limit periodic sets of the system. Theorem 1. There is γ0 such as for all γ < γ0 , there are ε0 , a decreasing real sequence (εi ) bounded from above by ε0 , and a decreasing real sequence (ci ) such that: ∀i, ci + ci+1 < εi − εi+1 (1) for any ε verifying: εi+1 + ci+1 < ε < εi − ci there exists a periodic orbit γε1 of system (Sε ) which is of saddle type. There is R1 = (1, 0, zR1 ) with 0 ≤ zR1 ≤ zT , such as, when ε tends to 0, the one– 1 parameter family (γε1 ) converges in the Hausdorff distance sense to ΓR 2 . (2) There exists a decreasing real sequence (di ) such that, for all i, di < ci , and for any ε verifying: εi − d i < ε < εi + d i there is a unique attractive periodic orbit γε2 of system (Sε ) passing in a neighborhood of P of O(ε2/3 )–diameter. There is R2 = (1, 0, zR2 ) with 0 ≤ zR2 < zT , such as, when ε tends to 0, the 2 one–parameter family (γε2 ) converges in the Hausdorff distance sense to ΓR 2 . For the values of ε different from those prescribed in the theorem, there are several attractive periodic orbits. Thus, one can define various one–parameter families (γε )0<ε<ε0 of periodic orbits. 4.2. Simulations. The figure 6 displays, for the same values of the parameters except ε, a periodic orbit of each type described above. As, for the chosen values of a1 , a2 , b1 , b2 and d2 , dhc 1 ' 0.08273891, we have chosen d1 = 0.08273890 to obtain zH1 ' 0.04727 and zH2 ' 0.04733. Then γ < 10−4 . We note that the two values of ε are close. This shows the high sensitivity to the parameters and the fact that the system is not structurally stable under such conditions. Finally, as explained in the theorem, we can exhibit each type of periodic orbit for values of ε < ε0 as small as we wish. But, according to Proposition 2, the intervals of values around εi for which the system admits a periodic orbit of the

10

ALEXANDRE VIDAL

second type become very small. Moreover, the basin of attraction of a such an orbit is exponentially small according to ε.

Figure 6. The two types of periodic orbits for same values of parameters except ε. (Simulations done using XPP, see http://www.math.pitt.edu/˜bard/xpp/xpp.html). REFERENCES [1] J.P. Fran¸coise, C.C. Pugh, Keeping track of limit cycles, Journal of Differential Equations, 65 (1986), 139–157. [2] J.P. Fran¸coise, Oscillations en biologie, Analyse qualitative et mod` eles, Coll.: Math´ ematiques et Applications (vol. 46), Springer, 2005. [3] K.D. Holl, B.B. Howarth, Paying for restoration, Restoration Ecology, 8 (2000), 260–267. [4] A. Klebanoff, A. Hastings, Chaos in three species food chains, Journal of Mathematical Biology, 32 (1994), 427–451. [5] Yu.A. Kutznetsov, S. Rinaldi, Remarks on food chain dynamics, Mathematical Biosciences, 134 (1996), 1–33. [6] J. Murray, Mathematical Biology, Springer, Berlin Heidelberg New York, 1993. [7] J. Rinzel, G.B. Ermentrout, Analysis of neuronal excitability and oscillations, In: Koch C and Segev I (eds) Methods in neuronal modeling MIT Press, Cambridge, MA., 1989 [8] D. Terman, E. Lee, Uniqueness and stability of periodic bursting solutions. Journal of Differential Equations, 158 (1999), no. 1, 48–78. [9] A. Vidal, Stable periodic orbits associated with bursting oscillations. Lecture Notes in Control and Information Sciences, 158 (2006), no. 341, 439–446.

Received September 2006; revised January 2007. E-mail address: [email protected]