Cooperative hunting in a discrete predator-prey system

We propose and investigate a discrete-time predator-prey system with cooperative hunting in the predator population. The model is constructed from the classical Nicholson-Bailey host-parasitoid system with density dependent growth rate. A sufficient condition based on the model parameters for which both populations can coexist is derived, namely that the predator's maximal reproductive number exceeds one. We study existence of interior steady states and their stability in certain parameter regimes. It is shown that the system behaves asymptotically similar to the model with no cooperative hunting if the degree of cooperation is small. Large cooperative hunting, however, may promote persistence of the predator for which the predator would otherwise go extinct if there were no cooperation.


Introduction
Cooperation between individuals of social animals is frequently observed and widespread in biological systems. For example, carnivores such as wolves, wild dogs and lions often work together to capture and kill their preys [10]. Other organisms such as spiders, birds and ants also seek and attack prey collaboratively [11]. However, there are only a few mathematical models constructed to study such a biological phenomenon.
Previous research on 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. It is concluded that cooperative hunting can improve persistence of the predator but may also promote a sudden collapse of the predator. In addition, this research suggests that cooperative hunting is a mechanism for inducing Allee effects in predators.
Ever since the pioneer work of May [9], mathematical models of difference equations have played important roles in the understanding of population interactions. There are many populations with non-overlapping generations and discrete-time models are more appropriate to describe such populations. Additionally, data of ecological studies are usually collected in discrete formats. Motivated by these, the goal of this study is to propose and investigate the effects of cooperative hunting among predators upon predator-prey interactions in the discrete-time setting. Our model derivation is built on the well-known Nicholson-Bailey model with density-dependent prey growth rate. Based on the stability of the boundary equilibria, we provide a set of sufficient conditions for population coexistence, where the conditions do not depend on the cooperative hunting. We show that the system has the same asymptotic dynamics as the model of no cooperative hunting if the degree of cooperation is small. If the degree of cooperative hunting is large, then the system may support two coexisting steady states for which the predator would otherwise go extinct if there were no cooperation in this parameter regime. Numerical simulations are presented to confirm our analytical findings and to further our understanding of the predator-prey interactions.
In the following section, model derivation and persistence of the populations are presented. Section 3 provides results on the existence and the number of interior steady states. Asymptotic dynamics and local stability of the interior steady states are given in Section 4. The final section summarizes results and provides conclusions.

