Cooperative hunting in a discrete predator–prey system II

ABSTRACT We investigate a discrete-time predator–prey system with cooperative hunting in the predators proposed by Chow et al. by determining local stability of the interior steady states analytically in certain parameter regimes. The system can have either zero, one or two interior steady states. We provide criteria for the stability of interior steady states when the system has either one or two interior steady states. Numerical examples are presented to confirm our analytical findings. It is concluded that cooperative hunting of the predators can promote predator persistence but may also drive the predator to a sudden extinction.


Introduction
Cooperation is frequently observed and widespread among individuals of social animals in many biological systems. For example, carnivores such as wolves, wild dogs and lions often work together to capture and kill their preys [8]. Other organisms such as spiders, birds and ants also seek and attack prey collaboratively to increase their hunting success [9].
Earlier research incorporating cooperative hunting includes Berec [3] who uses ordinary differential equations to model predator-prey interactions with a Holling type II functional response. Due to this functional response, Berec studies the effects of cooperative hunting relative to population oscillations. Cosner et al. [5] on the other hand propose models of partial differential equations to explore the effects of predator aggregation when predators encounter a cluster of prey. Recently, Alves and Hilker [2] use models of ordinary differential equations of predator-prey interactions with cooperative hunting in predators to investigate impacts of cooperative hunting. In the absence of predators, the prey population is governed by the logistic equation. They conclude that cooperative hunting can improve persistence of the predator but may also promote a sudden collapse of the predator. Further, their research suggests that cooperative hunting is a mechanism for inducing Allee effects in predators.
Since the pioneer work of May [7], mathematical models of difference equations have played important roles in the understanding of population interactions. There are many populations in nature with non-overlapping generations and continuous-time models are not adequate to describe such populations. In addition, data collected in ecological studies are usually in discrete formats. Consequently, discrete-time models are more appropriate to study such population interactions. Motivated by these, Chow et al. [4] propose and investigate a discrete-time predator-prey model with cooperative hunting among predators to study the effects of cooperation upon predator-prey interactions. Their model derivation is built on the well-known Nicholson-Bailey system with density-dependent prey growth rate and it is proven that both populations coexist indefinitely if the maximal reproductive number of the predator is larger than one. The interaction may support two coexisting steady states if predator's maximal reproductive number is less than one but with intense cooperative hunting among predators.
In this study, we investigate local stability of the coexisting steady states when predators cooperative intensively. In particular, we prove that one of the coexisting steady states is always a saddle point, while the other steady state may be a repeller or can change its stability from asymptotically stable to a repeller. Numerical investigations will be performed to verify our analytical findings and to provide further explorations.
In the following section, we briefly review the discrete-time model and its prior results. Section 3 presents our main results by determining local stability of the interior steady states. Numerical simulations are given in Section 4 and the final section provides a brief summary and discussion.

