Complex Dynamical Behavior of a Predator-Prey System with Group Defense

Abstract

A diffusive predator-prey system with prey refuge is studied analytically and numerically. The Turing bifurcation is analyzed in detail, which in turn provides a theoretical basis for the numerical simulation. The influence of prey refuge and group defense on the equilibrium density and patterns of species under the condition of Turing instability is explored by numerical simulations, and this shows that the prey refuge and group defense have an important effect on the equilibrium density and patterns of species. Moreover, it can be obtained that the distributions of species are more sensitive to group defense than prey refuge. These results are expected to be of significance in exploration for the spatiotemporal dynamics of ecosystems.

1. Introduction

Some observations [1–3] indicate that there exists a biological phenomenon called group defense in predator-prey interactions. Namely, the prey species exhibit group defense whereby predation is decreased or prevented since the preys could better defend or disguise themselves when their numbers are large enough. This biological phenomenon has been widely incorporated into many predator-prey systems, and it is shown that the group defense has a significant effect on the dynamics of these systems [4–6]. During the past decades, predator-prey system has been relatively well-studied because the predator-prey interaction is universal throughout nature [7–10]. Recently, a class of the so-called semi-ratio-dependent predator-prey systems based on the Leslie-Gower models have been extensively investigated, an important trait of which is that the carrying capacity of predator is proportional to the number of preys [11–14]. However, the Leslie-Gower model incorporating group defense is rarely investigated. So a modified Leslie-Gower predator-prey system with group defense as suggested by Sokol and Howell [15] is introduced, which can be expressed as follows: 𝑑𝐻𝑑𝑡=𝑟𝐻−𝑎𝐻²−𝑐𝑃𝐻𝑏+𝐻²,𝑑𝑃𝑑𝑡=𝑠𝑃−𝑑𝑃²,𝑒+𝐻 (1) where 𝐻 and 𝑃 correspond to the densities of prey and predator species at time 𝑡.

On the other hand, many studies [16–18] have indicated that the prey refuge has a positive effect on the dynamics of predator-prey models. In particular, mite predator-prey interactions [19] often exhibit spatial refuges that can afford the prey with some degree of protection from predation and reduce the likelihood of extinction due to predation. Thus, system (1) can be extended by including a protective refuge 𝑚𝐻 for prey, where 𝑚∈[0,1) is constant and (1−𝑚)𝐻 is available to be the predator. More precisely speaking, the spatial factor should be included into system (1) since all organisms inhabit the three-dimensional space. Moreover, a large number of studied models incorporating the spatial factor have shown very rich and complex dynamics. Accordingly, system (1) can be extended as follows: 𝜕𝐻𝜕𝑡=𝑑₁Δ𝐻+𝑟𝐻−𝑎𝐻²−𝑐(1−𝑚)𝑃𝐻𝑏+(1−𝑚)²𝐻²,𝜕𝑃𝜕𝑡=𝑑₂Δ𝑃+𝑠𝑃−𝑑𝑃²()𝐻,𝑒+1−𝑚 (2) where 𝐻(𝑥,𝑡) and 𝑃(𝑥,𝑡) correspond to the densities of prey and predator at time 𝑡 and position 𝑥. Δ is the Laplacian operator. Δ𝐻 and Δ𝑃 denote the populations moving randomly in the environment. 𝑑₁ and 𝑑₂ are diffusivities. 𝑟 and 𝑠 are the intrinsic growth rates of 𝐻 and 𝑃, respectively. 𝑎 measures the strength of intraspecific competition among the individual of prey. 𝑏 is interpreted as the half-saturation constant. 𝑐 stands for the predation rate. 𝑑/𝑠 measures the ration of prey to support one predator. 𝑒 reflects the extent to which the environment provides protection to prey. All the parameters of system (2) are considered to be positive.

This paper is arranged as follows. In the next section, mathematical analysis is given to obtain the critical wavenumber. In Section 3, on the basis of the critical wavenumber, numerical simulations are performed to explore the effect of prey refuge and group defense on the distributions of species. Furthermore, we also analyzed numerically the impact of prey refuge and group defense on the interior equilibrium densities of species. Finally, we end with conclusions.

