Compensatory and overcompensatory dynamics in prey–predator systems exposed to harvest

Density dependent prey–predator systems under the impact of harvest are considered. The recruitment functions for both the prey and predator belong to the Deriso–Schnute family which allow us to study how the dynamical behaviour of both populations changes when compensatory density dependence turns overcompensatory. Depending on the degree of overcompensation, we show in the case of no harvest that an increase of the fecundity of the prey always acts in a destabilizing fashion. If the degree of overcompensation becomes sufficiently large, such an increase can lead to large amplitude chaotic oscillations of the prey, which actually may drive the predator population to extinction. The impact of harvest also depends on the degree of overcompensatory density dependence. If only the prey is the target population, increased harvest in general seems to stabilize the dynamics. On the other hand, harvesting only the predator may in some cases tend to stabilize dynamics, but there are also parameter regions where this turns out to be a strong destabilizing effect.


Introduction
As summarized in [1], simple discrete population models in the form where x is population size, are perfectly capable of generating rich dynamical behaviour, ranging from stable equilibria, periodic orbits to chaotic oscillations. If f (x t ) is a declining function, the model is called density dependent, and if f (x t )x t eventually declines, the model is said to incorporate overcompensatory density depen-B Ørjan Kristensen orjan.f.kristensen@uit.no Arild Wikan arild.wikan@uit.no dence, as accounted for in [2]. If f (x t )x t lacks this property, the system is said to incorporate compensatory density dependence. Examples of population models that focus on dynamical consequences by use of overcompensatory density dependent recruitment functions, both linked to concrete species as well as pure theoretical studies, may be found in [3][4][5][6][7][8][9]. In the case of studies which incorporate compensatory density dependence, we refer to [10][11][12][13][14].
One-population models as described above may easily be extended. Indeed, prey-predator interactions may be scrutinized through models of the form x t+1 = f (x t , y t )x t , y t+1 = g(x t , y t )y t , where x is the prey and y the predator population. One common feature of several of these prey-predator models is that they possess three equilibria that correspond to (i) extinction of both species, (ii) survival of the prey and extinction of the predator, and (iii) survival of both species. Moreover, the interplay between the species may lead to dynamics of stunning complexity. There is a vast literature regarding such models and models which also include age or stage structure, see [15][16][17][18][19][20][21][22][23][24][25], and more recently, cf [26][27][28][29] and references therein, the question of how to control and stabilize chaotic behaviour in prey-predator systems has been addressed. The dynamics found in such systems strongly depends on compensatory or overcompensatory recruitment.
In cases where the prey and/or the predator populations are of commercial interest, harvest plays a crucial role with respect to the dynamical outcome and possible extinction of one or both species. Throughout the history there have been several reports of crisis due to low availability of economically important fish species. Today, the economics in the fishery industry is characterized by high fixed costs and high operational leverage, see [30], and is consequently vulnerable for severe variations in the fisheries regulations. Hence, one may argue that it is almost mandatory at some stage to include harvest in modelling such species. Excellent studies which address the role of harvest may be obtained in [31][32][33][34][35][36][37][38][39][40]. Now, following [41], there are in principle four factors which may explain oscillatory behaviour in (fish) populations. They are (A) environmental changes, (B) ecosystem dynamics (multispecies interactions, ...), (C) stock dynamics (recruitment, cannibalism, ...), and (D) harvest patterns. There is an ongoing debate on which of these factors may lead to overcompensatory or compensatory density dependence in populations. In [3] and [9] and references therein, it is argued that climate forcing or changes in productivity regions (factor A) can destabilize populations by introducing overcompensatory dynamics. On the other hand, the view of the authors of [12,14,42] appears to be that populations in communities of ecologically similar species fluctuate in a compensatory way, that is, increases of abundance of one species are accompanied by decreases of other ecologically equivalent species (factors B, C). The purpose of this paper is to study all these effects on the dynamics, and in doing so we shall consider and analyze a discrete density dependent prey-predator model where both species are exposed to harvest. In our analysis we shall use the general Deriso-Schnute relation [32,43] which covers compensatory as well as overcompensatory density dependence in the recruitment.
The plan of the paper is as follows: In Sect. 2 we present the model and analyze equilibria and stability. In Sect. 3 we investigate possible dynamics without harvest, while harvest is included in Sect. 4. Finally, in Sect. 5 we summarize and conclude.