Review of notations
In this section, we briefly review the model and its prior results obtained in [4]. In particular, global dynamics of the model are stated when the degree of cooperation is small and the number of interior steady states are determined when predators engage in cooperative hunting intensively.
We first go over the model derivation. A more detailed description is given in [4]. Let x(n) and y(n) denote, respectively, the prey and predator populations at generation n = 0, 1, 2, . . .. In the absence of cooperative hunting and by applying a similar argument as in the derivation of Nicholson-Bailey model [1], the number of encounters between prey and predators in generation n is assumed to follow the law of mass action, ax(n)y(n), where the constant a > 0 denotes searching efficiency of the predators. It is also assumed that the number of encounters is distributed randomly and follows a Poisson distribution with probability p(m) = e −μ μ m /m!, where m = 0, 1, 2, . . . is the number of encounters and μ is the average of encounters per prey per generation. Therefore, μ = ax(n)y(n)/x(n) = ay(n) and thus 1 − p(0) = 1 − e −ay(n) is the probability of an individual prey being preyed upon in generation n.
With cooperative hunting, the number of encounters between prey and predators at time n becomes ax(n)y(n)(1 + αy(n)), where α ≥ 0 denotes degree of cooperative hunting. There is no cooperation among predators if α = 0 and cooperation is stronger if α is larger. It follows that the probability of an individual prey escaped from being preyed upon at time n is e −ay(n)(1+αy(n)) . The probability is smaller due to cooperation among predators. Further, we assume that the prey in the absence of predators is modelled by the Beverton-Holt equation. The Beverton-Holt model is often considered as the discrete analogue of the continuous-time logistic equation which is used in the study of cooperative hunting by Alves and Hilker [2].
Putting all of these together, the interaction between the two populations proposed in [4] is described by the following system with nonnegative initial conditions, where λĝ(x) = λ/1 + kx, λ, k > 0, is the prey's per capita growth rate. The parameter β > 0 is the predator conversion for each prey consumed.
To analyse the model, we nondimensionalize system (1) by letting Ignoring the tildes, (1) is converted into the following system with only three parameters In the following, we briefly summarize prior results of system (3). In particular, we have shown in [4, Proposition 2.1] that the extinction steady state E 0 = (0, 0) is globally asymptotically stable if λ < 1 and it is globally attracting if λ = 1. Throughout this paper, we assume λ > 1. Then (3) has an additional boundary steady state E 1 = (x, 0), wherex = λ − 1 is the prey's carrying capacity. If 2α ≤ 1, Theorem 4.2 of [4] proves that the prey existence steady state E 1 = (x, 0) is globally asymptotically stable if βx < 1, which extends a previous result in [6] for α = 0. Notice that βx can be interpreted as the predator's maximal reproductive number as it is the predator's reproductive number when the prey is stabilized at its carrying capacitȳ x. If βx > 1, then (1) has a unique interior steady state E * = (x * , y * ) and there exists a unique β d > 0 such that E * is asymptotically stable if β < β d and a repeller if β > β d . The interior steady state undergoes a Neimark-Sacker bifurcation at β = β d . See Theorems 4.2 and 4.3 in [4]. Moreover, it is shown in [4, Theorem 2.2] that the system is uniformly persistent if βx > 1. Consequently, dynamics of the interaction are similar to the model of no cooperation if 2α ≤ 1. That is, cooperative hunting has no effects on the population dynamics if its magnitude is small.
When the degree of cooperation is in the middle range, 1 < 2α < (3λ − 1)/(λ − 1), it is shown in [4] that all of the above results for 2α ≤ 1 remain valid with one exception, namely the global asymptotic stability of E 1 = (x, 0) for βx < 1 is stated as a conjecture since we can only verify the statement in a narrow parameter region with βx √ 2αe (1−2α)/4α < 1. If the degree of cooperative hunting is large, 2α > (3λ − 1)/(λ − 1), the number of interior steady states can be either zero, one or two [4], Theorem 3.2]. In particular, the system may have two interior steady states when βx < 1 for which the predator would otherwise go extinct in the absence of cooperative hunting. Therefore, cooperation can mediate coexistence in the interaction. In this investigation, we will study local stability of the interior steady states when 2α > (3λ − 1)/(λ − 1), which is not explored in Chow et al. [4].
The nontrivial x-isocline is given by with f (y) = −λae −♦ (1 + 2αy) < 0. Define by solving f (y) = 0, i.e. e ♦ = λ. Then To be biologically feasible, we consider only y ∈ (0, y c ) for the existence of interior steady states since the prey population in the steady state would be negative by (7) if y > y c . For the convenience of the reader, the existence and number of the interior steady state derived in [4,Theorem 3.2] is restated as follows.
A key ingredient in the proof of Theorem 2.1 is that on the interval [0, y c ), forms a family of nonintersecting curves and for fixed 0 ≤ y < y c , h β (y) ↓ strictly to 0 when β ↑ ∞. It follows from (11) and (12) that y * 1 ↓ strictly and both y * 2 , y * ↑ strictly as β increases.
We use Figure 1(a) to illustrate Theorem 2.1 with λ = 10 and α = 30. Then 2α − (3λ − 1)/(λ − 1) ≈ 56.8 > 0. The two isoclines are tangent to each other at β = β * ≈ 0.0507, and there is a unique positive intersection when β = 1/x ≈ 0.111. It is clear that the system has two interior steady states if β is in (0.0507, 0.111) while there is a unique interior steady state if β > 0.111 and there is no interior steady state if β < β * = 0.0507. Once existence and the number of interior steady states are determined, we proceed to review notations relative to the study of local stability of the steady states. With fixed λ > 1, Theorem 2.1 means that any point (x, y) on x = f (y) with 0 < y < y c will be an interior steady state for a certain β. Using x = f (y) we can rewrite det(J) and tr(J) of the Jacobian matrix as functions of y, Since tr(J) > 0 on [0, y c ) and by the Jury condition, we need only to study det(J) and in order to determine local stability of (x, y) on x = f (y). By a straightforward calculation, we verified the following in [4].
As a result of Lemma 2.2, det(J(y)) < 1 for y ∈ (0, y d ) and det(J(y)) > 1 for y ∈ (y d , y c ). For the other inequality in the Jury condition, it is shown in [4] that V(y) satisfies the following.

