Hopf bifurcation and Turing instability in a diffusive predator-prey model with hunting cooperation

1 Introduction

The predator is a population that eats another population (known as prey). The prey is the food source of the predator, and the predator dies if it does not obtain food. Modeling the interaction between preys and predators has become one of the important methods to reveal their evolution law. The mathematical model describing this interaction was first proposed to explain the oscillatory levels of certain fish catches in the Adriatic by Lotka-Volterra [1,2]. Since then, based on different biological meanings, statistical data, and experiments, a lot of predator-prey models have been proposed and studied and have become an important tool for mathematicians and biologists to learn more about nature and better protect the environment.

Cooperative hunting is a social behavior intrinsically present in many predator populations. This is a common interest among predators for successfully capturing preys compared to individual efforts. Different predators employ different strategies based on their moving abilities, social skills, communication capabilities, cognitions, prey types, and availabilities. Bailey et al. [3] reviewed different strategies taken by predators during cooperative hunting of over 40 species, including carnivores [4] (lions [5], wolves [6], African wild dogs [7], and chimpanzees [8]), birds [9], ants [10], and spiders [11].

Although there are many cooperative predators in nature, there are only a few mathematical models incorporating such a biological mechanism. To the best of our knowledge, an earlier predator-prey model with hunting cooperation has been proposed by Berec [12]. Berec found that hunting cooperation has a destabilizing effect on predator-prey dynamics. Based on the research of Berec, Alves and Hilker [13] extended the classical Lotka-Volterra model by including hunting cooperation as follows:

where N and P are prey and predator densities, respectively, r is the intrinsic growth rate of the prey, K is the carrying capacity of the prey, e is the conversion efficiency, m is the mortality rate of the predator, β is the rate of predation without hunting cooperative, and α is the cooperative coefficient. r , K , e , m , β are positive constants, and, α is a non-negative constant. If α = 0 , model (1.1) becomes a classical Lotka-Volterra model with logistic growth of the prey. For simplicity, consider dimensionless variables with the following scaling:

model (1.1) is rewritten as follows:

with

Alves and Hilker [13] investigated the stability of equilibrium in phase plane and bifurcation diagrams of model (1.2) via numerical simulations. It is shown that for some σ and a Hopf or Bogdanov-Takens bifurcations occur. Recently, Wu and Zhao [14] studied the stability of equilibrium and the existence of Hopf bifurcations for model (1.2) via mathematical analysis.

As far as we know, reaction-diffusion equations have been widely employed to study the behavior of biological phenomena in nature. In particular, Turing [15] showed that a system of coupled reaction-diffusion equations can be used to describe patterns and forms in biological systems. Turing’s theory shows that diffusion could destabilize an otherwise stable equilibrium of the reaction-diffusion system and lead to nonuniform spatial patterns. This kind of instability is usually called Turing instability. Naturally, we then consider the reaction-diffusion model corresponding to model (1.2):

where Ω is a bounded domain in R n with smooth boundary ∂ Ω and ν is the outward unit normal on ∂ Ω . d 1 and d 2 are the diffusion coefficients of u and v , respectively. The initial data u 0 ( x ) , v 0 ( x ) are nonnegative smooth functions, which are not identically zero. Model (1.3) has been extensively investigated by many researchers and some interesting results have also been obtained. For example, Wu and Zhao [14] considered the stability of equilibrium and, the existence of Hopf bifurcations and Turing instability. Capone et al. in [16] proved the definitive boundedness of solutions via the existence of positive invariants and attractive sets and obtained the Turing instability. Singh and Dubey [17] investigated meticulousness in the spatial dynamics of the predator-prey model with hunting cooperation via diffusion-driven instability.

However, we note that the results in [14,16] provide no further information on the direction and stability of the Hopf bifurcation in models (1.2) and (1.3). Therefore, it may be interesting to investigate such cases. We point out that this topic has been largely discussed in the literature for the predator-prey model. Based on the aforementioned reasons, in the present article, we consider the following predator-prey model with hunting cooperation

Using the center manifold and the normal form theory, we derive conditions for determining the direction of Hopf bifurcation and the stability of the bifurcating periodic solution. Based on the above discussions, we further investigate the following reaction-diffusion model:

It is worth mentioning that we can immediately reach the direction and stability of the Hopf bifurcation in model (1.5) without any tedious calculations.

For other research studies of the predator-prey model with hunting cooperation, see [18] for cross-diffusion, see [19] for time delay, see [20,21] for discrete time; in particular, Pati et al. [20] explored bifurcation patterns, novel organized structures, and chaos in the bi-parameter space of a discrete-time predator-prey model with cooperative hunting.

