# Prey-Predator Interactions in Two and Three Species Population Models

Authors: Arild Wikan; Ørjan Kristensen
Journal: Discrete Dynamics in Nature and Society
 (2019)

---

## Abstract

Discrete nonlinear two and three species prey-predator models are considered. Focus is on stability and nonstationary behaviour. Regarding the two species model, depending on the fecundity of the predator, we show that the transfer from stability to instability goes through either a supercritical flip or a supercritical Neimark-Sacker bifurcation and moreover that there exist multiple attractors in the chaotic regime, one where both species coexist and another where the predator population has become extinct. Sizes of basin of attraction for these possibilities are investigated. Regarding the three species models, we show that the dynamics may differ whether both predators prey upon the prey or if the top predator preys upon the other predator only. Both the sizes of stable parameter regions as well as the qualitative structure of attractors may be different.


---

## Body

## 1. Introduction

In 1924 and 1926, respectively, Lotka [1] and Volterra [2] independently established a two species prey-predator model which today is known under the name ‘the Lotka-Volterra prey-predator model’. The model consists of a system of two coupled nonlinear differential equations and as it is well known; the dynamical outcome of such a system is either a stable equilibrium or a limit cycle. Unfortunately, the Lotka-Volterra model has an undesired property; namely, it is structurally unstable, which in turn implies that most attempts to apply the model on real world phenomena are likely to fail. Therefore, after the pioneer works in the 1920s, there has been a tremendous development of prey-predator models. At first, most of these models were formulated in continuous time; see for example the work by Rosenzweig and MacArthur 1963 [3] and Holling 1965 [4] and the study of equations of Kolmogorov type as presented by Freedman and Waltman [5]. The studies cited above, together with lots of other contributions, lead to a variety of functional responses for different species which are widely used in prey-predator interaction models.Regarding discrete population models, we find it fair to say that there was a major breakthrough in 1976 when Sir Robert May [6] published his influential Nature paper where he showed that a simple one-dimensional nonlinear difference equation model could generate dynamics of stunning complexity, ranging from stable fixed points, periodic orbits of even and odd periods, and chaotic behaviour. Later, the number of papers on discrete population models flourished (confer [7–11]), and it became clear that the dynamics found from these studies was much richer than from their continuous counterparts. Ergodic properties of discrete models may be obtained in [12, 13] while the question of permanence is addressed in [14]. Discrete harvest models, both with or without age structure, are studied in [15–17].Parallel to the development of discrete age and stage structured population models, it became also customary to analyze prey-predator models formulated in discrete time. Indeed, Neubert and Kot [18] showed that the equilibrium in a two species prey-predator model may undergo a subcritical flip bifurcation with a subsequent concomitant crash of the predator population. Other excellent studies may be obtained from [19–24] and, more recently, the dynamical behaviour of fractional order Lotka-Volterra and generalized Lotka-Volterra models together with its discretizations have been scrutinized in [25, 26].Unlike most of the papers quoted above, we shall in this paper assume interactions between the prey and predator species of exponential form, a choice which is inspired by the seminal work of Ricker [27], also cf. [21]. The purpose of this work is to analyze (A) a two species prey-predator model and (B) two versions of a three species models (two predators), where focus is on stability, nonstationary, and chaotic behaviour as well as on mechanisms which may lead to extinction of predators. Regarding the results, we prove in case (A) that the size of the region in parameter space where the equilibrium is stable strongly depends on the fecundity of the predator and moreover that the transfer from stability to instability may go through either a supercritical flip bifurcation or alternatively through a supercritical Neimark-Sacker bifurcation when the fecundity of the prey is increased. In the chaotic regime there may be two different attractors, one where both the prey and the predator coexist and another where only the prey survives. We investigate the size of basin of attraction for these possibilities. In case (B) focus is much on the same as in (A). One major result is that if the top predator preys upon both the prey and the predator or only on the predator, this has profound effects on the size of stable parameter regions and on possible nonstationary dynamics.The plan of the paper is as follows. In Section2 we formulate and analyze the two species prey-predator model. Section 3 deals with the case where there is one prey population and two predator populations. Finally, in Section 4 we summarize and discuss results.


## 2. The 2-Dimensional Model

