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

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.


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][8][9][10][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][16][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][20][21][22][23][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 Section 2 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.

The 2-Dimensional Model
Let   and   be the sizes of a prey and a predator population at time , respectively.The relation between the two species at two consecutive time steps (years) is assumed to be on the form Natural restrictions to impose are which 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) The capital letters  and  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 (, ) = (0, 0), the point ( x, ỹ) = (ln, 0) where  > 1, and the nontrivial one where  * satisfies the equation and ln  >  * .In order to investigate stability properties we linearize about the equilibrium.This gives birth to the eigenvalue equation where the coefficients are ( * ,  * ) 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 +  1 +  2 > 0, 1 −  1 +  2 > 0, and 1 − | 2 | > 0 hold.Following Murray [28], when the first of these inequalities fails (i.e., when 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.,  =  ± , the Neimark-Sacker case).Consequently, ( * ,  * ) is stable for those parameter combinations who satisfy Hence, (4) will be stable whenever where and in order for  1 <  2 we must have  * +  * < 4.
If  1 =  *  * the solutions of ( 6) are  1 = −1 and  2 = 3 − ( * +  * ).If  *  * =  2 is assumed, the solutions of (6) are and these eigenvalues are indeed located on the boundary of the unit circle since Note that when  * +  * → 4, all eigenvalues (both real and complex) approach  = −1 and the stable parameter region approaches zero.
Example .First we scrutinize the special case  =  = 1, i.e., By use of the Jury criteria, it is straightforward to show that (, ) = (0, 0) is stable for 0 <  < 1 and all values of .
The coefficients in (6) become and subsequently the Jury criteria (8a)-(8c) may be cast in the form The stable parameter region is characterized by  1 <  *  * <  2 .Moreover, (16b) and  1 <  2 imply  2 <  <  4 .Hence, depending on , the largest  interval where ( * ,  * ) may be stable is  2 <  <  4 .At threshold  1 =  *  * the equilibrium point may be expressed as where The solutions of ( 6) are  1 = −1 and  2 = 3 − ln .Note that when  → where and the corresponding complex modulus 1 eigenvalues become In Figure 1(a) we have visualized the stable parameter regions in the  −  plane.When 0 <  < 1, (.) = (0, 0) is stable for any value of .In the interval 1 <  <  2 , the point ( x, ỹ) = (ln , 0) is stable for  < ( − 1) −1 , and whenever  2 <  <  4 , the nontrivial equilibrium ( * .* ) given by ( 14) is stable for those combinations of  and  which satisfy )); i.e., the stable parameter region is located between the curves.Our next goal is to study the nonstationary dynamics as ( * ,  * ) becomes unstable.This is done by considering two representative fixed values of ,  = 4 and  = 5, see 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:
Proof.See Appendix.
Figure 2(a) displays a stable 2-cycle just beyond threshold in the case (, ) = (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.
Still assuming  = 4, we may describe the dynamics of the system for a range of  values by use of the Lyapunov exponent  of the orbit generated by (13).In Figure 2(b) we have computed  for  values between 30 and 60.  < 0 in the interval  2 <  < 41.438 where ( * ,  * ) is stable.Stable 2-cycles ( < 0) exist whenever 41.438 <  < 53.454 which is followed by stable 4-cycles when 53.454 <  < 56.173.The findings above are also displayed in the bifurcation diagram, Figure 2(c).Chaos is introduced as  becomes positive, i.e., when  exceeds 57.045.In Figure 2(d) we visualize the dynamics for the pair (, ) = (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 (, ) = (60, 4), Figure 2(e) shows  and  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 when  exceeds 63.61 and is visualized in Figure 2(f) where  = 65.For several iterations, the populations exhibit chaotic oscillations similar to what is shown in Figures 2(d) and 2(e) where  = 60, but once  falls below a critical value   we observe a dramatic change.When  <   , 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  small, map (13) degenerates to  +1 =   .Hence, the prey indeed manages to recover and may in fact be large before it is damped again by the factor  − .This mechanism explains the change of dynamics seen in Figure 2(f).Thus for  > 63.61 the only possibility is an attractor, which we from now on shall refer to as 2, where the prey shows highly oscillatory behaviour and  = 0.For all practical purposes we may regard 2 as generated by the pure prey map  →  exp(−) (after the original map (13) first has driven the predator to extinction).2 also exists for  ≤ 63.61.There 2 coexists with the attractors already accounted for.In Figure 3 we show the situation in somewhat more detail.For each  value in Figure 3 we have considered about 20000 different initial values ( 0 ,  0 ) and for each of them performed 10000 iterations (map ( 13)) to see where the corresponding orbit settles.The fraction of all orbits starting at ( 0 ,  0 ) which does not converge towards 2 is below the curve.Hence, for  < 36 (roughly), all orbits converge towards ( * ,  * ); thus, ( * ,  * ) is not locally but also globally stable in this region.As we continue to increase , 36 <  < 63.61 there is coexistence between the stable equilibrium ( * ,  * ) and 2, periodic orbits and 2, and chaotic orbits and A2, and we observe that the basin of attraction for 2 gradually increases as  becomes larger.The ultimate situation occurs when  = 63.61;then, all orbits converge towards 2.
Next, consider  > 4.59 and  2 <  <  4 .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: under the assumptions  ∈ ( 2 ,  4 ) and  > 4.59.en, for those combinations of  and G such that  *  * =  ln ( + ) −1 , the equilibrium undergoes a supercritical Neimark-Sacker bifurcation.
Proof.See Appendix.
Assuming  = 5, a supercritical Neimark-Sacker bifurcation will occur when  is increased to 40.18. Figure 4(a), where (, ) = (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  < 4.59 (the 'flip' cases).This is exemplified in Figure 4(b).
Still keeping  = 5 fixed, we have also in this case computed the Lyapunov exponent  and the bifurcation diagram for different values of .Results are presented in Figure 5 where we have used the initial value ( 0 ,  0 ) = (2.5, 1.25).The equilibrium ( * ,  * ) is stable when  2 <  < 40.18 and the invariant curve which is created at  = 40.18persists in the interval 40.18 <  <  4 .For higher values of  we find windows ( < 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  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  values as well.Chaos is introduced when  becomes positive, i.e., when  exceeds 65.6, which is interrupted by windows of different widths where the dynamics is periodic ( < 0).The large window corresponds to period 2 orbits.
The value of  jumps when  reaches 90.18.Depending on ( 0 ,  0 ) the ultimate fate of orbits generated by ( 13) is now 2.This is exemplified in Figure 6(a) when (, ) = (91, 5) and ( 0 ,  0 ) = (2.5, 1.25).Those initial conditions and associated orbits which do not settle on 2 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  from the '2 period window' in Figure 5.When  exceeds 94.6, all orbits are attracted by 2.
Example .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 where  = (4 + ) −1 .Moreover, note that the difference between  * given by ( 24) and ( 14) may be written as Hence, when  <  we conclude that the size of  * increases and the size of  * 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 threshold  1 =  *  * (confer ( 9)) the eigenvalues are  1 = −1 and  2 = 3 − ( * +  * ) = 3 − ln , while the eigenvalues at the Neimark-Sacker threshold  *  * =  2 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 when  = 4.623 the graphs of  1 ,  *  * and  2 intersect at  =  4 .This is displayed in Figure 7. Thus, the largest possible  interval ( 2 <  <  4 ) where ( * ,  * ) may be stable occurs for this particular value of .Note that  here is slightly larger than in the  =  = 1 case which suggests that a decrease of  acts in a stabilizing fashion.
For comparison reasons we have also investigated the cases  = 4 and  = 5 in somewhat more detail.Regarding the former, we find that equilibrium (24) is stable in the interval  2 <  < 37.At  = 37, (24) undergoes a  (supercritical) flip bifurcation and a stable period 2 orbit is established.Through further enlargement of  periodic orbits of period 2  ,  > 1 is the dynamical outcome.Chaos is introduced when  = 41 and the structure of the attractor is similar to the one shown in Figure 2(d).For higher values of  ( > 46), chaotic behaviour of 2 form dominates.
(Note that some of the  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 case  = 5, we find that  2 <  < 39 guarantees a stable equilibrium.At  = 39 a (supercritical) Neimark-Sacker bifurcation takes place, creating an attracting invariant curve which persists until  = 62 where it breaks up and chaotic oscillations are the outcome.Chaotic attractors of 2 form have also been identified just as in the  =  = 1 case.

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 .There is basically two ways of including : one way is to assume that  preys upon  only, and the other possibility is to assume that  preys upon both  and .In order to analyze the former case we consider the model In the latter case: Model (27) possesses the nontrivial equilibrium where  * satisfies the equation Regarding equilibrium ( * ,  * ,  * ) of model ( 28) we find and  * ,  * must be obtained from the system Stability analysis is performed by linearizing about the equilibrium in the same way as in Section 2. This gives birth to eigenvalue equations of form where the coefficients  1 −  3 in case of model (27) and equilibrium (29) become while the corresponding coefficients considering model (28) and equilibrium (31) may be written as Here Equilibria ( 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 hold.
The strategy we shall use in order to study the impact from top predator  on  and  is to fix values of fecundities  and  and then study the dynamic consequences of varying .
Example .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 and in order for  * > 0,  must satisfy the inequality Moreover, left hand side of 1 (see (37)) reduces to which is clearly positive for all nonzero equilibrium populations  * ,  * , and  * .
In order to investigate the impact from the top predator  on the dynamics in more detail we consider the case (, ) = (43, 5) (which means that for  * = 0, the dynamics occurs on an invariant curve in the  −  plane as already accounted for; confer Figure 4(a)).When  * (or ) is small, 2 > 0, 3 < 0, 4 > 0 and we find quasiperiodic behaviour restricted on invariant curves in the  −  −  plane; see Figure 8(a).If we continue to increase , the nontrivial fixed point (38) becomes attracting; confer Figure 8(b).Thus, in this part of parameter space an enlargement of  acts stabilizing.However, this region is small.Indeed, through further increase of , 2 becomes negative and stable periodic orbits of period 2  ,  ≥ 1 are established; see Figure 8(c).Eventually, the system exhibits chaotic behaviour (see Figure 8(d)), and when  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 2, 3, and 4 may be obtained in Figure 9. Turning to the parameter combination (, ) = (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  * , 2 < 0 and becomes even more negative as  grows.Hence, we observe qualitatively the same dynamical behaviour as in the (, ) = (43, 5) case.The only difference really is that it is not possible to obtain an interval where ( * ,  * ,  * ) is stable.
Example .Next, suppose that  preys upon both  and .Then, under the assumption  =  =  =  =  =  = 1 equilibrium (31) of model ( 28) is found to be and in order to ensure  * > 0. The equilibrium is stable as long as the Jury criteria (37) are satisfied.
In order to account for the dynamics we refer to Figure 10, where we have computed the Lyapunov exponent  of the orbit generated by (28) in the case (, ) = (43, 5) and 1.43 <  < 9.For  close to 1.4371,  = 0 and we observe quasiperiodic behaviour.When  exceeds 1.5, equilibrium (41) becomes stable and remains stable until  = 3.98 where 2 (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  preys upon  only.Whenever 3.98 <  < 6.10, there are stable period 2 orbits and at  = 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 2  ,  ≥ 3 as  is further increased.Instead at  = 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). becomes positive at  = 7.9 and the dynamics turns chaotic, as displayed in Figure 11(b).After the jump at  = 8.2 there are chaotic oscillations where  =  = 0, i.e., the same kind of dynamics as found in Example 5.

Discussion
In this paper we have analyzed two and three species preypredator 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 (, , and ), results from only selected values of  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 largest  interval where ( * ,  * ) may be stable is  2 <  <  4 .This interval exists for only one value of ,  =   (  = 4.59 when  =  = 1,   = 4.623 if  = 1,  = 1/2).When  <   , the region where ( * ,  * ) is stable is on the form  2 <  <  1 , while  >   implies the region  2 <  <  2 and  1 ,  2 <  4 .Moreover, for a fixed value of ,  <   the equilibrium will undergo a supercritical flip bifurcation at threshold  =  1 while a supercritical Neimark-Sacker bifurcation occurs at  =  2 ,  >   .Thus, the value of  decides what kind of nonstationary behaviour model (3) may generate as  is increased.In the former case an increase of  acts stabilizing while it acts in a destabilizing way in the latter case.For a given value of  an enlargement of  acts in a destabilizing fashion and when  becomes sufficiently large, chaotic behaviour is the outcome.However, depending  <   or  >   , 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  exceeds a critical threshold   , possibility (B) will always be the case.If  <   and | −   | is small, the ultimate fate of an orbit (extinction of the predator or not) strongly depends on the initial conditions ( 0 ,  0 ).
Turning to the three species models ( 27), (28) our finding is that the inclusion of a top predator  may act qualitatively different on the dynamics depending on whether  or both  and  are targets.If  preys upon  only there exists a tiny region in state space where  is small, and if increased, acts in a stabilizing fashion.However, if  becomes larger, further increase has a strong destabilizing effect.The transfer from a stable equilibrium ( * ,  * ,  * ) to nonstationary behaviour occurs when an eigenvalue  of (33) crosses the unit circle at −1 (2 = 0) and periodic dynamics of period 2  ,  ≥ 1 is established.Eventually, the dynamics becomes chaotic.The reason why an increase of  acts destabilizing may be understood along the following line: when  preys upon  only, the size of  becomes smaller which reduces the impact on  from .This allows  to perform large amplitude oscillations for smaller values of  than if  had been absent.Consequently, in a worst-case scenario,  may fall below a critical value   such that  does not manage to recover from a low population density and the result is extinction of both  and .
When  preys upon both species  and , 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 ( * ,  * ,  * ) is stable appears to be significantly larger than in the previous case where  preys upon  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  is on a level where it prevents large amplitude oscillations of the prey.However, if  is further increased, we observe just as in the former case, periodic dynamics of period 2 and 4 but not periodic dynamics of period 2  ,  ≥ 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  preys upon  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 , 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.
Regarding the four terms in (A.12), the two terms in the middle are positive.When  →  2 , the third term dominates so  becomes very large.When  →  4 , the second term dominates and  becomes very large as well.For values of  not close to  2 or  4 , the second term is by far the largest and we conclude that (A.12) is positive for all parameter combinations.Thus, the bifurcation is supercritical.The functions ℎ and  contain second and third order terms of  and V and may be expressed as
At threshold  *  * =  2 (, ) =  ln ( + ) −1 the Jacobian becomes we define the transformation matrix  this time as in  are the eigenvectors corresponding to  = (2 − ln )/2 + (/2) of  and  = √ln (4 − ln).Then we proceed in exactly the same way as in the proof of Theorem 2, which allows us to cast the map (13) into standard form as