Main results
In the previous section, we have briefly introduced four critical points y * , y e , y d and y t in (0, y c ). In this section, we will investigate their order relations from which local stability of an interior steady state (x, y) can be determined directly. For instance, Jury condition and Lemmas 2.2 and 2.3 imply (x, y) is asymptotically stable if y ∈ (y t , y d ), which makes sense only if y t < y d .
We will first prove in Section 3.1 that y * = y t . The comparison between y t and y d is given in Section 3.2 and the order of y d and y e is studied in Section 3.3. Local stability of the interior steady states is presented in Section 3.4. These proofs depend heavily on Lemmas 2.2 and 2.3. However, the technique explored in this investigation is very different from that employed in Chow et al. [4].

Comparison between y d and y *
In this subsection, we will compare y d with y * . By Proposition 3.1 and (11), Using Lemmas 2.2 and 2.3, we can then determine immediately the local stability of any interior steady states. To this effort, we take the advantage of the new variable s defined in Section 3.1. Define It is easy to check that on [1, ∞), In particular, f d (s) is strictly increasing and concave on [1, ∞). Moreover, By Lemma 2.2 and (21), For λ > 1, a straightforward calculation shows that 1 < λ * =: as the third root (2λ 2 + 3λ + λ √ 4λ 2 − 3)/2(λ + 1) > λ. We now prove that the order between y d and y * is equivalent to the order between f d (λ * ) and t(λ * ).
holds for λ large enough. The claim is proven.

Comparison between y d and y e
In contrast to (25), if y * < y d , then either y * = y t < y d < y e or y * = y t < y e ≤ y d .
So we need to compare y d with y e before using Lemmas 2.2 and 2.3 to determine local stability of any interior steady states. We shall mimic the method used in the proof of Theorem 3.2. For s ≥ 1, definẽ Hence,g (s) = 2 / s + (λ + 1)/(λ − 1)s 2 − 2λ/(λ − 1)s 3 and Moreover,g Note that by (17), We claim that Recall that E * = (x e , y e ) is the unique interior steady state when β = 1/x = 1/(λ − 1). That means y e is the unique positive solution to h(y) = f (y) when β = 1/(λ − 1).
Using (18) and (16), we get Note thatλ * depends only on λ. Therefore similar to (34), we have Using (38) and with the replacement of (24) by (39), we can repeat the same arguments as in Theorem 3.2 to obtain the following result. Computer simulation indicates that both y d > y e and y d < y e can happen depending on parameter regimes, which will be shown theoretically as follows.

Local stability of the interior steady states
In this subsection, we provide local stability of the interior steady states. Recall from Theorem 3.2 that there are two interior steady states E * 1 = (x * 1 , y * 1 ) and E * 2 = (x * 2 , y * 2 ) where y * 1 < y * 2 when β ∈ (β * , 1 /x ) and there is only one steady state E * = (x * , y * ) when β ≥ 1/x. In particular, (11) holds. We shall separate our discussion into y d ≤ y * and y d > y * .
We conclude from Theorem 3.6 that the interior steady state E * 1 that is close to the stable boundary steady state is always a saddle point when the model has two interior steady states.
It remains to consider y d > y * . That is f d (λ * ) > t(λ * ) by Theorem 3.2. We first use Theorem 3.4 to find out the order of y d and y e . Then Lemmas 2.2 and 2.3 can be applied to determine the local stability of E * i and E * .
is asymptotically stable or a repeller depending on y * 2 < y d or y * 2 > y d , while E * is a repeller, , then E * 2 is asymptotically stable and E * = (x * , y * ) is asymptotically stable or a repeller depending on whether y * < y d or y d < y * < y c .
Proof: By Theorem 3.4, we have 0 < y * < y d < y e < y c , 0 < y * < y d = y e < y c and 0 < y * < y e < y d < y c for cases (a), (b) and (c), respectively. Lemmas 2.2 and 2.3 can be applied as in Theorem 3.6. The detail is omitted We conclude from Theorem 3.7 that the unique interior steady state E * is either a repeller or asymptotically stable when β > β c , while E * 1 is always a saddle point and E * 2 can change its stability when β ∈ (β * , β c ).