Letxt and yt be the sizes of a prey and a predator population at time t, respectively. The relation between the two species at two consecutive time steps (years) is assumed to be on the form(1)xt+1=fxt,ytxtyt+1=gxt,ytytNatural restrictions to impose are(2)∂f∂x≤0,∂f∂y≤0,∂g∂x≥0,∂g∂y≤0which biologically means that intraspecific competition leads to a decrease in size of both populations while interspecific competition (predation) results in a decrease of the survival of the prey and an increase of the size of the predator population.In this section we shall consider the model (which satisfies the restrictions above)(3)xt+1=Fe-xte-αytxtyt+1=Ge-yt1-e-βxtytThe capital letters F and G denote density independent fecundity terms. α and β are positive interaction parameters and from a biological point of view it is natural to assume α≥β. When β→0+, (3) degenerates to a ‘pure’ prey map.Map (3) possesses three equilibria, the trivial one (x¯,y¯)=(0,0), the point (x~,y~)=(lnF,0) where F>1, and the nontrivial one(4)x∗,y∗=x∗,1αln⁡F-x∗where x∗ satisfies the equation(5)e1/αx∗-e1/α-βx∗=F1/αG-1and ln⁡F>x∗. In order to investigate stability properties we linearize about the equilibrium. This gives birth to the eigenvalue equation(6)λ2+a1λ+a2=0where the coefficients are(7)a1=x∗+y∗-2,a2=1-x∗-y∗+Q∗+1x∗y∗,Q∗=αβF-1/αGe1/α-βx∗(x∗,y∗) is a stable equilibrium as long as all eigenvalues of (6) are located on the inside of the unit circle and according to the Jury criteria this is satisfied whenever the inequalities 1+a1+a2>0, 1-a1+a2>0, and 1-|a2|>0 hold. Following Murray [28], when the first of these inequalities fails (i.e., when 1+a1+a2=0), it corresponds to λ=1. The second one fails when λ=-1 (the flip case) and the third fails when the solution of (6) is a pair of complex valued eigenvalues located on the boundary of the unit circle (i.e., λ=e±iθ, the Neimark-Sacker case). Consequently, (x∗,y∗) is stable for those parameter combinations who satisfy(8a)Q∗+1x∗y∗>0(8b)x∗y∗>2x∗+y∗-4Q∗+1(8c)x∗y∗<x∗+y∗Q∗+1Hence, (4) will be stable whenever(9)f1<x∗y∗<f2where(10)f1=2x∗+y∗-4Q∗+1,f2=x∗+y∗Q∗+1and in order for f1<f2 we must have x∗+y∗<4.Iff1=x∗y∗ the solutions of (6) are λ1=-1 and λ2=3-(x∗+y∗). If x∗y∗=f2 is assumed, the solutions of (6) are(11)λ1,2=122-x∗-y∗±i4x∗+y∗-x∗+y∗2and these eigenvalues are indeed located on the boundary of the unit circle since(12)λ1,2=122-x∗+y∗2+4x∗+y∗-x∗+y∗22=1Note that when x∗+y∗→4, all eigenvalues (both real and complex) approach λ=-1 and the stable parameter region approaches zero.Example 1.
First we scrutinize the special caseα=β=1, i.e.,(13)xt+1=Fe-xte-ytxtyt+1=Ge-yt1-e-xtytBy use of the Jury criteria, it is straightforward to show that (x¯,y¯)=(0,0) is stable for 0<F<1 and all values of G. Regarding (x~,y~)=(ln⁡F,0) it is stable whenever 1<F<e2 and G<F(F-1)-1. Note that when F→1+, then G becomes arbitrary and when F→e2, then G→e2(e2-1)-1. At threshold G=F(F-1)-1, λ1=1 and λ2=1-ln⁡F which satisfies -1<λ2<1. The nontrivial equilibrium point may be expressed as(14)x∗,y∗=ln⁡F+GG,ln⁡FGF+GNote that x∗+y∗=ln⁡F and in order to have a feasible equilibrium we must assume FG>F+G  (or F>G(G-1)-1).The coefficients in (6) become(15)a1=ln⁡F-2,a2=1-ln⁡F+GF+1x∗y∗and subsequently the Jury criteria (8a)-(8c) may be cast in the form(16a)F+GFx∗y∗>0(16b)x∗y∗>2Fln⁡F-2F+G=f1F,G(16c)x∗y∗<Fln⁡FF+G=f2F,GThe stable parameter region is characterized by f1<x∗y∗<f2. Moreover, (16b) and f1<f2 imply e2<F<e4. Hence, depending on G, the largest F interval where (x∗,y∗) may be stable is e2<F<e4.At thresholdf1=x∗y∗ the equilibrium point may be expressed as(17)x∗,y∗=b,c=12ln⁡F+a,12ln⁡F-awhere(18)a=12ln⁡F2-2Fln⁡F-2F+GThe solutions of (6) are λ1=-1 and λ2=3-ln⁡F. Note that when F→e2, then (b,c)→(2,0) and G=e2(e2-1)-1≈1.157. When F→e4, the quantities b, c, and G may be expressed by use of the Lambert function W0 as (b,c)=(2(1+W0(e-1)),2(1-W0(e-1)))≈(2.55,1.44) and G=-e4(W0(e-1))2((W0(e-1))2-1)-1≈4.58956.Fromx∗y∗=f2 we find(19)x∗.y∗=p,q=12ln⁡F+d,12ln⁡F-dwhere(20)d=12ln⁡F2-Fln⁡FF+Gand the corresponding complex modulus 1 eigenvalues become(21)λ=2-ln⁡F2±iln⁡F4-ln⁡F4In Figure 1(a) we have visualized the stable parameter regions in the F-G plane. When 0<F<1, (x¯.y¯)=(0,0) is stable for any value of G. In the interval 1<F<e2, the point (x~,y~)=(ln⁡F,0) is stable for G<F(F-1)-1, and whenever e2<F<e4, the nontrivial equilibrium (x∗.y∗) given by (14) is stable for those combinations of F and G which satisfy f1<x∗y∗<f2 (confer (16b), (16c)); i.e., the stable parameter region is located between the curves.Figure 1
(a) Stability regions for equilibria:F∈〈0,1〉: (x¯,y¯)=(0,0) is stable for any G. F∈〈1,e2〉: (x~,y~)=(ln⁡F,0) is stable under the curve G=F(F-1)-1. F∈〈e2,e4〉: (x∗,y∗)=(ln⁡F+GG-1,ln⁡FG(F+G)-1) is stable between the curves. (b) The points (F,G)=(41.438,4) and (F,G)=(40.18,5) in parameter space are located on the flip and Neimark-Sacker bifurcation curve, respectively.
(a) (b)Our next goal is to study the nonstationary dynamics as(x∗,y∗) becomes unstable. This is done by considering two representative fixed values of G, G=4 and G=5, see Figure 1(b), and vary F. Assuming G=4, (x∗,y∗) is stable as long as e2<F<41.438 and undergoes a flip bifurcation when F is increased to F=41.438 where f1=x∗y∗. G=5 implies that (x∗,y∗) is stable in the interval e2<F<40.18 and undergoes a Neimark-Sacker bifurcation when F becomes 40.18 where x∗y∗=f2.As it is well known, a flip bifurcation may be of both supercritical or subcritical nature. In the former case, there exists a stable 2-cycle just beyond instability threshold. In the latter case such a stable 2-cycle does not exist. Regarding our model we have the following result:Theorem 2.
Consider the map(22)x,y→Fe-xe-yx,Ge-y1-e-xyunder the assumptions F∈(e2,e4) and G∈(1.157,4.59). Then for all combinations of F and G such that x∗y∗=2F(ln⁡F-2)(F+G)-1, the equilibrium (x∗,y∗) undergoes a supercritical flip bifurcation.Proof.
S e e  Appendix.Figure2(a) displays a stable 2-cycle just beyond threshold in the case (F,G)=(45,4). The figure clearly demonstrates how the orbits of the prey and predator are synchronized. The amplitude of the predator follows one time unit (year) after the amplitude of the prey.Figure 2
(a) 2-Periodic dynamics generated by (13). Parameter values (F,G)=(45,4). (b) Values of Lyapunov exponent L of orbits generated by (13). G=4 and 30<F<60. (c) Bifurcation diagram generated by (13). Same parameter values as in (b). (d) Chaotic dynamics generated by (13). (F,G)=(60,4). (e) x and y as functions of time. (F,G)=(60,4). (f) The route to attractor A2. (F,G)=(65,4).
(a) (b) (c) (d) (e) (f)Still assumingG=4, we may describe the dynamics of the system for a range of F values by use of the Lyapunov exponent L of the orbit generated by (13). In Figure 2(b) we have computed L for F values between 30 and 60. L<0 in the interval e2<F<41.438 where (x∗,y∗) is stable. Stable 2-cycles (L<0) exist whenever 41.438<F<53.454 which is followed by stable 4-cycles when 53.454<F<56.173. The findings above are also displayed in the bifurcation diagram, Figure 2(c). Chaos is introduced as L becomes positive, i.e., when F exceeds 57.045. In Figure 2(d) we visualize the dynamics for the pair (F,G)=(60,4). Although the dynamics occurs in the chaotic regime, one may still argue that a certain kind of two-periodicity is preserved. This is due to the fact that each of the two subsets of the attractor shown in Figure 2(d) are visited only once every second iteration. Finally, still considering (F,G)=(60,4), Figure 2(e) shows x and y as functions of time. Evidently, the synchronization found in the 2-periodic case (Figure 2(a)) is necessarily not present in the chaotic regime.There is also another kind of dynamics that may occur. It dominates completely whenF exceeds 63.61 and is visualized in Figure 2(f) where F=65. For several iterations, the populations exhibit chaotic oscillations similar to what is shown in Figures 2(d) and 2(e) where F=60, but once x falls below a critical value xc we observe a dramatic change. When x<xc, the predator population becomes very small as well, actually so small that it does not manage to recover and consequently goes extinct. At the same time, in case of x small, map (13) degenerates to xt+1=Fxt. Hence, the prey indeed manages to recover and may in fact be large before it is damped again by the factor e-x. This mechanism explains the change of dynamics seen in Figure 2(f). Thus for F>63.61 the only possibility is an attractor, which we from now on shall refer to as A2, where the prey shows highly oscillatory behaviour and y=0. For all practical purposes we may regard A2 as generated by the pure prey map x→Fexp⁡(-x)x (after the original map (13) first has driven the predator to extinction). A2 also exists for F≤63.61. There A2 coexists with the attractors already accounted for. In Figure 3 we show the situation in somewhat more detail. For each F value in Figure 3 we have considered about 20000 different initial values (x0,y0) and for each of them performed 10000 iterations (map (13)) to see where the corresponding orbit settles. The fraction of all orbits starting at (x0,y0) which does not converge towards A2 is below the curve. Hence, for F<36 (roughly), all orbits converge towards (x∗,y∗); thus, (x∗,y∗) is not locally but also globally stable in this region. As we continue to increase F, 36<F<63.61 there is coexistence between the stable equilibrium (x∗,y∗) and A2, periodic orbits and A2, and chaotic orbits and A2, and we observe that the basin of attraction for A2 gradually increases as F becomes larger. The ultimate situation occurs when F=63.61; then, all orbits converge towards A2.Figure 3
Below the curve is the fraction of all orbits starting at20000 different initial conditions (x0,y0) which do not converge towards attractor A2. G=4, 10<F<80.Next, considerG>4.59 and e2<F<e4. If an equilibrium shall undergo a supercritical Neimark-Sacker bifurcation, then a pair of complex valued eigenvalues must cross the unit circle outwards at instability threshold and an attracting quasiperiodic orbit restricted to an invariant curve is created beyond the threshold. Regarding (13) we have the following:Theorem 3.
Consider the map(23)x,y→Fe-xe-yx,Ge-y1-e-xyunder the assumptions F∈(e2,e4) and G>4.59. Then, for those combinations of F and G such that x∗y∗=Fln⁡F(F+G)-1, the equilibrium undergoes a supercritical Neimark-Sacker bifurcation.Proof.
See Appendix.AssumingG=5, a supercritical Neimark-Sacker bifurcation will occur when F is increased to 40.18. Figure 4(a), where (F,G)=(43,5), shows an attracting invariant curve as predicted by Theorem 3. There is also a clear tendency that the predator population has a peak one unit of time later than the prey, but the tendency is not as profound as in cases where G<4.59 (the ‘flip’ cases). This is exemplified in Figure 4(b).Figure 4
(a) An invariant curve just beyond Neimark-Sacker bifurcation threshold. Parameter values(F,G)=(43,5). (b) x and y as functions of time in the case (F,G)=(43,5).
(a) (b)Still keepingG=5 fixed, we have also in this case computed the Lyapunov exponent L and the bifurcation diagram for different values of F. Results are presented in Figure 5 where we have used the initial value (x0,y0)=(2.5,1.25). The equilibrium (x∗,y∗) is stable when e2<F<40.18 and the invariant curve which is created at F=40.18 persists in the interval 40.18<F<e4. For higher values of F we find windows (L<0) where the dynamics is periodic as well as new invariant curves. The rationale behind this dynamics is as follows: once the invariant curve is established, map (13) is basically topological equivalent to a circle map. Associated with a circle map is a rotation number σ and whenever there are specific F values such that σ becomes rational, the dynamics is periodic. Moreover, the implicit function theorem guarantees that the periodicity is maintained in an interval about the specific F values as well. Chaos is introduced when L becomes positive, i.e., when F exceeds 65.6, which is interrupted by windows of different widths where the dynamics is periodic (L<0). The large window corresponds to period 2 orbits.Figure 5
(a) The Lyapunov exponentL of orbits generated by (13). G=5, 30<F<100. Initial values (x0,y0)=(2.5,1.25). (b) Bifurcation diagram generated by (13). Same parameter values as in (a).
(a) (b)The value ofL jumps when F reaches 90.18. Depending on (x0,y0) the ultimate fate of orbits generated by (13) is now A2. This is exemplified in Figure 6(a) when (F,G)=(91,5) and (x0,y0)=(2.5,1.25). Those initial conditions and associated orbits which do not settle on A2 end up on a chaotic attractor as displayed in Figure 6(b). This attractor has been established through a series of period doubling bifurcations as a result of increasing F from the ‘2 period window’ in Figure 5. When F exceeds 94.6, all orbits are attracted by A2.Figure 6
(a) Convergence towards attractorA2. (F,G)=(91,5). Initial condition (x0,y0)=(2.5,1.25). (b) Chaotic dynamics where both species survive. (F,G) as in Figure 6(a) but (x0,y0)=(2.5,1.7).
(a) (b)Example 4.
Next, consider the caseα=1 and β=1/2. This means that we are now entering a part of a parameter space where the predator gains less from the prey than the prey is able to ‘give’, so in many respects one may argue that this example is of more biological relevance than the previous one.The nontrivial equilibrium point is found to be(24)x∗,y∗=ln⁡1+a22,ln⁡4F1+a2where a=(4F+G)G-1. Moreover, note that the difference between x∗ given by (24) and (14) may be written as(25)x24∗-x14∗=ln⁡1+1+a-12F+G>0and further(26)y24∗-y14∗=ln⁡2FF+G2FF+G+FGa-1<0Hence, when β<α we conclude that the size of x∗ increases and the size of y∗ decreases compared to the case α=β. Biologically, this makes sense. β<α means that the predator will benefit less by eating which in turn will lead to a decrease of the size of the predator population. Subsequently, the predation pressure on the prey becomes smaller which again leads to an increase of the prey population.At bifurcation thresholdf1=x∗y∗ (confer (9)) the eigenvalues are λ1=-1 and λ2=3-(x∗+y∗)=3-ln⁡F, while the eigenvalues at the Neimark-Sacker threshold x∗y∗=f2 are given by (11) and (24). Consequently, we do not expect any qualitative changes of the dynamics compared to the symmetric case α=β=1.Numerically, we have found that whenG=4.623 the graphs of f1, x∗y∗ and f2 intersect at F=e4. This is displayed in Figure 7. Thus, the largest possible F interval (e2<F<e4) where (x∗,y∗) may be stable occurs for this particular value of G. Note that G here is slightly larger than in the α=β=1 case which suggests that a decrease of β acts in a stabilizing fashion.Figure 7
Graphs off1, x∗x∗, and f2, F<100. The F region where (x∗,y∗) is stable is located to the left of the vertical line. G=4.623.For comparison reasons we have also investigated the casesG=4 and G=5 in somewhat more detail. Regarding the former, we find that equilibrium (24) is stable in the interval e2<F<37. At F=37, (24) undergoes a (supercritical) flip bifurcation and a stable period 2 orbit is established. Through further enlargement of F periodic orbits of period 2k, k>1 is the dynamical outcome. Chaos is introduced when F=41 and the structure of the attractor is similar to the one shown in Figure 2(d). For higher values of F  (F>46), chaotic behaviour of A2 form dominates. (Note that some of the F thresholds above strongly depend on the initial conditions.) Thus, compared with the symmetric case α=β=1, see the Lyapunov exponent diagram in Figure 2(b); the case under consideration appears to be somewhat more unstable in the sense that the road from stability to chaos appears to be shorter here.Considering the caseG=5, we find that e2<F<39 guarantees a stable equilibrium. At F=39 a (supercritical) Neimark-Sacker bifurcation takes place, creating an attracting invariant curve which persists until F=62 where it breaks up and chaotic oscillations are the outcome. Chaotic attractors of A2 form have also been identified just as in the α=β=1 case.


