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.
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 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.
2. 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 formNatural restrictions to impose arewhich 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 , (3) degenerates to a ‘pure’ prey map.
Map (3) possesses three equilibria, the trivial one (, the point where , and the nontrivial onewhere satisfies the equationand . In order to investigate stability properties we linearize about the equilibrium. This gives birth to the eigenvalue equationwhere 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 , , and hold. Following Murray [28], when the first of these inequalities fails (i.e., when ), it corresponds to . The second one fails when (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 satisfyHence, (4) will be stable wheneverwhereand in order for we must have .
If the solutions of (6) are and . If is assumed, the solutions of (6) areand these eigenvalues are indeed located on the boundary of the unit circle sinceNote that when , all eigenvalues (both real and complex) approach and the stable parameter region approaches zero.
Example 1. First we scrutinize the special case , i.e.,By use of the Jury criteria, it is straightforward to show that is stable for and all values of . Regarding it is stable whenever and . Note that when , then becomes arbitrary and when , then . At threshold , and which satisfies . The nontrivial equilibrium point may be expressed asNote that and in order to have a feasible equilibrium we must assume (or ).
The coefficients in (6) becomeand subsequently the Jury criteria (8a)-(8c) may be cast in the formThe stable parameter region is characterized by . Moreover, (16b) and imply . Hence, depending on , the largest interval where may be stable is .
At threshold the equilibrium point may be expressed aswhereThe solutions of (6) are and . Note that when , then and . When , the quantities , , and may be expressed by use of the Lambert function as and .
From we findwhereand the corresponding complex modulus eigenvalues becomeIn Figure 1(a) we have visualized the stable parameter regions in the plane. When , is stable for any value of . In the interval , the point is stable for , and whenever , the nontrivial equilibrium given by (14) is stable for those combinations of and which satisfy (confer (16b), (16c)); i.e., the stable parameter region is located between the curves.
(a)
(b)
Our next goal is to study the nonstationary dynamics as becomes unstable. This is done by considering two representative fixed values of , and , see Figure 1(b), and vary . Assuming , is stable as long as and undergoes a flip bifurcation when is increased to where . implies that is stable in the interval and undergoes a Neimark-Sacker bifurcation when becomes where .
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 mapunder the assumptions and . Then for all combinations of and such that , the equilibrium undergoes a supercritical flip bifurcation.
Proof. Appendix.
Figure 2(a) displays a stable 2-cycle just beyond threshold in the case . 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.
(a)
(b)
(c)
(d)
(e)
(f)
Still assuming , 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 and . in the interval where is stable. Stable 2-cycles () exist whenever which is followed by stable 4-cycles when . The findings above are also displayed in the bifurcation diagram, Figure 2(c). Chaos is introduced as becomes positive, i.e., when exceeds . In Figure 2(d) we visualize the dynamics for the pair . 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 , Figure 2(e) shows and as functions of time. Evidently, the synchronization found in the -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 and is visualized in Figure 2(f) where . For several iterations, the populations exhibit chaotic oscillations similar to what is shown in Figures 2(d) and 2(e) where , 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 . 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 the only possibility is an attractor, which we from now on shall refer to as , where the prey shows highly oscillatory behaviour and . For all practical purposes we may regard as generated by the pure prey map (after the original map (13) first has driven the predator to extinction). also exists for . There 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 different initial values and for each of them performed iterations (map (13)) to see where the corresponding orbit settles. The fraction of all orbits starting at which does not converge towards is below the curve. Hence, for (roughly), all orbits converge towards ; thus, is not locally but also globally stable in this region. As we continue to increase , there is coexistence between the stable equilibrium and , periodic orbits and , and chaotic orbits and , and we observe that the basin of attraction for gradually increases as becomes larger. The ultimate situation occurs when ; then, all orbits converge towards .
Next, consider and . 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 mapunder the assumptions and . Then, for those combinations of and G such that , the equilibrium undergoes a supercritical Neimark-Sacker bifurcation.
Proof. See Appendix.
Assuming , a supercritical Neimark-Sacker bifurcation will occur when is increased to . Figure 4(a), where , 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 (the ‘flip’ cases). This is exemplified in Figure 4(b).
(a)
(b)
Still keeping 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 . The equilibrium is stable when and the invariant curve which is created at persists in the interval . For higher values of we find windows () 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 , which is interrupted by windows of different widths where the dynamics is periodic (). The large window corresponds to period orbits.
(a)
(b)
The value of jumps when reaches . Depending on the ultimate fate of orbits generated by (13) is now . This is exemplified in Figure 6(a) when and . Those initial conditions and associated orbits which do not settle on 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 , all orbits are attracted by .
(a)
(b)
Example 4. Next, consider the case and . 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 bewhere . Moreover, note that the difference between given by (24) and (14) may be written asand furtherHence, 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 (confer (9)) the eigenvalues are and , while the eigenvalues at the Neimark-Sacker threshold are given by (11) and (24). Consequently, we do not expect any qualitative changes of the dynamics compared to the symmetric case .
Numerically, we have found that when the graphs of , and intersect at . This is displayed in Figure 7. Thus, the largest possible interval where may be stable occurs for this particular value of . Note that here is slightly larger than in the case which suggests that a decrease of acts in a stabilizing fashion.
For comparison reasons we have also investigated the cases and in somewhat more detail. Regarding the former, we find that equilibrium (24) is stable in the interval . At , (24) undergoes a (supercritical) flip bifurcation and a stable period orbit is established. Through further enlargement of periodic orbits of period , is the dynamical outcome. Chaos is introduced when and the structure of the attractor is similar to the one shown in Figure 2(d). For higher values of (), chaotic behaviour of form dominates. (Note that some of the thresholds above strongly depend on the initial conditions.) Thus, compared with the symmetric case , 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 , we find that guarantees a stable equilibrium. At a (supercritical) Neimark-Sacker bifurcation takes place, creating an attracting invariant curve which persists until where it breaks up and chaotic oscillations are the outcome. Chaotic attractors of form have also been identified just as in the 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 . 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 modelIn the latter case:Model (27) possesses the nontrivial equilibriumwhere satisfies the equationRegarding equilibrium of model (28) we findand , must be obtained from the systemStability analysis is performed by linearizing about the equilibrium in the same way as in Section 2. This gives birth to eigenvalue equations of formwhere the coefficients in case of model (27) and equilibrium (29) becomewhile the corresponding coefficients considering model (28) and equilibrium (31) may be written asHereEquilibria (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 inequalitieshold.
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 5. We start by considering model (27). The parameter space is huge so at first we limit the discussion and assume that . Then, the nontrivial equilibrium (29) may be expressed asand in order for , must satisfy the inequalityMoreover, left hand side of (see (37)) reduces towhich 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 (which means that for , the dynamics occurs on an invariant curve in the plane as already accounted for; confer Figure 4(a)). When (or ) is small, , , 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 , becomes negative and stable periodic orbits of period , are established; see Figure 8(c). Eventually, the system exhibits chaotic behaviour (see Figure 8(d)), and when exceeds , both predators are driven to extinction and the prey performs chaotic oscillations similar to what was reported in Section 2. The graphs of , , and may be obtained in Figure 9. Turning to the parameter combination we know from the analysis in Section 2 that the dynamics is a stable period orbit (confer Figure 2(a)). For small values of , and becomes even more negative as grows. Hence, we observe qualitatively the same dynamical behaviour as in the case. The only difference really is that it is not possible to obtain an interval where is stable.
(a)
(b)
(c)
(d)
Example 6. Next, suppose that preys upon both and . Then, under the assumption equilibrium (31) of model (28) is found to beandin order to ensure . 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 and . For close to , and we observe quasiperiodic behaviour. When exceeds , equilibrium (41) becomes stable and remains stable until where (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 , there are stable period orbits and at a stable period 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 , as is further increased. Instead at 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 and the dynamics turns chaotic, as displayed in Figure 11(b). After the jump at there are chaotic oscillations where , i.e., the same kind of dynamics as found in Example 5.
(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 or , 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 . This interval exists for only one value of , ( when , if , ). When , the region where is stable is on the form , while implies the region and . Moreover, for a fixed value of the equilibrium will undergo a supercritical flip bifurcation at threshold while a supercritical Neimark-Sacker bifurcation occurs at . 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 .
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 () and periodic dynamics of period , 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 and but not periodic dynamics of period , . Instead, the fourth iterate of (28) undergoes a Neimark-Sacker bifurcation and 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.
Appendix
A. Proof of Theorem 2 in the Main Text
At threshold the Jacobian may be expressed asNext, define the matrixwhere the columns are the eigenvectors corresponding to and of , respectively. Then after expanding and up to third order, applying the change of coordinates (in order to translate the bifurcation to the origin), together with the transformationwe may cast map (13) into standard form aswhereandThe next step involves the restriction of (A.4) to the center manifold. To do this we first seek (approximate) the center manifold as a graphand by inserting (A.7) into the second component of (A.4) we arrive atfrom which we obtain and . Finally, by inserting (A.7) into the first component of (A.4), neglecting terms of order and higher, this results in the restricted mapwhereNow, appealing to Theorem 3.5 in [29] the bifurcation will be of supercritical nature whenever the quantity (evaluated at threshold)which is equivalent toRegarding the four terms in (A.12), the two terms in the middle are positive. When , the third term dominates so becomes very large. When , the second term dominates and becomes very large as well. For values of not close to or , 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.
B. Proof of Theorem 3 in the Main Text
At threshold the Jacobian becomesand we define the transformation matrix this time aswhere the columns in are the eigenvectors corresponding to of and . 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 asThe functions and contain second and third order terms of and and may be expressed asNow, following Wan [30], the Neimark-Sacker bifurcation will be of supercritical nature ifis negative at instability threshold and that (which ensures that the eigenvalues leave the unit circle at threshold). The quantities in areNext, assuming , we have computed all derivatives in (B.6) and the results of are then subsequently substituted back into (B.5). Then has been computed numerically in the interval . The graph of is displayed in Figure 12, and clearly . Finally, (at threshold)( in the special cases or ). Hence, we conclude that the bifurcation is supercritical.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The publication charges for this article have been covered by a grant from the publication fund of UiT, The Arctic University of Norway.