Model derivation and persistence of populations
Let N (t) and P (t) denote the hosts and parasitoids in generation t respectively, t = 0, 1, · · · . In the classical Nicholson-Bailey model [1], the number of encounters between hosts and parasitoids in generation t is assumed to follow the law of mass action, aN (t)P (t), where the constant a > 0 denotes searching efficiency of the parasitoids. It is also assumed that the number of encounters is distributed randomly and follows a Poisson distribution with probability p(n) = e −µ µ n n! , where n = 0, 1, 2, · · · is the number of encounters and µ is the average of encounters per host per generation. It follows that µ = aN (t)P (t) N (t) = aP (t) and thus 1 − p(0) = 1 − e −aP (t) is the probability of an individual host being parasitized in generation t since only the first encounter results in parasitism. The well known Nicholson-Bailey model is given by N (t + 1) = rN (t)e −aP (t) where all of the parameters are positive constants. Notice that the host population grows exponentially in (2.1) and the unique interior steady state is always unstable when it exists [1].
In the context of predator-prey interactions, we let x(n) and y(n) denote respectively the prey and predator populations at time n = 0, 1, 2, · · · . In the absence of cooperative hunting and by applying a similar argument as in the derivation of Nicholson-Bailey model, the probability of a prey that escaped from predation at time n is e −ay(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 the 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.
To avoid the perpetual instability of the interior steady state in (2.1), we modify the density-independent growth rate given in (2.1) by assuming that the per capita growth rate of the prey is density-dependent and is modeled by the Beverton-Holt equation. Putting these together, the interaction between prey and predator populations is described by the following difference equations: x(n + 1) = λx(n)g 0 (x(n))e −ay(n)(1+αy(n)) with nonnegative initial conditions, where λg 0 (x) = λ/(1 + kx), λ, k > 0, is the prey's per capita growth rate. The parameter α ≥ 0 denotes cooperative hunting of the predator if α > 0, and a > 0 is the searching efficiency of the predator. Further, β > 0 is the predator conversion for each prey consumed. We nondimensionalize system (2.2) by letting Ignoring the tildes, (2.2) is converted into the following system with only three parameters e −y(n)(1+αy(n)) y(n + 1) = βx(n) 1 − e −y(n)(1+αy(n)) .

(2.4)
We first observe that solutions of (2.4) remain nonnegative and are bounded for n ≥ 0. The trivial steady state E 0 = (0, 0) exists for all feasible parameters and the Jacobian matrix of (2.2) at (x, y) is given by and hence E 0 is asymptotically stable if λ < 1 and it is a saddle point with the stable manifold lying on the nonnegative y-axis if λ > 1.
Proposition 2.1 If λ ≤ 1, then E 0 = (0, 0) is globally asymptotically stable for (2.4). Proposition 2.1 implies that if the intrinsic growth rate λ of the prey population is not greater than one, then the prey population goes extinct and so does the predator population.
We assume λ > 1 for the remainder discussion so that the prey population can persist in the absence of predator. It follows that (2.2) has another boundary steady state Notice thatx can be viewed as the carrying capacity of the prey population. The Jacobian matrix of (2.4) evaluated at E 1 is Observe that 0 < 1 + λxg ′ 0 (x) < 1. Therefore E 1 is asymptotically stable if βx < 1 and it is a saddle point with its stable manifold lying on the positive x-axis if βx > 1.
Since there are only two boundary steady states for which their stability is known, we prove that system (2.4) is uniformly persistent when E 1 is unstable. That is, there exists η > 0 such that lim inf n→∞ x(n) ≥ η and lim inf n→∞ y(n) ≥ η for all solutions of (2.4) with x(0) > 0 and y(0) > 0. Our proof is based on the boundary dynamics of (2.4) using Theorem 4.1 of [7]. Theorem 2.2 Let λ > 1 and βx > 1. Then system (2.4) is uniformly persistent.
and ǫ <z is shown. We then have lim inf n→ y(n) ≥z > ǫ and obtain a contradiction. Consequently, {E 1 } is isolated in R 2 + . It is straightforward to see that the stable sets of E 0 and E 1 lie on Y and hence system (2.4) is uniformly persistent by Theorem 4.1 of [7]. Theorem 2.2 indicates that both prey and predator populations can coexist if βx > 1, where βx can be viewed as the maximal reproductive number of the predator since the prey population is stabilized at the carrying capacityx level. The parameter α plays no role in the sufficient condition for coexistence. On the other hand if βx < 1, then since E 1 = (x, 0) is asymptotically stable, the system is not uniformly persistent.
To study the effects of cooperative hunting, we need to understand the dynamics of the model when there is no cooperative hunting in the predator. The dynamics of such a model are given in Theorem 3.1 of [8] and are restated as follows.
Although it is not proved analytically in [8], it is observed that when α = 0 the unique interior steady state loses its stability via a Neimark-Sacker bifurcation as β increases.
3 Interior steady states for α > 0 In this section, we study existence and number of interior steady states. These are achieved by investigating geometry of the isoclines.
has exactly one interior steady state in case βx ≥ 1 and at most two interior steady states in case βx < 1.
Otherwise, w ′′ ≥ 0 by (3.16) and (3.13), which leads to a contradiction to (3.10) and (3.12). Using w ′′ (y 1 ) ≤ 0 and u(∞) = ∞ again, there exist 0 < y 3 < y 4 with y 1 ∈ [y 3 , y 4 ] such that Now we are ready to prove the conclusions of the theorem. Consider If y c ≤ y 0 , then w is concave on the whole interval [0, y c ] by (3.17). Since w(0) = q(0) = 0, it follows easily that system (3.11) has none or one intersection on (0, y c ) depending on whether In case y c > y 0 , w is concave on the interval [0, y 0 ] and convex on [y 0 , y c ]. We will count the number of intersection points on each interval above. Under (3.20), there is no interior intersection on (0, y 0 ] as before. There is no intersection on (y 0 , y c ) either due to q(y 0 ) > w(y 0 ) or by (3.12), q ↑ strictly on [y 0 , y c ] and w ↓ strictly on [y 0 , y c ]. (3.21) If βx > 1, then (3.19) implies that k = q ′ (0) < w ′ (0). We need to compare q(y 0 ) with w(y 0 ). By In view of (3.18), we need to compare y 4 with y c . In case y 4 < y c , we will count separately the number of intersection points on (0, y 3 ], (y 3 , y 4 ] and (y 4 , y c ). Denote them by N 1 , N 2 and N 3 respectively. First note that by (3.18) and (3.12), In fact N 2 + N 3 = 1 + 0 or 0 + 1 depending on whether q(y 4 ) ≥ w(y 4 ).
Recall that βx is the maximal reproductive number of the predator since the prey population is stabilized at its carrying capacityx. If this reproductive number exceeds one, then the predator-prey interaction can support a unique coexisting steady state. If this reproductive number is smaller than one and the degree of cooperation is also small, then Theorem 3.1 implies that the predator-prey interaction has no coexisting steady state while the interaction may support two coexisting steady states if the degree of predator cooperation is large.
to emphasize its dependence on β. For fixed y > 0, h β (y) decreases strictly to 0 as β increases to ∞. Therefore on the interval [0, y c ), We let β increase from 0. Apparently, both isoclines do not intersect for β small. At a certain β * , both isoclines become tangent to each other. The tangent point is unique. Otherwise by increasing β a little over β * , two isoclines will have four intersection points, which is contrary to Theorem 3.1. Denote the tangent point by E * = (x * , y * ). In particular, (3.26) holds. For . Two isoclines will have exactly two interior steady states in view of Theorem 3.1. Denote them by E * i = (x * i , y * i ), i = 1, 2, with 0 < y * 1 < y * < y * 2 < y c . Note that as β increases, E * 1 moves to the left and E * 2 moves to the right along the the x-isocline. When β ↑ 1/x, E * 1 becomes (x, 0), no longer an interior steady state, and lim E * 2 exists. Denote it by E * = (x e , y e ). This is consistent with Theorem 3.1 which indicates that there is exactly one interior steady state for β ≥ 1/x. The discussion is summarized as follows.
(3.27) Theorem 3.2 provides a criterion in terms of β for which the predatorprey interaction can support two coexisting steady states, where β is the predator conversion for each prey consumed. If the maximal reproductive number of the predator is smaller than one, then system (2.4) has two coexisting steady states if β is larger than the critical value β * .
This property holds as long as λ > 1 and β > 1 x for which Theorem 3.1 guarantees that there is exactly one interior steady state E * = (x * , y * ).

Stability and dynamics of the model
To study asymptotic dynamics of (2.4) for λ > 1 and α > 0, we separate our discussion into two cases: 2α ≤ 1, and 2α > 1. Since local stability can be determined by the Jury conditions, we first provide a result on the determinant of the Jacobian matrix J at an interior steady state E = (x, y). Using (3.1) and (3.6) we may rewrite J as follow Recall ♦ = y(1 + αy) by (3.2). Following Theorem 3.2, we will let E vary along the x-isocline and thus J in (4.1) becomes a function of y.  Under this assumption, Theorem 3.1 implies that (2.4) has a unique interior steady state E = (x, y) if βx > 1, and there is no interior steady state if βx < 1. Our analysis is more complete in this parameter regime. In particular, E 1 = (x, 0) is globally asymptotically stable if βx < 1, which is similar to the system with α = 0 as illustrated in Proposition 2.3.
Proof. Under the given assumptions, E 1 is locally asymptotically stable and there is no interior steady state. It is straightforward to verify that 1 − e −y(1+αy) < y for y > 0.
Theorem 4.2 states that if the intrinsic growth rate λ of the prey population is greater than one, the degree of predator's cooperative hunting is small, 2α ≤ 1, and the maximal reproductive number of the predator is less than one, βx < 1, then the predators will go extinct and the prey population will stabilize at its carrying capacityx.
Since magnitude of the interior steady state is monotone with respect to β by (4.4) and β is the prey conversion to predator, we use β as the bifurcation parameter. Let β d > 0 be the corresponding β value of y d given in Lemma 4.1. The following result follows from Lemma 4.1 and (4.8). Let C * be defined in (4.14). To study the Neimark-Sacker bifurcation, we let the unique interior steady state be denoted by E * = (x * , y * ). Theorem 4.3 Let λ > 1, 2α ≤ 1 and βx > 1. Then the unique interior steady state E * = (x * , y * ) is asymptotically stable if β < β d and a repeller if β > β d . Moreover, E * undergoes a Neimark-Sacker bifurcation at β = β d . The bifurcation is supercritical if C * > 0 and the bifurcation is subcritical if C * < 0.
If C * > 0, then the system has an attracting closed invariant circle for β > β d and near β d . If C * < 0, then the system has an unstable closed invariant circle for β < β d and near β d .
It is not easy to determine analytically whether the bifurcation is supercritical or subcritical since C * cannot be computed analytically. Numerical investigation does indicate that the bifurcation is supercritical so that the model has an attracting invariant closed curve when β > β d and near β d .
We conclude from Theorems 4.2, 4.3 and Proposition 2.3 that cooperative hunting does not affect dynamical interactions of the prey and predator if the degree of cooperative hunting is small, 2α ≤ 1.

Dynamics of the model when 2α > 1
Recall Theorem 3.1 indicates that (2.4) has a unique interior steady state if βx > 1, (2.4) has no interior steady state if βx < 1 and 1 < 2α ≤ 3λ − 1 λ − 1 , and there are either zero, one or two interior steady states if βx < 1 and We shall determine stability of an interior steady state when it exists. In the case when (2.4) has no interior steady state, we suspect that E 1 is globally asymptotically stable. The following result provides a restriction on the parameter α for which E 1 is a global attractor.
Theorem 4.4 provides a sufficient condition for which the predators go extinct when the degree of cooperation is large. That is, large cooperation among predators drive the predators to extinction under the condition given by the theorem. We now study local stability of an interior steady state E = (x, y) when 1 < 2α. Recall from Lemma 4.1 and (4.6) that in order to determine local stability at E, we have to study V (y) defined in (4.7). Note that y c is defined in (3.7). The result is as follows.
Since d ds  19), we obtain G ′ > 0, G > 0, F ′ > 0, F > 0 and at last V > 0 on (0,ỹ 1 ). Together with (4.18), that V (y) > 0 It remains to consider the case 2α > 3λ − 1 λ − 1 . By (4.24) and (4.25), L(y) < 0 on (0, δ) and L(y) > 0 on (δ,ỹ 1 ) for some δ ∈ (0,ỹ 1 ). The same holds for G ′ by (4.23). Since   When the maximal reproductive number of the predators exceeds one and the degree of cooperation α is neither too small nor too large, Theorem 4.6 implies that the predator-prey interaction can support a unique coexisting steady state. Consequently, both populations can coexist indefinitely as a steady state if the predator conversion β is smaller than a critical value β d . Otherwise, coexistence of both populations may be more complicated if β is larger than β d . x .
See Fig 1(a). We present some numerical investigations for the asymptotic dynamics of the system. Using the parameter values λ = 5 and α = 1/2.1, then a unique interior steady state exists if β > 0.25 and a Neimark-Sacker bifurcation occurs when β is close to 0.6. Figure 3(a) presents an invariant closed curve for β = 0.609. We next increase α to 3/2.1, Figure 3(b) provides an invariant closed curve for the case when βx > 1 and 2α < 3λ − 1 λ − 1 . The initial conditions are chosen near the unstable unique interior steady state for both plots.
We next study the scenario when 2α > 3λ − 1 λ − 1 , βx < 1 and the system has two interior steady states. The parameter values used are λ = 5, β = 4.2/20, and α = 20/2.1. We choose an initial condition (2.3, 0.2) which is close to E * 2 . The solution converges to a closed invariant circle. If initial condition (3.9, 0.1) is used, then the solution converges to the boundary steady state E 1 = (x, 0) = (4, 0) as shown in Fig 3(c). Notice that in this parameter regime, both E * 1 and E * 2 are unstable. We wish to demonstrate the stability of E * 2 and thus we decrease β to 3.76/20 while keep all other parameter values the same. Then the two isoclines have two positive intersections which results in two interior steady states. We use the same initial conditions as in Fig 3(c). In this circumstance, one solution converges to the stable interior steady state while the other solution converges to the boundary steady state E 1 . See Fig 3(d). Therefore, bistability occurs and the predator may survive depending on initial conditions while the predator would go extinct if there is no cooperative hunting.
We summarize conditions for the existence of interior steady states in Table 1 and a list of notations is given in Table 2. Table 1 Existence of interior steady states Parameter regime parameter regime number of interior steady states βx > 1 α ≥ 0 1 depends on the location of y relative to y d and y t . In addition, there are also y * and y e involved. See (3.27), Lemmas 4.1 and 4.5. It is hard to compare the order of these quantities theoretically. We postpone our investigation to a future study.