## 3. The 3-Dimensional Models

In this section we shall extend the two-dimensional model (3) by including one more predator which we shall refer to as the top predator z. There is basically two ways of including z: one way is to assume that z preys upon y only, and the other possibility is to assume that z preys upon both x and y. In order to analyze the former case we consider the model(27)xt+1=Fe-xte-αytxtyt+1=Ge-yt1-e-βxte-γztytzt+1=He-zt1-e-δytztIn the latter case:(28)xt+1=Fe-xte-αyte-βztxtyt+1=Ge-yt1-e-γxte-δztytzt+1=He-zt1-e-ϵxt1-e-σytztModel (27) possesses the nontrivial equilibrium(29)x∗,y∗,z∗=ln⁡Fe-αy∗,y∗,ln⁡H1-e-δy∗where y∗ satisfies the equation(30)eαβ-1y∗1-Fβe-αβy∗1-e-δγy∗=-FβHγGRegarding equilibrium (x∗,y∗,z∗) of model (28) we find(31)x∗,y∗,z∗=ln⁡Fe-αy∗e-βz∗,y∗,z∗and y∗, z∗ must be obtained from the system(32)He-z∗Fϵe-αϵy∗e-βϵz∗-11-e-σy∗-Fϵe-αϵy∗e-βϵz∗=0Ge-y∗Fγe-αγy∗e-βγz∗-1e-δz∗-Fγe-αγy∗e-βγz∗=0Stability analysis is performed by linearizing about the equilibrium in the same way as in Section 2. This gives birth to eigenvalue equations of form(33)λ3+a1λ2+a2λ+a3=0where the coefficients a1-a3 in case of model (27) and equilibrium (29) become(34)a1=x∗+y∗+z∗-3a2=3-2x∗+y∗+z∗+1+αβGe-y∗-1x∗y∗+x∗z∗+1+γδHe-z∗-1y∗z∗a3=x∗+y∗+z∗-1+αβGe-y∗-1x∗y∗-x∗z∗-1+γδHe-z∗-1y∗z∗+γδHe-z∗-1+αβGe-y∗-1x∗y∗z∗-1while the corresponding coefficients considering model (28) and equilibrium (31) may be written as(35)a1=x∗+y∗+z∗-3a2=3-2x∗+y∗+z∗+1+αγAx∗y∗+1+βϵBx∗z∗+1+δσCy∗z∗a3=x∗+y∗+z∗-1+αγAx∗y∗-1+βϵBx∗z∗-1+δσCy∗z∗+1+δσC+αγA-αδϵB+βσγAC+βϵBx∗y∗z∗-1Here(36)A=Ge-y∗e-δz∗-1B=He-z∗1-e-σy∗-1C=He-z∗1-e-ϵx∗-1Equilibria (29) and (31) are stable as long as all eigenvalues of (33) are located inside the unit circle in the complex plane. According to the Jury criteria [28] this is satisfied as long as the inequalities(37)S1=1+a1+a2+a3>0S2=1-a1+a2-a3>0S3=1-a3>0S4=1-a32-a2-a3a1>0hold.The strategy we shall use in order to study the impact from top predatorz on x and y is to fix values of fecundities F and G and then study the dynamic consequences of varying H.Example 5.
We start by considering model (27). The parameter space is huge so at first we limit the discussion and assume that α=β=γ=δ=1. Then, the nontrivial equilibrium (29) may be expressed as(38)x∗,y∗,z∗=ln⁡HF+GG+H,ln⁡FG+HHF+G,ln⁡HGF-1FG+Hand in order for z∗>0, H must satisfy the inequality(39)H>FGFG-F-GMoreover, left hand side of S1 (see (37)) reduces to(40)H+GF-1e-z∗x∗y∗z∗which is clearly positive for all nonzero equilibrium populations x∗, y∗, and z∗.In order to investigate the impact from the top predatorz on the dynamics in more detail we consider the case (F,G)=(43,5) (which means that for z∗=0, the dynamics occurs on an invariant curve in the x-y plane as already accounted for; confer Figure 4(a)). When z∗ (or H) is small, S2>0, S3<0, S4>0 and we find quasiperiodic behaviour restricted on invariant curves in the x-y-z plane; see Figure 8(a). If we continue to increase H, the nontrivial fixed point (38) becomes attracting; confer Figure 8(b). Thus, in this part of parameter space an enlargement of H acts stabilizing. However, this region is small. Indeed, through further increase of H, S2 becomes negative and stable periodic orbits of period 2k, k≥1 are established; see Figure 8(c). Eventually, the system exhibits chaotic behaviour (see Figure 8(d)), and when H exceeds 2.46, both predators are driven to extinction and the prey performs chaotic oscillations similar to what was reported in Section 2. The graphs of S2, S3, and S4 may be obtained in Figure 9. Turning to the parameter combination (F,G)=(45,4) we know from the analysis in Section 2 that the dynamics is a stable period 2 orbit (confer Figure 2(a)). For small values of z∗, S2<0 and becomes even more negative as H grows. Hence, we observe qualitatively the same dynamical behaviour as in the (F,G)=(43,5) case. The only difference really is that it is not possible to obtain an interval where (x∗,y∗,z∗) is stable.Figure 8
(a) An invariant curve in3-space in case of H small. Parameter values (F,G,H)=(43,5,1.3). (b) Convergence towards the stable equilibrium (x∗,y∗,z∗). Parameter values (F,G,H)=(43,5,1.32). (c) A 2 cycle generated by model (27). Parameter values (F,G,H)=(43,5,1.64). (d) Chaotic dynamics. Parameter values (F,G,H)=(43,5,2.38).
(a) (b) (c) (d)Figure 9
Graphs of Jury criteriaS2, S3, and S4, cf. (34) and (37).Example 6.
Next, suppose thatz preys upon both x and y. Then, under the assumption α=β=γ=δ=ϵ=σ=1 equilibrium (31) of model (28) is found to be(41)x∗,y∗,z∗=ln⁡F+GG,ln⁡G+HH,ln⁡FGHF+GG+Hand(42)H>GF+GFG-F-Gin order to ensure z∗>0. The equilibrium is stable as long as the Jury criteria (37) are satisfied.In order to account for the dynamics we refer to Figure10, where we have computed the Lyapunov exponent L of the orbit generated by (28) in the case (F,G)=(43,5) and 1.43<H<9. For H close to 1.4371, L=0 and we observe quasiperiodic behaviour. When H exceeds 1.5, equilibrium (41) becomes stable and remains stable until H=3.98 where S2 (confer (37)) equals zero and (41) undergoes a flip bifurcation. Note that the stable interval here is much larger than the corresponding interval found in Example 5. This finding suggests that if the top predator preys upon both populations, it leads to more stable dynamics compared to the situation where z preys upon y only. Whenever 3.98<H<6.10, there are stable period 2 orbits and at H=6.10 a stable period 4 orbit is established through yet another flip. However, in contrast to all cases discussed earlier, we do not experience the flip bifurcation sequence and orbits of period 2k, k≥3 as H is further increased. Instead at H=7.3 the fourth iterate of (28) undergoes a Neimark-Sacker bifurcation and the outcome is four disjoint invariant curves that are visited once every fourth iteration; see Figure 11(a). L becomes positive at H=7.9 and the dynamics turns chaotic, as displayed in Figure 11(b). After the jump at H=8.2 there are chaotic oscillations where y=z=0, i.e., the same kind of dynamics as found in Example 5.Figure 10
Values of Lyapunov exponentL of orbits generated by (28). Parameters: (F,G)=(43,5) and 1.43<H<9.Figure 11
(a)4 invariant curves just beyond threshold where the fourth iterate of (28) has gone through a Neimark-Sacker bifurcation. Parameter values (F,G,H)=(43,5,7.5). (b) Chaotic dynamics generated by (28). Parameter values (F,G,H)=(43,5,8).
(a) (b)