The model, equilibria and stability
Let x t and y t be the sizes of a prey and a predator population at time t respectively, and let h 1 and h 2 (0 ≤ h 1 , h 2 ≤ 1) be the fractions of the populations (harvest) which are removed from the populations at each time. The relation between the two species at two consecutive time steps (years) is assumed to be in the form Moreover, natural restrictions to impose on the functions are Biologically, (2) implies that intraspecific competition leads to a decrease in size of both populations x and y, while interspecific competition (predation) results in a decrease of the survival of the prey x and an increase of the size of the predator population y. By intraspecific competition we mean crowdening effects, possible cannibalistic behaviour or any other self-regulatory mechanism such that ∂ f /∂x ≤ 0, ∂ g /∂ y ≤ 0. In this paper we shall assume that f and g belong to the Deriso-Schnute family [32,43], hence we shall scrutinize the model where −1 ≤ γ < 0. The 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 that α ≥ β, which means that we are considering the part of the parameter space where the predator gains less (or equal) from the prey than the prey is able to "give". Note that when γ = −1,  (3), boundedness is ensured and (depending on parameters), the forward trajectory (x n , y n ) of any initial point (x 0 , y 0 ) where x 0 > 0, y 0 > 0, will remain in the first quadrant x ≥ 0, y ≥ 0 for all times n = 0, 1, 2, . . .. Model (3) possesses three equilibrium points. The trivial one (x, y) , 0), and the nontrivial fixed point where x * satisfies the equation In order to study stability properties we linearize about each of the equilibria, which in all three cases results in an eigenvalue equation in the form As it is well known, an equilibrium (x,ŷ) of a discrete two-dimensional system like (3) is locally asymptotic stable if all eigenvalues of the linearization are located on the inside of the unit circle in the complex plane, and moreover, this is ensured whenever the Jury criteria 1 + a 1 + a 2 > 0, 1 − a 1 + a 2 > 0 and 1 − |a 2 | > 0 hold. If the first of these inequalities fails, as a result of increasing a parameter (F or G in the cases we shall consider), such that 1 + a 1 + a 2 = 0, then λ becomes equal to unity, which in the generic case means that (x,ŷ) will undergo a saddle node bifurcation at instability threshold. If 1 − a 1 + a 2 becomes zero then λ = −1 and a flip (also known as a period doubling) bifurcation occurs. Hence, in the supercritical case, when (x,ŷ) loses its hyperbolicity and fails to be stable, a stable 2-cycle is established. In the subcritical case the 2-cycle is not stable. Finally, when 1 − |a 2 | becomes zero the solution of (6) is on the form λ = e 2πiθ and subsequently (x,ŷ) faces a Neimark-Sacker bifurcation at threshold which in the supercritical case implies that when (x,ŷ) fails to be stable, an invariant attracting curve is established about (x,ŷ).