This article is organized as follows. In Section 2, we investigate the existence, direction, and stability of the Hopf bifurcation in model (1.4) by the Poincaré-Andronov-Hopf bifurcation theorem. In Section 3, we first consider the Turing instability of the reaction-diffusion model (1.5) when the spatial domain is a bounded interval, and then we study the existence and direction of Hopf bifurcation and the stability of the bifurcating periodic solution. Moreover, numerical simulations are presented to verify and illustrate the theoretical results above.

2 Dynamics of the local model

In this section, we mainly discuss the existence, direction, and stability of the Hopf bifurcation in model (1.4).

It is easy to verify that (1.4) has a unique positive equilibrium ( u c , v c ) = 1 a σ , a σ − 1 a if and only if a > 1 σ . The Jacobian matrix of (1.4) at ( u c , v c ) is as follows:

and the corresponding characteristic equation is

where Tr ( a ) = a σ − 1 − σ a σ and Det ( a ) = 2 ( a σ − 1 ) 2 a > 0 . By the standard linearization method, we can obtain that if a < ( 1 + σ ) 2 σ holds, then the equilibrium ( u c , v c ) of (1.4) is locally stable.

Next, we analyze the Hopf bifurcation occurring at ( u c , v c ) by choosing a as the bifurcation parameter. Denote a 0 = ( 1 + σ ) 2 σ , then when a = a 0 , the Jacobian matrix J ( a ) has a pair of purely imaginary eigenvalues λ = ± i ω 0 , where ω 0 = 2 σ 3 ( 1 + σ ) 2 . Let λ = β ( a ) ± i ω ( a ) be the roots of λ 2 − Tr ( a ) λ + Det ( a ) = 0 , then

Furthermore, we can verify

By the Poincaré-Andronov-Hopf bifurcation theorem (see [22, Theorem 3.1.3]), we know that (1.4) undergoes a Hopf bifurcation at ( u c , v c ) as a passes through a 0 .

However, the detailed nature of the Hopf bifurcation needs further analysis of the normal form (see [23,24,25]) for model (1.4). To this end, we translate the equilibrium ( u c , v c ) to the origin by the translation u ˜ = u − u c , v ˜ = v − v c . For the sake of convenience, we still denote u ˜ and v ˜ by u and v , respectively. Thus, the local model (1.4) is transformed into

Rewrite (2.2) as

where

Choose ξ = ( i ω 0 − a u c v c , ( 1 + a v c ) v c ) T as one of the eigenvectors of J ( a ) corresponding to i ω 0 . Let

then

The transformation u v = P x y changes (2.3) into

where

So, the stability of Hopf bifurcation for model (1.4) at ( u c , v c ) is determined by the sign of the following quantity:

where all the partial derivatives are evaluated at the bifurcation point ( x , y , a ) = ( 0 , 0 , a 0 ) .

After calculation, we have

Then

By [26, Chapter 3] and β ′ ( a 0 ) > 0 , δ ( a 0 ) < 0 , we summarize our results as follows.

3 Dynamics of the diffusion model

For simplicity, in this section, we take the spatial domain Ω of model (1.5) as the one-dimensional interval Ω = ( 0 , π ) .

4 Numerical simulations

In this section, we present some numerical simulations to illustrate our theoretical analysis.

The local model (1.4) involves two parameters a , σ . First, choosing parameter σ = 1 , we have the critical point a 0 = 4 . By Theorem 2.1(i), we know the equilibrium 1 3 , 3 − 1 3 is locally stable when a = 3 < a 0 (see Figure 1). From Theorem 2.1(ii), a Hopf bifurcation occurs at a = 4.1 > a 0 , the direction of the bifurcation is supercritical and the bifurcating periodic solutions are asymptotically stable (Figure 2).

Figure 1

Phase portraits of model (1.4) with parameter σ = 1 , the equilibrium 1 3 , 3 − 1 3 is locally stable, where a = 3 < a 0 .

Figure 2

Phase portraits of model (1.4) with parameter σ = 1 , the bifurcating periodic solution is stable, where a = 4.1 > a 0 .

Take parameters as a = 2 , σ = 1 , d 1 = 1 , d 2 = 0.5 , then a < a 0 = 4 and (H2) are satisfied; the equilibrium ( u c , v c ) = 1 2 , 2 − 1 2 is still stable under diffusive effects (Figure 3). Initial values are chosen as u ( x , 0 ) = 0.75 + 0.001 rand ( 1 ) , v ( x , 0 ) = 0.2 + 0.001 rand ( 1 ) .

Figure 3

Phase portraits of (1.5) with parameters a = 2 , σ = 1 , d 1 = 1 , d 2 = 0.5 , the equilibrium 1 2 , 2 − 1 2 is stable.