2. Mathematical Analysis

Assume that system (2) satisfies the following initial conditions: 𝐻(0,𝑥)=𝐻₀(𝑥)>0,𝑃(0,𝑥)=𝑃₀(𝑥)>0,𝑥∈Ω, (3) and zero-flux boundary conditions: 𝜕𝑃=𝜕𝑛𝜕𝐻𝜕𝑛=0,𝑥∈𝜕Ω,𝑡>0. (4)

Here, Ω⊂𝑅² is bounded, the boundary 𝜕Ω is smooth, and 𝜕/𝜕𝑛 denotes the differentiation in the direction of the outward unit normal vector to 𝜕Ω. The zero-flux conditions imply that there are no fluxes of species through the boundary; that is, system (2) has no input from the outside.

It is interesting to note that the local dynamics of system (2) are the basis of better understanding pattern formations under the condition of Turing instability. Thus, without diffusion, the local system of (2) is 𝑑𝐻𝑑𝑡=𝑟𝐻−𝑎𝐻²−𝑐(1−𝑚)𝑃𝐻𝑏+(1−𝑚)²𝐻²=𝐻𝜑(𝐻,𝑃),𝑑𝑃𝑑𝑡=𝑠𝑃−𝑑𝑃²()𝐻().𝑒+1−𝑚=𝑃𝜙𝐻,𝑃 (5)

From the biological angle, the steady state 𝐸_∗=(𝐻_∗,𝑃_∗) of system (2) or (5) is most significant, which implies coexistence of both predator and prey, and (𝐻_∗,𝑃_∗)=(𝐻_∗,(𝑠/𝑑)[𝑒+(1−𝑚)𝐻_∗]) is the positive solution of 𝜑(𝐻,𝑃)=0, 𝜙(𝐻,𝑃)=0, where 𝐻_∗ is the positive root of 𝐻³−𝑟𝑎𝐻²+𝑎𝑏𝑑+𝑐𝑠(1−𝑚)²𝑎𝑑(1−𝑚)²()𝐻−𝑟𝑏𝑑−𝑐𝑒𝑠1−𝑚𝑎𝑑(1−𝑚)²=0. (6) Based on the fundamental theorem of Algebra, the cubic equation (6) has at least one real root. Then the sufficient conditions are necessary to ensure at least one positive real root of the cubic equation (6). From the first equation of the system (5), it is easy to know that 0<𝐻_∗<𝑟/𝑎. Assuming that 𝑓(𝐻)=𝐻³−(𝑟/𝑎)𝐻²+((𝑎𝑏𝑑+𝑐𝑠(1−𝑚)²)/𝑎𝑑(1−𝑚)²)𝐻−((𝑟𝑏𝑑−𝑐𝑒𝑠(1−𝑚))/𝑎𝑑(1−𝑚)²), the function 𝑓(𝐻) is continuous in [0,𝑟/𝑎], 𝑓(0)=−(𝑟𝑏𝑑−𝑐𝑒𝑠(1−𝑚))/𝑎𝑑(1−𝑚)² and 𝑓(𝑟/𝑎)=(𝑐𝑠𝑟(1−𝑚)+𝑐𝑠𝑎𝑒)/𝑎²𝑑(1−𝑚)². From the Zero Theorem, there must exist at least one positive real number 𝐻_∗∈(0,𝑟/𝑎) such that 𝑓(𝐻_∗)=0 if 𝑓(0)𝑓(𝑟/𝑎)<0; that is to say 𝑟𝑏𝑑−𝑐𝑒𝑠(1−𝑚)>0.