Summary and conclusions
Mathematical models of predator-prey interactions are interesting dynamical systems. There are many populations in nature with non-overlapping generations. Consequently, continuous-time models are not appropriate to describe such population interactions and discrete-time systems can be used to explore such populations. Cooperation among individuals of the same predator species is frequently observed in nature and it can change dynamical interactions of biological systems [10,11]. Motivated by the recent research of Alves and Hilker [2] on continuous-time models of predator-prey interactions with cooperative hunting in predators, we propose and investigate a parallel discrete-time system. The model derivation is based on the classical Nicholson-Bailey system but with density-dependent growth rate in the prey population. Similar to [2], cooperative hunting of the predator is modeled via the attack rate of the predator. Due to this cooperation, the probability of an individual prey escaped from being preyed upon is decreased. In order to investigate the effects of cooperative hunting, the dynamics of the system with no cooperation among predators are summarized first.
Comparing the system of cooperation with that of no cooperation, several similar dynamical results are obtained. Indeed, both populations go extinct if the intrinsic growth rate of prey is smaller than one while both populations can coexist under the same sufficient conditions, namely that the prey's intrinsic growth rate and the predator's maximal reproductive number are both greater than. Further, it is proven that asymptotic dynamics of the model are similar to the system with no cooperation if 2α ≤ 3λ − 1 λ − 1 , where α is the degree of cooperation. Consequently, if the degree of cooperation α is small, then cooperative hunting does not change dynamical interactions between the prey and predators. On the other hand, if the degree of cooperation α is large, i.e., 2α > 3λ − 1 λ − 1 , then cooperative hunting becomes critical for the survival of the predator in the case that βx < 1. The lumped parameter βx can be interpreted as the maximal reproductive number of the predator. Without cooperation, the predator population goes extinct if this reproductive number is smaller than one. See Proposition 2.3. With cooperative hunting, the predator-prey interactions may support two interior steady states when this reproductive number is less than one as illustrated in Theorem 3.1. As a result, the predator and prey may coexist even if the maximal reproductive number of predator is smaller than one. Therefore, cooperation between predators can promote survival of the predator which would otherwise go extinct in the absence of this mechanism. Comparing our results with those of the continuous-time model studied by Alves and Hilker [2], first notice that in the absence of predator's cooperation the unique interior steady state in [2] is globally asymptotically stable whenever it exists. This is not true for our system since the unique interior steady state can undergo a Neimark-Sacker bifurcation for system (2.4) when predators do not engage in cooperation. However, using the concept of uniform persistence we prove that both populations can coexist indefinitely as long as the maximal reproductive of predators exceed one independent of whether predators cooperate or not. This coexistence is not proved in [2] when predators engage in hunting cooperation. On the other hand, the number of interior steady states for both of the continuous and discrete-time models is the same. In particular, both systems can have two coexisting steady state if the degree of cooperation is large and the predator's maximal reproductive number of predator is less than one. In our study, however, we are able to quantify this degree of cooperation explicitly in terms of the prey's intrinsic growth rate. Furthermore, for small degree of cooperation, the asymptotic dynamics of the continuous-time model are the same as the model with no cooperation. The discrete-time model proposed in this study also possess this property, namely that the asymptotic dynamics of the system with small magnitude of predator cooperation behave asymptotically the same as the model of no cooperation. (c) (d) Figure 3: Invariant closed curves and bistability are presented. In (a), 2α ≤ 1 and βx > 1. In (b), 2α > 1 and βx > 1. In (c) where 2α > 1 and βx < 1, the system has two interior steady states and two attractors are shown using two different initial conditions. One solution converges to the boundary steady state E 1 = (4, 0) and the other converges to the closed invariant circle. We decrease β to β = 3.76/20 so that system (2.4) still has two interior steady states, where one is unstable and the other is asymptotically stable. Using the same initial conditions as in (c), one solution converges to E 1 while the other converges to the interior steady state as shown in (d).