If the parameters are taken as a = 2 , σ = 1 , d 1 = 2 , d 2 = 0.2 , then a < a 0 = 4 and (H3) are satisfied, the equilibrium ( u c , v c ) = 1 2 , 2 − 1 2 is still stable under diffusive effects (Figure 4). Initial values are chosen as u ( x , 0 ) = 0.75 + 0.001 rand ( 1 ) , v ( x , 0 ) = 0.2 + 0.001 rand ( 1 ) .

Figure 4

Phase portraits of (1.5) with parameters a = 2 , σ = 1 , d 1 = 2 , d 2 = 0.2 , the equilibrium ( 1 2 , 2 − 1 2 ) is stable.

However, when a = 2 , σ = 1 , d 1 = 10 , d 2 = 0.2 , then a < a 0 = 4 and (H1) are satisfied; the equilibrium ( u c , v c ) = 1 2 , 2 − 1 2 becomes unstable under diffusive effects (Figure 5). Initial values are chosen as u ( x , 0 ) = 0.75 + 0.001 rand ( 1 ) , v ( x , 0 ) = 0.2 + 0.001 rand ( 1 ) .

Figure 5

Phase portraits of (1.5) with parameters a = 2 , σ = 1 , d 1 = 10 , d 2 = 0.2 , the equilibrium 1 2 , 2 − 1 2 becomes unstable.

Taking a = 4.1 , σ = 1 , d 1 = 1 , d 2 = 0.5 , then a > a 0 = 4 and (H2) are satisfied. In this case, theoretical analysis showed that the limit cycle from Hopf bifurcation is stable for (1.5). For numerical results, see Figure 6. Initial values are taken as u ( x , 0 ) = 0.75 + 0.2 rand ( 1 ) , v ( x , 0 ) = 0.2 + 0.2 rand ( 1 ) , close to the limit cycle.

Figure 6

Phase portraits of (1.5) with parameters a = 4.1 , σ = 1 , d 1 = 1 , d 2 = 0.5 . Bifurcated homogeneous periodic solution is stable.

The similar phenomenon occurs when parameters are taken as a = 4.1 , σ = 1 , d 1 = 2 , d 2 = 0.2 , then a > a 0 = 4 and (H3) are satisfied. In this case, the limit cycle from Hopf bifurcation is stable for (1.5) (Figure 7). Initial values are taken as u ( x , 0 ) = 0.75 + 0.2 rand ( 1 ) , v ( x , 0 ) = 0.2 + 0.2 rand ( 1 ) , close to the limit cycle.

Figure 7

Phase portraits of (1.5) with parameters a = 4.1 , σ = 1 , d 1 = 2 , d 2 = 0.2 . Bifurcated homogeneous periodic solution is stable.

When parameters are taken as a = 4.1 , σ = 1 , d 1 = 10 , d 2 = 0.2 , then a > a 0 = 4 and (H1) are satisfied. In this case, the limit cycle from Hopf bifurcation is unstable for (1.5). Numerical results are shown in Figure 8. Initial values are taken as u ( x , 0 ) = 0.75 + 0.2 rand ( 1 ) , v ( x , 0 ) = 0.2 + 0.2 rand ( 1 ) .

Figure 8

Phase portraits of (1.5) with parameters a = 4.1 , σ = 1 , d 1 = 10 , d 2 = 0.2 . Bifurcated homogeneous periodic solution is unstable.

5 Conclusion

In this article, we study the diffusive predator-prey model (1.5) with hunting cooperation under Neumann boundary conditions. The predator cooperation in hunting rate a plays a key role in determining the dynamics of the model. We provide detailed analyses of the Hopf Bifurcation and Turing instability in model (1.5) via three possible mechanisms: (i) For the local model of (1.5), we analyzed the stability of the equilibrium and derived conditions for determining the direction of Hopf bifurcation and the stability of the bifurcating periodic solution by the center manifold and the normal form theory (see Theorem 2.1). Our result showed that hunting cooperation could be beneficial to the predator population. (ii) We studied the Turing instability of the reaction-diffusion model (1.5) when the spatial domain is a bounded interval (see Theorem 3.1). Our result showed that hunting cooperation might be one of the determining factors in producing Turing patterns. In fact, we can find that Turing instability never occurs in model (1.5) without hunting cooperative. (iii) Combining the analysis of (i) and (ii), we immediately reach the direction and stability of the Hopf bifurcation for model (1.5) without any tedious calculations (see Theorem 3.2). Moreover, numerical simulations are also carried out to illustrate theoretical analysis, from which the theoretical results are verified. More interesting and complex behavior of such models will further be explored.