## 4. Discussion

In this paper we have analyzed two and three species prey-predator models. Focus has been on stability properties and dynamic behaviour in unstable and chaotic parameter regions. As accounted for, the parameter space is huge. Therefore, we have mainly concentrated on special cases, for example, the casesα=β=1 or α=1, β=1/2 in the two species model. Regarding fecundity terms (F, G, and H), results from only selected values of G have been presented. However, these values indeed seem to be representative as lots of numerical experiments with other values suggest.First, let us comment on the two species models. The largestF interval where (x∗,y∗) may be stable is e2<F<e4. This interval exists for only one value of G, G=GC (GC=4.59 when α=β=1, GC=4.623 if α=1, β=1/2). When G<GC, the region where (x∗,y∗) is stable is on the form e2<F<C1, while G>GC implies the region e2<F<C2 and C1,C2<e4. Moreover, for a fixed value of G,G<GC the equilibrium will undergo a supercritical flip bifurcation at threshold F=C1 while a supercritical Neimark-Sacker bifurcation occurs at F=C2,G>GC. Thus, the value of G decides what kind of nonstationary behaviour model (3) may generate as F is increased. In the former case an increase of G acts stabilizing while it acts in a destabilizing way in the latter case. For a given value of G an enlargement of F acts in a destabilizing fashion and when F becomes sufficiently large, chaotic behaviour is the outcome. However, depending G<GC or G>GC, the routes to chaos will be different. Scrutinizing the chaotic regime, there are two possibilities: (A) there may be chaotic oscillations where both populations survive or (B) the prey exhibits chaotic large amplitude oscillations while the predator faces extinction. If F exceeds a critical threshold FC, possibility (B) will always be the case. If F<FC and |F-FC| is small, the ultimate fate of an orbit (extinction of the predator or not) strongly depends on the initial conditions (x0,y0).Turning to the three species models (27), (28) our finding is that the inclusion of a top predator z may act qualitatively different on the dynamics depending on whether y or both x and y are targets. If z preys upon y only there exists a tiny region in state space where z is small, and if increased, acts in a stabilizing fashion. However, if z becomes larger, further increase has a strong destabilizing effect. The transfer from a stable equilibrium (x∗,y∗,z∗) to nonstationary behaviour occurs when an eigenvalue λ of (33) crosses the unit circle at -1 (S2=0) and periodic dynamics of period 2k, k≥1 is established. Eventually, the dynamics becomes chaotic. The reason why an increase of z acts destabilizing may be understood along the following line: when z preys upon y only, the size of y becomes smaller which reduces the impact on x from y. This allows x to perform large amplitude oscillations for smaller values of F than if z had been absent. Consequently, in a worst-case scenario, x may fall below a critical value xC such that y does not manage to recover from a low population density and the result is extinction of both y and z.Whenz preys upon both species x and y, we find changes regarding stability properties as well as changes concerning the nonstationary dynamics. As long as the interaction parameters are equal, the region where (x∗,y∗,z∗) is stable appears to be significantly larger than in the previous case where z preys upon y only. Thus, one may argue that the inclusion of small and intermediate sized top predator populations acts stabilizing compared to the findings in the former case. The total predation pressure on x is on a level where it prevents large amplitude oscillations of the prey. However, if H is further increased, we observe just as in the former case, periodic dynamics of period 2 and 4 but not periodic dynamics of period 2k, k≥3. Instead, the fourth iterate of (28) undergoes a Neimark-Sacker bifurcation and 4 invariant closed curves are established which again implies that the route to chaos is different from the route when z preys upon y only. Consequently, there are large regions in parameter space where the dynamical outcomes are different between the two models (27), (28) under consideration. However, if we continue to increase H, we will eventually arrive at the same situation as we experienced through model (27); namely, that both the predator and the top predator will die.



---

*Source: 9543139-2019-02-17.xml*