On the other hand, the stability behavior of the homogeneous state 𝐸_∗ with small perturbations is of importance. Thus, let us linearize system (5) around 𝐸_∗. This yields the following system: 𝑑𝐗𝑑𝑡=𝐉𝐗, (7) where 𝐗=(𝐻,𝑃)^𝑇 and 𝐉=𝜕𝜑𝜕𝐻𝜕𝜑𝜕𝑃𝜕𝜙𝜕𝐻𝜕𝜙𝜕𝑃_𝐸∗=−𝑎𝐻_∗+2𝑐(1−𝑚)³𝑃_∗𝐻_∗²𝑏+(1−𝑚)²𝐻²_∗²−𝑐(1−𝑚)𝐻_∗𝑏+(1−𝑚)²𝐻²_∗𝑠²(1−𝑚)𝑑=.−𝑠𝐴𝐵𝑅𝑄 (8) Then, the characteristic equation of system (7) is 𝜆²−tr(𝐉)𝜆+det(𝐉)=0, (9) where tr(𝐉)=𝐴+𝑄 and det(𝐉)=𝐴𝑄−𝐵𝑅.

It is obvious that if tr(𝐉)<0 and det(𝐉)>0, system (7) has two distinct eigenvalues with negative real parts. Thus, 𝐸_∗ is locally asymptotically stable for system (5). Assuming that 𝑟=0.6, 𝑠=0.49, 𝑎=0.28, 𝑏=0.3, 𝑐=0.57, 𝑑=0.7, 𝑒=0.2, and 𝑚=0.15, we can get that 𝑟𝑏𝑑−𝑐𝑒𝑠(1−𝑚)=0.0785>0, and thus there exists 𝐸_∗=(0.6167,0.5069). Further, we can obtain that tr(𝐉)=−0.2540<0 and det(𝐉)=0.0359>0. Therefore, 𝐸_∗ is locally asymptotically stable for system (5) shown in Figure 1.

Figure 1 

Phase portrait of system (5) corresponding to different initial values which shows that 𝐸∗ is locally asymptotically stable.