The dynamics, no harvest
The dynamics generated by model (3) strongly depends on the values of γ (i.e. the degree of compensatory or overcompensatory recruitment). We shall now explore this further and in doing so we will in this section assume no harvest, thus we consider For comparative purposes we start by discussing the case α = β = 1, (the case α > β will be considered later), and finally we also find it convenient to define γ = − 1 /n, n ≥ 1, and use n instead of γ . (n = 1 corresponds to the Beverton-Holt case, n → ∞ to the Ricker case).
Under these restrictions it follows directly that (x, y) = (0, 0) is stable for 0 < F < 1 and any value of G. (x,ỹ) can be expressed as (x,ỹ) = (n( n √ F − 1), 0) where F > 1, and from (7a, 7c) we find that it is stable whenever G < F(F − 1) −1 and F < (n(n − 2) −1 ) n , n > 2. When n ≤ 2 there is no restriction on F other than F > 1. When n = 3 and n = 4 the restrictions are 1 < F < 27, 1 < F < 16 respectively, and n → ∞ implies 1 < F < e 2 . Hence an increase of n (or γ ) leads to a smaller stability region of (x,ỹ). Regarding the nontrivial equilibrium, confer (4) and (5), we may through the definitions and G > F(F −1) −1 is necessary for (15) to be feasible. Stability is ensured whenever (10) holds, i.e. as long as For those values of n (or γ ) where L 1 = x * y * has a solution, (x * , y * ) = ( p, q) satisfies and for those parameter values where x * y * = L 2 , the threshold solution (x * , y * ) = (r , s) may be obtained from Note that when F → (n(n − 2) −1 ) n it follows from (17) and (18) (17) and (19), the point of intersection between the curves L 1 = x * y * and x * y * = L 2 must satisfy from which we conclude that there is no intersection when 1 ≤ n ≤ 2. When n = 1 or 2 the left inequality of (16) is trivially satisfied. Assuming n = 1, the right inequality may be written as which holds as well. The corresponding inequality in the n = 2 case becomes Here we observe that c → 1 (i.e. y * → 0) implies b → F, so (23) degenerates to When n > 2, (γ > − 1 /2), there is a profound change of dynamics. Here, in contrast to the n = 1 and n = 2 cases, parameter values exist where both L 1 = x * y * and x * y * = L 2 . This is displayed in Fig. 1, in case of n = 3, n = 4 and n → ∞. On the lower curve we find the combinations of F and G such that L 1 = x * y * , and the dominant root of the associated eigenvalue equation is λ = −1. The upper curve shows the corresponding case x * y * = L 2 where λ is complex valued, |λ| = 1 and given by (13). The region where (x * , y * ) is locally asymptotic stable is located above the lower curve and below the upper curve, and evidently, the larger the F, the smaller the stable parameter region. Hence, an increase of F acts in a destabilizing fashion. Depending on G, the largest interval where (x * , y * ) is stable is given by 27 < F < 885 when n = 3, and 16 < F < 288.5 when n = 4. When n → ∞ (or In order for (x * , y * ) to be stable we must have L 1 < L 2 which implies x * + y * < 4, and since x * + y * = ln F it follows that F < e 4 . Moreover, (9b) may be expressed as Our next goal is to describe the dynamics in somewhat more detail, and we start by the n = 3 case. When 27 < F < 885, (x * , y * ) is unstable for G values below the lower curve. On the curve, (x * , y * ) undergoes a (supercritical) flip bifurcation and the dynamics in the unstable parameter region is stable period 2 orbits. In Fig. 2 where (n, F, G) = (3, 500, 7) we show an orbit starting at (x 0 , y 0 ) = (8, 2) which converges to a 2-cycle. For the parameter values at hand the basin of attraction for the cycle appears to be R 2 + . We have not detected any stable periodic orbits of period 2 k where k > 1.
Above the curve, in the F interval 27 < F < 294, we find that (x * , y * ) may be both locally asymptotic stable as well as globally stable. Indeed, assume F fixed. For G values below a critical value G C , (x * , y * ) is globally asymptotic stable. However, in a nonlinear system like the one we are scrutinizing here one may not rule out the possibility of multiple attractors. What we find is that when G = G C , the third iterate of (14) undergoes a saddle node bifurcation where both a stable and an unstable 3cycle is established. Consequently, there exists a parameter interval G > G C where there is coexistence between the stable fixed point and a stable 3-cycle, so the ultimate fate of an orbit depends on the initial condition. If we continue to increase G the 3-cycle eventually turns to a chaotic attractor, see Fig. 3a, but still a certain kind of 3periodicity also persists in the chaotic regime. Evidently, there also exists a parameter interval where (x * , y * ) coexists with the chaotic attractor as well. The trapping region for (x * , y * ) is shown in Fig. 3b when (n, F) = (3, 150) and 24 < G < 38. In the remaining part of the figure we find points (x 0 , y 0 ), depending on G, which converge towards the 3-cycle or the chaotic attractor. Through further enlargement of G the chaotic attractor disappears. Our conjecture is that this happens when the attractor and the unstable branches of the 3-cycle collide. If we continue to increase G, (x * , y * ) is the only stable attractor. These findings are summarized in the bifurcation diagram, see Fig. 4. Next, confer Fig. 1a, and assume 294 < F < 885. Then (x * , y * ) is stable between the curves. For a given value of G < G M = 14.077, an increase of F leads to an intersection between L 1 and x * y * , and the dynamics (2-periodical behaviour) in the unstable region below the lower curve has already been accounted for. If G > G M is kept fixed, an increase of F makes x * y * and L 2 intersect, and on the curve x * y * = L 2 , |λ| = 1 and complex valued. Consequently, (x * , y * ) experiences a Neimark-Sacker bifurcation. Hence, beyond stability threshold, the dynamics is quasiperiodic and restricted to an invariant curve, see Fig. 5a.
As we continue to increase F, the invariant curve turns through frequency locking into a stable period 5 orbit, cf. Fig. 5b ,and the corresponding trapping region (shaded region) is displayed in Fig. 5c (in the non-shaded part we find initial points which settle on a 3-cycle. Note that for most of these points, y 0 > x 0 which biologically is not very relevant). At even higher F values we find chaotic dynamics, cf. Fig. 5c, as well as stable periodic orbits of long period. The dynamics reported above may also be described by use of the Lyapunov exponent L, and as it is well known, cf [44],  L < 0 corresponds to a stable equilibrium or a stable periodic orbit. L = 0 means that we have a quasiperiodic orbit restricted to an invariant curve while L > 0 indicates a chaotic orbit. In Fig. 5d we show the values of L in the case of 350 < F < 550 and G = 22.5. When F < 390, L < 0 and (x * , y * ) is stable. Whenever 390 < F < 440, L = 0 and the dynamics is restricted to an invariant curve. The stable 5-cycle is found when 440 < F < 490, (L < 0). Chaotic dynamics occurs in tiny parameter intervals about F = 500 where L > 0. Finally, let us scrutinize γ → 0 (n → ∞). In this case model (14) may be expressed as which we have partly studied in a previous paper, see [21]. Therefore, for the sake of completeness and comparative purposes, we shall present some of our findings from [21], but several computational details will be omitted. The fixed point (x,ỹ) may be written as (x,ỹ) = (ln F, 0), and it is stable whenever 1 < F < e 2 and G < F(F − 1) −1 , which is the smallest possible stable region as long as −1 < γ < 0. The nontrivial fixed point becomes (x * , y * ) = (ln b, ln c) where G > F(F −1) −1 and x * + y * = ln F. Regarding stability properties we refer to Fig. 1c which should be compared to Fig. 1a, b. Depending on G, the largest F interval where (x * , y * ) may be stable is e 2 < F < e 4 . The stable region is located between the curves, and just as in Fig. 1a and 1b, the lower curve consists of those parameter combinations such that L 1 = x * y * , while the upper curve consists of parameter values such that x * y * = L 2 . We also find that an increase of F, G fixed, results in a smaller parameter region where (x * , y * ) is stable.
At threshold L 1 = x * y * the solution of the eigenvalue equation is λ 1 = −1, λ 2 = 3 − ln(bc) and |λ 2 | < 1. The corresponding threshold values become Note that when F → e 2 then ( p, q) → where T is a 2 × 2 matrix whose columns are the eigenvectors corresponding to λ 1 = −1 and λ 2 = 3 − ln F of the Jacobian evaluated at threshold L 1 = x * y * , it is possible to rewrite (24) in "standard form" as H and Q are polynomials which contain second and third order terms. Moreover, the restriction of (26) to the center manifold gives and since 1 2 cf [21], we conclude that the bifurcation is supercritical. A similar proof in case of γ arbitrary, γ ∈ (− 1 /2, 0), seems out of reach.
Regarding the dynamics in the unstable parameter region below the curve where L 1 = x * y * , it is much richer here than in the n = 3 case. Indeed, when n = 3 we find period 2 orbits only. Here, there are additional periodic orbits of period 2 k , k = 2, 3, . . ., as well as chaotic dynamics. An example of a chaotic attractor is shown in Fig. 6a. In the case of F sufficiently large there is also a possibility that the predator population may be driven to extinction, and the larger the F the larger the possibility. The mechanism is visualized in Fig. 6b where (F, G) = (66, 4). Actually, this may happen also for other values of γ given that they are close to zero, as shown in Fig. 6c. For several iterations both populations perform chaotic oscillations similar to what is shown in Fig. 6a, but once x falls below a critical value x C we observe a profound change of dynamics. When x < x C the predator population becomes very small, actually so small that it does not manage to recover, and subsequently, in case of x small, (24) degenerates to x t+1 = F x t , from which we conclude that the prey indeed manages to recover and may in fact be large before it is damped again by the factor e −x . Hence, in case of F larger than a critical value F C , the only attractor is the one where the prey shows highly oscillatory behaviour and y = 0. For all practical purposes we may consider this attractor as being generated by x → Fe −x x.
At threshold x * y * = L 2 (x * , y * ) = (r , s) = 1 2 ln(bc) + √ d, In [21] it was proved that (x * , y * ) undergoes a supercritical Neimark-Sacker bifurcation when x * y * = L 2 . Regarding the nonstationary dynamics we find just beyond instability threshold quasiperiodic orbits restricted to invariant curves, just as in the n = 3 (γ = − 1 /3) case. Through further enlargement of F, keeping G fixed, the dynamics is captured by use of Lyapunov exponent calculations and bifurcation diagrams, see Fig. 7. As both figures demonstrate, periodical dynamics is not as pronounced here compared to the corresponding n = 3 case. However, a stable period 2 orbit, created through frequency locking (73 < F < 87), is clearly visible. There are two qualitative different possibilities of chaotic dynamics just as we found in the parameter region below the curve L 1 = x * y * . We find orbits where both populations perform chaotic oscillations, but also large amplitude chaotic prey oscillations generated by x → Fe −x x where y = 0. These two scenaria are displayed in Fig. 8a, b respectively. In Fig. 8a where (F, G) = (72, 5) both species survive. In Fig. 8b where (F, G) = (95, 5) the predator has died and the prey exhibits chaotic oscillations. Finally, let us close this section by briefly considering the impact of different interaction parameters α and β. Assuming γ = −1 (or n = 1) and (α, β) = (1, 1), the nontrivial equilibrium (x * , y * ) (1,1) is given by (15). If (α, β) and straightforward calculations show that In a similar way, if γ → 0 we find, by use of the abbreviation u = (4F + G)G −1 , that From a biological point of view the results above make sense. Indeed, α > β means that the predator will benefit less by eating compared to the symmetric case α = β, which in turn will lead to a decrease of the size of the predator population. Consequently, the predation pressure on the prey becomes smaller, allowing the prey population to increase. Moreover, as already proved, cf. Eqs. (12,13), the spectral radii are −1 and e iθ at instability thresholds L 1 = x * y * and x * y * = L 2 respectively. Hence, regarding the nonstationary behaviour, there are no reasons to expect severe qualitative dynamical changes between the cases α = β and α > β. However, there are some changes. In Fig. 9a, where n = 3 (γ = − 1 /3), we show the graphs of the curves L 1 = x * y * and x * y * = L 2 when (α, β) = (1, 1), (1, 3 /4) and (1, 1 /2), and in Fig. 9b the same curves are displayed as γ → 0. Considering x * y * = L 2 in the n = 3 (γ = − 1 /3) case, there is a clear tendency that the size of the stable parameter region shrinks as β becomes smaller. This is hardly not the case when γ → 0. Thus, one may argue that a reduction of β acts destabilizing for small values of n (n ≥ 3). On the other hand, the curves L 1 = x * y * are quite similar in all cases, which clearly suggests that the values of F and G at lower instability threshold are relatively unaffected by the values of β.

The impact of harvest
As documented in the previous section, there is a profound variation of dynamical behaviour when γ varies in the interval [−1, 0). In this section our intention is to study the impact of harvest and we start by scrutinizing the case γ → 0, i.e. the model Applying the notationsF = F (1 − h 1 , the nontrivial equilibrium of (31) may be written as where x * obeys the equation Moreover, in order to relate results here to results obtained in the last section, we find it natural to assume α = β = 1. Then, by use ofb = (F +Ĝ)Ĝ −1 ,ĉ =FĜ(F +Ĝ) −1 , (x * , y * ) = (lnb, lnĉ) (33) and stability of (33) is guaranteed whenever where We find it natural to consider these three cases: Assuming (A), then First, suppose F and G values such that x * y * = L 2 in case of no harvest, which implies (ln b)(ln c) = cG −1 ln(bc). Then, as harvest is introduced, the term ln(1 − h) ln b in x * y * becomes much more negative (recall that ln b ≥ 2, cf. the text after Eq. (25)) than the term cG −1 ln(1 − h) in L 2 . Thus, an increase of h acts in a stabilizing fashion. Indeed, if (F, G) = (43, 5), the dynamics generated by (31) in case of no harvest is a quasistationary orbit restricted to an invariant curve, while the dynamics is nothing but a stable fixed point when h = 0.1. On the other hand, if F and G are chosen in such a way that L 1 = x * y * in absence of harvest, we find that an increase of h acts destabilizing. This is due to the fact that the term 2cG −1 in L 1 is less than 2, while ln b > 2 in x * y * . Consequently, as h is increased, instability threshold is shifted towards a smaller parameter region. Example: If (F, G, h) = (40, 4, 0) the dynamical outcome of (31) is a stable fixed point, while (F, G, h) = (40, 4, 0.1) leads to a stable 2-cycle. Finally, recall from Sect. 3 (no harvest) that whenever G < G C the transfer from stability to instability goes through a supercritical flip bifurcation, while an increase of F leads to a supercritical Hopf bifurcation when G > G C . Moreover, as shown, in both cases (G < G C , G > G C ) large F values result in chaotic oscillations, and if F is sufficiently large the predator population may actually face extinction. In this part of parameter space an increase of h acts in the same destabilizing fashion, as exemplified in Fig. 10. In Fig. 10a, where h = 0, both populations x and y exhibit chaotic oscillations. In Fig. 10b, where h = 0.2, the predator dies and subsequently the prey starts to perform large amplitude chaotic oscillations. Thus, under the assumption that the populations are located in the chaotic regime, an enlargement of F or introduction to harvest may both drive the predator to extinction.
Regarding case (B), The impact of harvesting only the prey population is the following: (i) Assume that F is large, then h 1 > 0 will reduce the size of x, which in turn will prevent chaotic prey oscillations and consequently reduce the possibility that the predator population will fall below a critical value from where it does not manage to recover. (ii) Suppose parameter values F and G such that x * y * = L 2 in case of no harvest. Then, by comparing terms as accounted for in (A), we find that harvest acts stabilizing here too.
(iii) Finally, assume F and G such that L 1 ≈ x * y * . SinceF(F +G) −1 < F(F +G) −1 both terms in L 1 contribute to a decrease of L 1 when h 1 > 0. On the other hand, the product of the two first logarithmic terms in x * y * increases while the product of the last becomes small as h 1 is increased. Thus, harvesting the prey population only acts stabilizing. Considering (C), Comparing the limits L 1H , L 2H above with the corresponding limits in case of no harvest, it is immediately clear that L 1H > L 1 and L 2H > L 2 when h 2 is increasing, and that x * y * is decreasing. In consequence, if a population exhibits nonstationary oscillations in case of no harvest, the dynamics will shift to a stable fixed point when h 2 is increased. If we continue to increase h 2 , x * y * becomes smaller than L 1H , hence a stable 2-cycle is soon the dynamical outcome. Thus, depending on the dynamics without harvest, the introduction to harvest may act in a stabilizing as well as a destabilizing fashion. There is a rapid change of dynamics when h 2 is increased. Indeed, if (F, G, h 2 ) = (43, 5, 0) → (43, 5, 0.1) → (43, 5, 0.2) → (43, 5, 0.4) the dynamics changes from an invariant curve to a stable fixed point, a 2-cycle and eventually to chaotic dynamics where only the prey survives, as visualized in the bifurcation diagram, Fig. 11. Fig. 12 shows attractions in the chaotic regime. In Fig. 12a where h 2 = 0.35 both species survive, while Fig. 12b where h 2 = 0.4 displays the "ultimate" case where the predator has died. When γ = − 1 /n the expressions corresponding to (33) and (34) may be cast in the form x * y * = n 2 1 + n F − n b + n √ĉ (37) Next, let us focus on the n = 3 (γ = − 1 /3) case. Assuming equal harvest rates h, we find that an increase of h acts in a stabilizing way, provided x * y * > L 2 in case of no harvest. Indeed, assuming (F, G) = (500, 22.5) cf. Fig. 5d, the populations exhibit chaotic oscillations when h = 0. If h = 0.1 the chaotic behaviour disappears and the dynamics is restricted to an invariant curve (quasiperiodic orbits), similar to what is displayed in Fig. 5a. A stable equilibrium (x * , y * ) is established when h = 0.2. Thus, in this part of parameter space, harvest acts as a strong stabilizing effect. On the other hand, if F and G have values such that there is 2-periodic behaviour in case of no harvest, we actually find that the same qualitative picture persists when h is increased. Hence, harvest is neither stabilizing nor destabilizing. If the populations perform chaotic oscillations as shown in Fig. 3 (h = 0), we observe that an enlargement of h is stabilizing as well. h = 0.1 results in a stable 3-cycle and when h = 0.2 the 3-cyclic attractor disappears and the only attractor is the equilibrium (x * , y * ).
Under the assumption x * y * > L 2 when there is no harvest, we find that the introduction of harvesting only the prey (h 1 > 0, h 2 = 0) leads to the same qualitative dynamical picture as increasing h = h 1 = h 2 . If map (3) generates 2-periodic dynamics in case of no harvest, an enlargement of h 1 will reduce the amplitude of the oscillations and eventually (x * , y * ) will become stable, given that h 1 is sufficiently large. Regarding the chaotic attractor shown in Fig. 3, it is necessary to use h 1 = 0.3 in order to establish a 3-cycle and h 1 = 0.4 in order to arrive at a stable equilibrium. Hence, in this part of parameter space, an increase of h 1 acts in a rather weak stabilizing manner.
Finally, consider h 1 = 0 and h 2 > 0. Again, assuming parameter values F and G as in Fig. 5a, an increase of h 2 acts stabilizing whenever h 2 is small, but after a stable equilibrium (x * , y * ) is established, further increase of h 2 results in both 2-an 4-periodic dynamics. During this process the size of of x increases while y decreases, and the predator faces extinction. Turning to the chaotic attractor displayed in Fig. 3, it disappears immediately when h 1 = 0 and h 2 > 0, the only dynamics is a stable equilibrium. The rationale behind this is that the third iterate of (3) will not undergo a saddle node bifurcation when h 1 = 0 and h 2 > 0.

Discussion
In this paper we have analyzed a discrete prey-predator model where both compensatory and overcompensatory density dependence as well as harvest are included. The "degree" of compensatory or overcompensatory density dependence is measured by the parameter γ (or n) and −1 ≤ γ < 0 (1 ≤ n < ∞). Moreover, the parameter space of the full model is huge, therefore results are mostly presented by use of selected parameter values in combination with numerical experiments.
Considering F and G values such that x * y * < L 1 (i.e. in the region below the lower curves in Fig. 1a, c), there is 2-periodic dynamics in case of γ (or n) small, but through an enlargement of γ (or n), it is possible to find orbits of period 2 k as well as chaotic oscillations with large prey amplitudes and small predator amplitudes respectively. For parameter values F and G such that x * y * > L 2 (i.e. in the region above the upper curves in Fig. 1a, c), the dynamics is quasiperiodic and restricted to invariant curves when x * y * > L 2 and |x * y * − L 2 | is small. When F is increased such that the difference |x * y * − L 2 | becomes larger, one may by use of Lyapunov exponent calculations detect chaotic dynamics, in tiny intervals only when γ (or n) is small, (n = 3), in larger intervals when γ → 0. In [45], where a series of species ranging from trees, plants, moths, fish, small mammals, rodents and butterflies have been under consideration, the authors concluded that compensatory dynamics are rare in natural ecological communities. The results from our study support their conclusion. The highly oscillatory behaviour detected in several fish species [46], small rodent population [47], and other species do not seem to be a result of compensatory density dependence. We should also stress our finding that sufficiently large values of F in the overcompensatory case may actually lead to extinction of the predator population. This will not occur assuming compensatory density dependence.
Finally, let us turn to the impact of harvest. Assuming no harvest of the predator population, an increase of h 1 acts in a stabilizing fashion exclusively whenever γ ∈ (− 1 /2, 0) (or n > 2). In particular, if F is large, small values of h 1 will effectively prevent chaotic prey oscillations in the γ → 0 case, which in turn will reduce the risk that the predator population will fall below a critical limit from where it does not manage to recover. On the other hand, harvesting only the predator population may act both stabilizing and destabilizing. For any γ ∈ (− 1 /2, 0) (or n > 2), we find that as long as x * y * > L 2 and |x * y * − L 2 | is small in case of no harvest, an increase of h 2 leads to a stable equilibrium. It acts destabilizing in all other cases, and it is indeed possible to drive the predator population to extinction as accounted for in Sect. 4. Considering the symmetric case h = h 1 = h 2 , the dynamical consequences of increasing h on the whole is quite similar to what we found by increasing h 2 only. In parameter regions where an enlargement of h 1 tends to stabilize the dynamics while an increase of h 2 acts in a destabilizing way, we find that an increase of h qualitatively matches an increase of h 2 . The main difference really is that the changes from stability to nonstationary and chaotic dynamics are somewhat slower. A similar phenomenon has also been observed in density dependent one-population models with age-structure, cf. [34]. Indeed, considering two age-classes, applying harvest rates h 1 , h 2 at the age-classes respectively, there are parameter regions where an increase of h 1 is stabilizing, h 2 destabilizing, while an increase of h = h 1 = h 2 turns out to have a weak destabilizing effect.
Summary: The main findings of our analysis is as follows. The parameter regions where equilibria are stable shrink as density dependence shifts from being compensatory to overcompensatory. Moreover, an increase of the fecundity of the prey will in general lead to nonstationary dynamics. Depending on the degree of overcompensation, such an increase may actually introduce chaotic dynamics where oscillations become so severe that the predator population faces extinction. Turning to harvest, if only the prey is the target, we show, independent of compensatory/overcompensatory density dependence, that increased harvest acts in a stabilizing fashion. On the other hand, if we harvest from the predator population only, this may in certain parameter regions act as a strong destabilizing effect. In the interplay between these two harvest strategies the latter seems to dominate.

Compliance with ethical standards
Conflicts of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.