Numerical simulations
In this section, we use numerical examples to investigate model (3) and to confirm the established analytical findings. Our parameter values adopted here are for demonstration only and are not from any field data.
Recall that y d denotes the positive root of det(J) = 1 while y t is the positive solution of det(J) + 1 − tr(J) = 0. These two roots are clearly smaller than y c . We illustrate numerically that either y d < y t or y d > y t can occur. However, y d > y t is more likely unless λ > 1 is small as shown in Corollary 3.3.
Let λ = 1.5 and α = 12 where λ is small. Thenx = 0.5, β c := 1/x = 2, and 2α − (3λ − 1)/(λ − 1) = 17 > 0. Moreover, y c ≈ 0.1468, y d ≈ 0.04163 < y t = y * ≈ 0.04230 and y e ≈ 0.082675. When β = β c + 0.1, there exists a unique interior steady state E * ≈ (0.25379, 0.08748) which is a repeller with real eigenvalues larger than one. Notice y d < y * and this confirms the result given in Theorem 3.6(b). Figure 1(c) plots E * and a solution of the system with initial condition (0.25479, 0.08848) chose to E * after the first 8000 iterations are removed and the next 5000 iterations are plotted. Observe that β > β c implies βx > 1 and the boundary steady state E 1 = (x, 0) is a saddle point with its stable manifold lying on the positive x-axis. There is no stable steady state in the system and the asymptotic dynamics provided in Figure 1(c) indicates, however, that both populations coexist. This coexistent result has been established theoretically in [4] since βx > 1.
We next decrease β to β c − 0.1. Since β c − 0.1 ≈ 1.9 > β * ≈ 1.680903, there are two interior steady states E * 1 ≈ (0.48932, 0.006622) and E * 2 ≈ (0.294864, 0.076617), where E * 1 is a saddle point and E * 2 is a repeller with real eigenvalues greater than one. The stability results of E * 1 and E * 2 are consistent with Theorem 3.6(a). We use the same initial condition as that in Figure 1(c) and the solution converges to the boundary steady state E 1 = (0.5, 0) as shown in Figure 1(d) where the first 8000 iterations are removed and the next 5000 iterations and E * 2 are plotted. Since β < β c , we have βx < 1 and the boundary state E 1 is asymptotically stable. The only stable steady state in the system is E 1 and it is suspected that solutions except those of E * 1 and E * 2 converge to E 1 . To study the case where y * < y d numerically, we let λ = 55 and α = 10/3. Then 2α − (3λ − 1)/(λ − 1) = 3.6296 > 0, β c = 0.01852 andx = 54. Moreover, y c ≈ 0.95667, y e ≈ 0.205285, y d ≈ 0.312369 > y t = y * ≈ 0.098958, and β * ≈ 0.0169961. With β = β c + 0.003, system (3) has a unique interior steady state E * = (x * , y * ) ≈ (31.14712, 0.278489) which is asymptotically stable with real eigenvalues less than one and E 1 = (54, 0) is a saddle point with stable manifold lying on the positive x-axis. Although not presented, solutions with initial conditions close to E * do converge to E * . Notice that stability of E * confirms Theorem 3.7(c) as y * < y d .
We next let β = β c + 0.0051. Then there is a unique interior steady state E * = (x * , y * ) ≈ (28.05156, 0.3125779), a repeller, which confirms Theorem 3.7(c) since y * > y d > y e . Using an initial condition close to E * , eliminating the first 3000 iterations and plotting the next 2000 iterations, the ω-limit set of the solution is an invariant closed curve as shown in Figure 2(a) where E * is also plotted. We anticipate that a Neimark-Sacker Figure 2. Simulations for the case where y * < y d are presented. In (a) there is a unique interior steady state which is a repeller and the system has an invariant closed curve. In (b)-(d), the system has two attractors and 60,000 randomly generated initial conditions are plotted. Initial conditions convergent to the boundary steady state are denoted by red while initial conditions convergent to the stable interior steady state are marked by blue colour.
bifurcation occurs when the stable interior steady state becomes unstable. Numerical simulations given for Figure 3(a),(b) to be presented below reconfirm the bifurcation.
Finally, we present bifurcation diagrams using β as the varying parameter. Here λ = 55 and α = 10/3, that is, parameter values of those in Figure 2(a ,b) are used. The first 1000 iterations are discarded and the next 1000 iterations are plotted in Figure 3(a ,b) so that only the attractors are shown. Notice that β lies in [0.01, 0.05] containing both β * ≈ 0.0169961 and β c ≈ 0.01852. It is clear that when β is in (0.02, 0.03) the system has a unique interior steady state and a Neimark-Sacker bifurcation occurs when the steady state becomes unstable. Since Figure 3(a),(b) provide only attractors, to illustrate the saddle-node bifurcation at β = β * , we plot steady states when β is in (0.0169962, 0.018525) in Figure 3(c,d). The solid and dashed curves denote stable and unstable steady states, respectively. It is clear in these plots that the system has two interior steady states when β ∈ (β * , β c ) and there is a unique interior steady state when β > β c and β ≈ β c .

