- About this Journal
- Abstracting and Indexing
- Aims and Scope
- Annual Issues
- Article Processing Charges
- Articles in Press
- Author Guidelines
- Bibliographic Information
- Citations to this Journal
- Contact Information
- Editorial Board
- Editorial Workflow
- Free eTOC Alerts
- Publication Ethics
- Reviewers Acknowledgment
- Submit a Manuscript
- Subscription Information
- Table of Contents
Mathematical Problems in Engineering
Volume 2013 (2013), Article ID 910349, 8 pages
Complex Dynamical Behavior of a Predator-Prey System with Group Defense
1School of Mathematics and Information Science, Wenzhou University, Wenzhou, Zhejiang 325035, China
2School of Life and Environmental Science, Wenzhou University, Wenzhou, Zhejiang 325035, China
Received 10 March 2013; Revised 8 June 2013; Accepted 13 June 2013
Academic Editor: Hai-Feng Huo
Copyright © 2013 Jianglin Zhao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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.
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  is introduced, which can be expressed as follows: 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  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 is constant and 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: 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: and zero-flux boundary conditions:
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
From the biological angle, the steady state of system (2) or (5) is most significant, which implies coexistence of both predator and prey, and is the positive solution of , , where is the positive root of 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 . Assuming that , the function is continuous in , and . From the Zero Theorem, there must exist at least one positive real number such that if ; that is to say .
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: where and Then, the characteristic equation of system (7) is where and .
It is obvious that if and , system (7) has two distinct eigenvalues with negative real parts. Thus, is locally asymptotically stable for system (5). Assuming that , , , , , , , and , we can get that , and thus there exists . Further, we can obtain that and . Therefore, is locally asymptotically stable for system (5) shown in Figure 1.
In order to investigate the temporal stability of the uniform state , the nonuniform perturbations for any solution of system (2) are introduced as follows: 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 where and .
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 , . According to (13), we can obtain . 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 for , where . 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 , , , , , , , , and 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 or , the Turing instability is no longer excited.
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.
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 . 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 . The homogeneous distribution , 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 or , 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 , , , and , 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 , 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 , 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 , 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 , 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.
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 . 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 , the distribution of prey is the same type as Figure 4(d). When , some irregular connected strips prevail and some small holes are incorporated into the large strips. When and , 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.
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.
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.
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).
- J. S. Tener, Muskoen, Queen's Printer, Ottawa, Canada, 1965.
- J. C. Holmes and W. M. Bethel, “Modification of intermediate host behaviour by parasites,” Behavioural Aspects of Parasite Transmission, vol. 51, pp. 123–149, 1972.
- P. Davidowwicz, Z. M. Gliwicz, and R. D. Gulati, “Can Daphnia prevent a blue-green algal bloom in hypertrophic lakes? A laboratory test,” Limnologica, vol. 19, no. 1, pp. 21–26, 1988.
- H. I. Freedman and S. G. Ruan, “Hopf bifurcation in three-species food chain models with group defense,” Mathematical Biosciences, vol. 111, no. 1, pp. 73–87, 1992.
- Y. Liu, “Geometric criteria for the non-existence of cycles in predator-prey systems with group defense,” Mathematical Biosciences, vol. 208, no. 1, pp. 193–204, 2007.
- S. Y. Li, Z. L. Xiong, and X. Wang, “The study of a predator-prey system with group defense and impulsive control strategy,” Applied Mathematical Modelling, vol. 34, no. 9, pp. 2546–2561, 2010.
- M. Zhao and L. Zhang, “Permanence and chaos in a host-parasitoid model with prolonged diapause for the host,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 12, pp. 4197–4203, 2009.
- M. Zhao and S. J. Lv, “Chaos in a three-species food chain model with a Beddington-DeAngelis functional response,” Chaos, Solitons & Fractals, vol. 40, no. 5, pp. 2305–2316, 2009.
- M. Zhao, X. T. Wang, H. G. Yu, and J. Zhu, “Dynamics of an ecological model with impulsive control strategy and distributed time delay,” Mathematics and Computers in Simulation, vol. 82, no. 8, pp. 1432–1444, 2012.
- H. G. Yu and M. Zhao, “Seasonally perturbed prey-predator ecological system with the Beddington-DeAngelis functional response,” Discrete Dynamics in Nature and Society, vol. 2012, Article ID 150359, 12 pages, 2012.
- A. F. Nindjin, M. A. Aziz-Alaoui, and M. Cadivel, “Analysis of a predator-prey model with modified Leslie-Gower and Holling-type II schemes with time delay,” Nonlinear Analysis. Real World Applications. An International Multidisciplinary Journal, vol. 7, no. 5, pp. 1104–1118, 2006.
- M. X. Wang, P. Y. H. Pang, and W. Y. Chen, “Sharp spatial patterns of the diffusive Holling-Tanner prey-predator model in heterogeneous environment,” IMA Journal of Applied Mathematics, vol. 73, no. 5, pp. 815–835, 2008.
- S. L. Yuan and Y. L. Song, “Bifurcation and stability analysis for a delayed Leslie-Gower predator-prey system,” IMA Journal of Applied Mathematics, vol. 74, no. 4, pp. 574–603, 2009.
- X. N. Guan, W. M. Wang, and Y. Cai, “Spatiotemporal dynamics of a Leslie-Gower predator-prey model incorporating a prey refuge,” Nonlinear Analysis. Real World Applications. An International Multidisciplinary Journal, vol. 12, no. 4, pp. 2385–2395, 2011.
- W. Sokol and J. A. Howell, “Kinetics of phenol oxidation by washed cells,” Biotechnology and Bioengineering, vol. 23, no. 9, pp. 2039–2049, 1981.
- F. Chen, L. Chen, and X. Xie, “On a Leslie-Gower predator-prey model incorporating a prey refuge,” Nonlinear Analysis. Real World Applications. An International Multidisciplinary Journal, vol. 10, no. 5, pp. 2905–2908, 2009.
- Z. H. Ma, W. L. Li, Y. Zhao, W. T. Wang, H. Zhang, and Z. Z. Li, “Effects of prey refuges on a predator-prey model with a class of functional responses: the role of refuges,” Mathematical Biosciences, vol. 218, no. 2, pp. 73–79, 2009.
- J. Rebaza, “Dynamics of prey threshold harvesting and refuge,” Journal of Computational and Applied Mathematics, vol. 236, no. 7, pp. 1743–1752, 2012.
- T. K. Kar, “Stability analysis of a prey-predator model incorporating a prey refuge,” Communications in Nonlinear Science and Numerical Simulation, vol. 10, no. 6, pp. 681–691, 2005.
- G. Q. Sun, Z. Jin, Q. X. Liu, and L. Li, “Pattern formation in a spatial S-I model with non-linear incidence rates,” Journal of Statistical Mechanics, vol. 2007, no. 11, Article ID P11011, 2007.