In order to investigate the temporal stability of the uniform state 𝐸_∗, the nonuniform perturbations for any solution of system (2) are introduced as follows: 𝐻(𝑥,𝑡)=𝐻_∗)),𝑃()+𝜀exp((𝐤𝑥𝑖+𝜆𝑡𝑥,𝑡=𝑃_∗+𝜏exp((𝐤𝑥)𝑖+𝜆𝑡), (10) where 𝜀, 𝜏 are small enough, 𝐤 is the wave number, and 𝜆 is the growth rate of fluctuation in time. Substituting (10) into system (2) and picking up the linear terms, then there is 𝜕𝐘𝜕𝑡=𝐃Δ𝐘+𝐉𝐘, (11) where 𝐘=(𝐻,𝑃)^𝑇 and 𝐃=diag(𝑑₁,𝑑₂).

Therefore, the characteristic equation of (11) is 𝜆²−tr_𝐤𝜆+Δ_𝐤=0, (12) where tr_𝐤=tr(𝐉)−𝑑₁+𝑑₂𝐤²,Δ_𝐤=det(𝐉)−𝐤²𝐴𝑑₂+𝑅𝑑₁+𝐤²²𝑑₁𝑑₂. (13)(14) Thus, the roots of (12) are 𝜆_1,2(𝐤)=tr_𝐤±tr_𝐤²−4Δ_𝐤2. (15)

The Turing instability implies that the homogeneous state is stable in the local system and will become unstable due to diffusion of populations. From (7), the interior equilibrium point 𝐸_∗ of system (5) is stable if tr(𝐉)<0, det(𝐉)>0. According to (13), we can obtain tr_𝐤<tr(𝐉). Thus, the stability of the interior equilibrium point corresponding to system (2) is only determined by the sign of Δ_𝐤. From (14), it is clear to show that Δ_𝑘<0 for 𝜅²₁<𝐤²<𝜅²₂, where 𝜅²_1,2=(𝐴𝑑₂+𝑅𝑑₁∓(𝐴𝑑₂+𝑅𝑑₁)²−4𝑑₁𝑑₂(𝐴𝑅−𝐵𝑄))/2𝑑₁𝑑₂. Furthermore, if 𝜅₁ and 𝜅₂ have positive values, the range of instability for the steady state 𝐸_∗ is limited, which is called the Turing space.

To well see the Turing space, the dispersion relation corresponding to several values of the biological parameter 𝑚 is plotted while the other parameter values are 𝑟=0.6, 𝑠=0.49, 𝑎=0.28, 𝑏=0.3, 𝑐=0.57, 𝑑=0.7, 𝑒=0.2, 𝑑₁=0.08, and 𝑑₂=0.12 in Figure 2. It can be seen from Figure 2 that there exists an interval of 𝑚 such that when 𝑚 lies in this interval, the available modes occur. It should be stressed that when 𝑚 is beyond this interval, the Turing instability does not happen and instead the stable state is excited. It is interesting to note that the upper and lower bounds of this interval are easy to be found. From Figure 2, it is clear to show that when 𝑚<0.01 or 𝑚>0.27, the Turing instability is no longer excited.

Figure 2 

Variation of dispersion relation of system (2) around the interior equilibrium point for different 𝑚. The red line corresponds to 𝑚=0.01, the green line corresponds to 𝑚=0.12, and the blue corresponds to 𝑚=0.27.

3. Simulation Analysis

For better understanding the impact of a prey refuge on the dynamical behavior of system (2), the values of the equilibrium density are first plotted as the functions of the prey refuge factor 𝑚. It can be seen from Figure 3 that the prey refuge has a significant effect on the densities of prey and predator at the uniform steady state of system (2). It is shown in Figure 3(a) that when 𝑚 becomes large, the prey density of the steady state will be increased. Moreover, it is interesting to be observed that there exists a quasilinear function between the prey density of the steady sate of system (2) and the parameter 𝑚 of prey refuge. Figure 3(b) clearly shows that when 𝑚 is large enough, the density of predator at steady state will be declined as increase of 𝑚. It is to note from Figure 3(b) that there is a critical value 𝑚_𝑐 such that when 𝑚<𝑚_𝑐, the predator density 𝑃_∗ will be always increase and when 𝑚 is beyond the critical values, 𝑃_∗ will be dramatically decreased. Obviously it is a nonlinear relationship between the predator density of the equilibrium state and prey refuge parameter. This implies that the effect of prey refuge on the dynamical behavior of system (2) is complex and profound.

(a)

(b)

Figure 3 

(a) The prey density. (b) The predator density. Both are the functions of 𝑚 at steady state. The other parameter values are 𝑟=0.6, 𝑠=0.49, 𝑎=0.28, 𝑏=0.3, 𝑐=0.57, 𝑑=0.7, 𝑒=0.2, 𝑑1=0.08, and 𝑑2=0.12.

From Figure 3, it can be only obtained the effect of prey refuge on the densities of prey and predator at steady state. However, it does not reflect the spatial information of species. Thus, to well investigate effect of prey refuge on the spatiotemporal dynamics of system (2), the spatial distributions of prey are explored by numerical simulations. For this purpose, system (2) with zero-flux boundary conditions is employed on a two-dimensional domain. Furthermore, the continuous problem determined by system (2) is solved in this domain with 200 × 200 lattice sites and the step between each lattice point is 𝜔=0.25. For the time evolution of system (2), the Euler method is chosen to approximate the value of density at the next time step on the basis of the density at the previous time step. Here, the time step is 𝜌=0.01. The homogeneous distribution 𝐻(𝑥,𝑦,0)=𝐻_∗, 𝑃(𝑥,𝑦,0)=𝑃_∗ is used as the initial condition and a small random perturbation is added to them. As the random perturbation propagates, system (2) will be evolved into the steady patterns, which are stationary in time and oscillatory in space. It should be pointed out that in our numerical simulations the spatial distributions of prey and predator are always the same type with change of 𝑚. It is because for system (2) the steady state predator is equal to its carrying capacity and this carrying capacity is assumed to be proportional to the number of prey. It implies that there is no difference between the effect of prey refuge on the patterns of prey and predator. Hence, spatial patterns of prey are only explored with change of 𝑚.

It should be stressed that the numerical simulations are run until the stationary states emerge, and in this way some snapshots shown in Figure 4 have been taken of numerical simulations when 𝑚 increases from 0.03 to 0.25. It is interesting to note from Figure 4 that the range of the changing densities of prey is given by the enclosed color bar in each snapshot, in which higher values denote higher prey densities. According to the dispersion relationship of system (2) shown in Figure 2, it is clear to know that when 𝑚≤0.01 or 𝑚≥0.27, system (2) gives rise to the stable state around the equilibrium point. Thus, the spatial distributions of prey under the condition of Turing instability are considered. For 𝑚=0.03, 𝑚=0.12, 𝑚=0.2, and 𝑚=0.25, the stationary snapshots are picked up to show the spatiotemporal dynamics of system (2) around the interior equilibrium point. Figure 4 reflects oscillations in space and the stationary state in time. Comparing the four diagrams shows that system (2) owns very rich and complex patterns which correspond to different habitats in real world. When 𝑚=0.03, the distribution of prey is mainly the irregularly elliptic patches, which exhibit a higher prey density in the centers of these patches. Moreover, the distribution of these patches is almost uniform in space. When 𝑚=0.12, the spatial pattern of prey is some stripes and a few blocks. It can be observed that on long strips there exist black dots which result from the congregation of population. Assuming that the value of 𝑚 continues to be increased, it shows that when 𝑚=0.2, the snapshot is obviously different from Figures 4(a) and 4(b). In the snapshot, there exist some holes on the long strips in which the prey density is very low. Furthermore, it can be observed that when 𝑚=0.25, the kind of holes prevails uniformly the whole domain. From a biological viewpoint, this type of irregularly elliptic structures shown in Figure 4(a) leads to a higher efficiency on reducing the pressure of predation by means of group defence than these structures shown in Figures 4(b)–4(d). From Figures 4(a)–4(d), the different distributions of one species denote the optimal section of self-organization in different habitats which is corresponding to the change of prey refuge. If the effect of prey refuge is relatively small, the prey will be evolved into the spatial pattern which can best reduce the attack rate of predator. Whereas the effect of prey refuge is so much strong that the prey can easily refuge from predation, then the population starts diffusion from the crowded patches to the low and tends to be uniform in spatial patterns. On the other hand, it is interesting to note from Figures 4(a)–4(d) that the maximum concentration of color in these snapshots shows the decreasing states as increase of 𝑚. This changing process also provides some information to explore the pattern formations. When the maximum density of prey is relatively high, the distribution of prey tends to be the irregularly elliptic patches for defending from predation. But the maximum density of prey is relatively low and the distribution of prey exhibits some holes in order to reduce the pressure among the prey individuals. In fact, these results show that the prey refuge not only determines the densities of prey and predator at steady state but also is beneficial to be evolved into adaption for their inhabits. More precisely, these results derive from the optimal evolution of species.

(a)

(b)

(c)

(d)

Figure 4 

Spatial distributions of prey as change of 𝑚. (a) 𝑚=0.03, (b) 𝑚=0.12, (c) 𝑚=0.2, and (d) 𝑚=0.25. The other parameter values are 𝑟=0.6, 𝑠=0.49, 𝑎=0.28, 𝑏=0.3, 𝑐=0.57, 𝑑=0.7, 𝑒=0.2, 𝑑1=0.08, and 𝑑2=0.12.

On the other hand, it is clear to be known from the definition of a prey refuge that the prey refuge factor incorporated in system (2) only denotes the individual behavior of species. However, according to the assumptions of system (2), the group defense stands for the cooperative behavior among the prey individuals. Therefore, it is necessary to investigate the effect of group defense on the distributions of prey. As a result, it can be conductive to understand that the population optimizes the strength of group defense when the population faces the pressures of both predation and intraspecific competition. Here, the strength of group defense is exhibited by means of the predation rate. If the predation rate is small, it indicates that the strength of group defense increases. Conversely, if the predation rate becomes large, the group defense is weak. Thus, the predation rate 𝑐 is chosen to report the effect of group defense on distributions of prey.

For further analysis of the effect of group defense on the dynamical behavior of system (2), it is first to consider the effect of group defense on the densities of prey and predator at steady state. Here, the density of prey is only plotted as a function of 𝑐 since from the expression of 𝑃_∗ it is clear to be seen that 𝑃_∗ as the function of 𝑐 has the same properties. According to Figure 5(a), it shows that as 𝑐 increases, the prey density at steady state will be decreased. It is also interesting to note from Figure 5(a) that there exist two inflection points. Assume that 𝑐₁, 𝑐₂ correspond to the horizontal values of two inflection points, respectively. More importantly, from Figure 5(b) the range of Turing instability lies in the interval (𝑐₁,𝑐₂). Furthermore, when 𝑐 is beyond this interval, then the states of system around the steady state become locally stable, and when these cases occur, the spatial distributions of prey will be chaotic [20]. Thus, the spatial distributions of prey under the condition of Turing instability are studied in detail. In order to improve the understanding of group defense, the spatiotemporal evolutions of prey corresponding to several values of 𝑐 are plotted and the snapshots are taken as the previous case. When 𝑐=0.51, the distribution of prey is the same type as Figure 4(d). When 𝑐=0.54, some irregular connected strips prevail and some small holes are incorporated into the large strips. When 𝑐=0.59 and 𝑐=0.61, the irregularly elliptic patches are formed. Comparing with Figure 4, the spatial distributions of prey as shown Figure 6 tend to be the elliptic patches and such distribution is helpful to reduce the pressure of predation as increase of 𝑐. In fact, when the predation rate becomes large, the prey tends to evolve into Figure 6(d) since this mode owns the higher strength of group defense and as a result it reduces the probability of extinction of species. Finally, from our numerical simulations, it shows that system (2) is more sensitive to the prey group defense than the prey refuge.

(a)

(b)

Figure 5 

(a) The predator density as a function of 𝑐 at steady state. (b) Variation of dispersion relation of system (2) around the interior equilibrium point for different 𝑐. The red line corresponds to 𝑐=0.5, the green line corresponds to 𝑐=0.57, and the blue corresponds to 𝑐=0.62. The other parameter values are 𝑟=0.6, 𝑠=0.49, 𝑎=0.28, 𝑏=0.3, 𝑑=0.7, 𝑒=0.2, 𝑚=0.12, 𝑑1=0.08, and 𝑑2=0.12.

(a)

(b)

(c)

(d)

Figure 6 

Spatial distributions of prey as change of 𝑐. (a) 𝑐=0.51, (b) 𝑐=0.54, (c) 𝑐=0.59, and (d) 𝑐=0.61. The other parameter values are 𝑟=0.6, 𝑠=0.49, 𝑎=0.28, 𝑏=0.3, 𝑑=0.7, 𝑒=0.2, 𝑚=0.12, 𝑑1=0.08, and 𝑑2=0.12.

From the previous analysis, it is obvious that the prey refuge and group defense have a significant influence on pattern formations of prey. Moreover, it is interesting to note that the distributions of prey are not only determined by the individual behavior but also result from the cooperation of prey. These results exhibit that species tend to be the optimal evolution to adapt to the environment and it is meaningful to understand the dynamical complexity of ecosystems in the real world.

4. Conclusions

In summary, a modified Leslie-Gower predator-prey with a prey refuge and group defense was investigated analytically and numerically. Initially, we performed the stability of system (5) at the interior equilibrium point and further obtained the range of the wavenumber such that system (2) leaded to the Turing bifurcation. Under the condition of Turing instability, the effects of prey refuge and group defense on the distributions of species were discussed in detail. Finally, the results of all simulations showed that the prey refuge and group defense had a highly significant effect on the spatiotemporal dynamics of system (2). It is expected to improve our understanding of the spatiotemporal dynamics of ecosystems in real world.

Acknowledgments

The authors would like to thank the editor and the anonymous referees for their valuable comments and suggestions on this paper. This work was supported by the National Key Basic Research Program of China (973 Program, Grant no. 2012CB426510), by the National Natural Science Foundation of China (Grant nos. 31170338 and 11226256), by the Key Program of Zhejiang Provincial Natural Science Foundation of China (Grant no. LZ12C0300), and by the Zhejiang Provincial Natural Science Foundation of China (Grant no. LY010010).