Summary and discussion
This investigation is motivated by the recent work of Alves and Hilker [2] who use continuous-time models of predator-prey interactions with cooperative hunting in predators to study the impacts of cooperation upon population interactions. There are many populations with non-overlapping generations and discrete-time models are more appropriate to describe such populations. Further, ecological data collected are usually in discrete formats. Consequently, Chow et al. [4] derive a discrete-time predator-prey system with hunting cooperation within the predators. The discrete model is based on the classical Nichelson-Bailey system with density dependent host growth rate and cooperative hunting of the predator is modelled via the attack rate of the predator. An earlier research [6] proves that predators go extinct when there is no cooperative hunting and βx < 1, where βx may be interpreted as the predator's maximal reproductive number since it is the reproductive number of predators when the prey population is stabilized at its carrying capacitȳ x = λ − 1. The study of [4] further concludes that the predators also go extinct if the maximal reproductive number is less than one and the degree of cooperative hunting is not large, i.e. 2α ≤ (3λ − 1)/(λ − 1).
Under the condition that the maximal reproductive number of predators is less than one, the study carried out in [4] shows that the discrete model (3) can have either zero, one or two interior steady states when the degree of cooperation is large, i.e. 2α > (3λ − 1)/(λ − 1). The present investigation continues the study of the model by establishing local stability of the interior steady states analytically and by exploring asymptotic dynamics of the system when the parameters are within this regime.
We find out with intense cooperative hunting, Theorems 3.6 and 3.7 of the present study suggest that strong cooperative hunting in predators can promote predator persistence when its maximal reproductive number is less than one. Indeed, in this circumstance the interaction has two coexisting steady states if the predator's conversion β exceeds a critical value β * , i.e. β > β * . If predator's conversion is less than β * , then cooperative hunting cannot prevent predator extinction. We prove analytically that the interior steady state that is close to the predator extinction steady state is always a saddle point while the other one can change stability from asymptotically stable to a repeller. Our numerical simulations indicate that the change of this stability is through a supercritical Neimark-Sacker bifurcation. As a consequence, both populations can coexist indefinitely if initial populations are in the basin of stable coexisting steady state or of an attracting invariant closed curve. It follows that cooperative hunting can mediate predator persistence.
To compare our results with those of Alves and Hilker [2], first noticing both studies conclude that cooperative hunting can promote predator persistence. While the conclusion in [2] is drawn from numerical investigations, our result is based on mathematical analysis. In particular, the degree of cooperative hunting in this study is quantified explicitly via α > (3λ − 1)/2(λ − 1). The inequality provides a relation between the degree α of predator cooperation and the growth rate λ of prey population. In addition, we prove that one of the coexisting steady states is always a saddle point, which was not shown analytically in [2]. However, since solutions of continuous-time models are continuous and unique for initial value problems, the stable manifold of the saddle interior steady state in [2] is a continuous decreasing curve separating the positive coordinate plane into two positively invariant regions with the region that lies below the stable manifold resulting in predator extinction. This stable manifold is termed as the predator's Allee threshold by Alves and Hilker [2]. The predator can survive as long as population distributions are above the threshold. This finding is no longer true for the discrete-time system as illustrated numerically in Figure 2. Specifically, the region for predator extinction includes not only those small predator populations but also large predator populations. Therefore, the predator Allee threshold in the discrete-time setting is more complicated than a strictly decreasing curve demonstrated numerically in [2].
This study concludes that cooperative hunting of predators may provide a mechanism for predator persistence when the environment is not favourable. Our result may have significant implications in biological control using predators since cooperative hunting can induce Allee effects in predators, which may cause biological control failure. Additionally, from the geometry of the basin of attraction of the predator extinction steady state, predators may go extinct even if the population sizes seem large at observation.