Effects of Behavioral Tactics of Predators on Dynamics of a Predator-Prey System

Abstract

A predator-prey model incorporating individual behavior is presented, where the predator-prey interaction is described by a classical Lotka-Volterra model with self-limiting prey; predators can use the behavioral tactics of rock-paper-scissors to dispute a prey when they meet. The predator behavioral change is described by replicator equations, a game dynamic model at the fast time scale, whereas predator-prey interactions are assumed acting at a relatively slow time scale. Aggregation approach is applied to combine the two time scales into a single one. The analytical results show that predators have an equal probability to adopt three strategies at the stable state of the predator-prey interaction system. The diversification tactics taking by predator population benefits the survival of the predator population itself, more importantly, it also maintains the stability of the predator-prey system. Explicitly, immediate contest behavior of predators can promote density of the predator population and keep the preys at a lower density. However, a large cost of fighting will cause not only the density of predators to be lower but also preys to be higher, which may even lead to extinction of the predator populations.

1. Introduction

The investigation of phenomena involving complex geometry, patterns, and scaling has undergone a substantial development in the past decades. In natural reality, many objects are “self-similar” at different scales, which is the essence of fractal. We see fractal patterns almost everywhere in nature, even the more random-appearing patterns of clouds or coastlines. Fractal has been applied into a variety of areas as a new concept and method, scaling is another important concept with it. In ecosystem, an ecological process including some small parts which happen on different scales is commonly observed; it is possible to take advantage of these time scales in order to reduce the dimension of the initial complete model and to build a simplified system that describes the dynamics of a small number of global variables. Aggregation methods are extended to systems with different time scales to reduce the number of variables [1–6].

Game theory as a very useful branch of mathematics has long been widely used to investigate the interaction between rational players in a game, the core problem in game theory is what kind of strategies will be adopted to improve player’s payoff as much as possible. Smith and Price [7] ingeniously linked its concept of payoff functions and players to the biological terminology of fitness and individuals, respectively, in order to explain the evolution of ritualised conflicts within species, for example, when individuals compete for mates or territory. Their seminal work thus inspired an entirely new approach to behavior ecology; evolutionary game theory is born to investigate problems involving frequency-dependent selection. It deals with entire populations of individuals, all programmed to use some strategies (or types of behavior); strategies with higher payoff will spread within population, through learning, copying, or inheriting, even by infection. The payoffs depend on the actions of coplayers and hence on the frequencies of the strategies within the populations.

In population biology, understanding the effects of interacting species at the population and community level is one of the most important issues. The classical predator-prey system is considered as a basic aspect and has long been investigated since Lotka in 1920 [8] followed by Volterra [9] from a theoretical point of view, respectively, where a simple coupled set of differential equations is devised to describe an autocatalytic reaction and the statistics of fish catches in the Adriatic. This first model (to be called Lotka-Volterra model (LVM)) has since become one of the central paradigms to elucidate the main features of species interactions and population cycles in predator-prey systems, where predator depends on the prey for food, and each species has an inherent birth and death rate. The two populations end up oscillating, with the predator population crests slightly offset from that of prey. Based on LVM, a proliferation of models has been employed to qualitatively capture and qualitatively understand the primary features of predator-prey systems as well as competition system (e.g., [10–12]), by considering factors which influence the stable coexistence of ecosystem [13] and maintenance of biodiversity [14–17].

Individual’s behaviors play an important role in population dynamics [18, 19]. Individuals have the ability of learning and can change tactics along their life according to the environment conditions and the results of past contests. Behavioral plasticity [20] allows an individual to be more flexible to adopt the behavior that can optimize its survival in the present environment [21]. Typically, aggressive behavior is commonly observed in animal kingdom when the prey resource is rare, thus making predators food-deprived to monopolize the limited resources, for example, in vertebrates [22], in birds [23], and in mammals [24]. Particularly, in a real world, predators are expected to choose their strategy both in relation to the density of prey and in relation with the strategies adopted by the other competitors. Inspired by this point, Auger et al. [6] studied a predator-prey system incorporating predators using the classical hawk and dove tactics. In their work, individual behavior structure of predators is firstly incorporated into the classical predator-prey model.

In this paper, a traditionally defined Lotka-Volterra interaction system is used to investigate the dynamics of predator-prey system, where predators use different strategies when disputing a prey. The canonical rock-paper-scissors (R-P-S) strategy [17] is assumed, with the second out-competing the first, the third out-competing the second, and the first out-competing the third, which has been invoked to capture the remarkable coexistence of three types of cyclically competition. The principle aim of this study is to explore the effect of individual behavior on the dynamics of the predator-prey system. In particular, we will look at the existence of coexistence equilibria between the prey and predator and the corresponding internal structure of the predator population: monomorphic or polymorphic state. We are interested in the coevolution at the individual and community level.

At first, we present the game dynamic model representing the fast scale part of the complete system in Section 2.1. Then the slow scale part is built up based on the classical LVM with prey-density dependence. Combining the fast and slow scale parts together, Sections 2.3 and 2.4 are devoted to the construction and analysis of the ensuing aggregated model. The dynamics of the system can be studied by mathematical analyses and computer simulations, some details of the analysis are given in the appendix.

2. Model and Analysis

Suppose that all predators hunt the same prey, after a predator catches a prey, this predator has to fight with others to keep the prey. Simplifying, we assume that once a prey is killed, a second predator in the vicinity arrives and is disputing the prey to the predator that has caught the prey. Predators can use three behavior tactics rock, paper, and scissors to dispute a prey when they meet. After a certain number of encounters, on average, the tactic that gives predators the better payoff will be carried out. This dynamic process can be described by replicator equations [25–30].

Because predators encounter and fight frequently (for example, each day) to keep a captured prey or to dispute a prey to another predator, we suppose the game dynamics correspond to a fast time scale. However, a few preys are killed each day in comparison to the total reservoir of prey, so the predator-prey interaction terms correspond to a slow time scale, involving birth, death, and predation.

Denote 𝑛(𝑡) and 𝑝(𝑡) as the densities of prey and predator population at time 𝑡, respectively. Next, we will construct models to describe the interaction of predators and prey.

2.1. Replicator Dynamics of Predators on a Fast Time Scale

Let us denote 𝑎_𝑖𝑗 (𝑖,𝑗∈{1,2,3}, 1, 2, 3 represent rock, paper, scissors, respectively) as the element of payoff matrix 𝐴, which corresponds to the gain that is obtained by an individual playing tactic 𝑖 against an individual playing tactic 𝑗. When two individuals with the same behavior meet, they share the gain 𝐺 equally but also the cost 𝐶 (𝐶>0) due to fighting, 𝑎_𝑖𝑖=(𝐺−𝐶)/2, and vice versa, the winner gets 𝐺 while the loser gets nothing. For example, when a rock strategy individual meets a paper strategy one, the first one gets nothing, but another gets 𝐺. Consequently, the payoff matrix 𝐴 takes the following form: 𝐀=𝐺−𝐶2𝐺0𝐺𝐺−𝐶200𝐺𝐺−𝐶2. (1) The gain 𝐺 of the game corresponds to the prey amount that two predators dispute each day. To simplify, a classical type I linear functional response is taken into account, that is, 𝐺=𝑎𝑛, (2) where 𝑎 is a positive parameter.

Denote 𝑝_𝑟(𝑡), 𝑝_𝑝(𝑡), and 𝑝_𝑠(𝑡) as the densities of individuals with rock, paper, and scissors strategies at time 𝑡, respectively. The total density of predators will be 𝑝(𝑡)=𝑝_𝑟(𝑡)+𝑝_𝑝(𝑡)+𝑝_𝑠(𝑡). (3) Furthermore, let 𝑥₁(𝑡), 𝑥₂(𝑡), and 𝑥₃(𝑡) be respectively the rock, paper, and scissor proportions in the population of predators at time 𝑡 as follows: 𝑥₁(𝑡)=𝑝_𝑟(𝑡)𝑝(𝑡),𝑥₂(𝑡)=𝑝_𝑝(𝑡)𝑝(𝑡),𝑥₃(𝑡)=1−𝑥₁(𝑡)−𝑥₂(𝑡). (4) The replicator equations describing the evolutionary dynamics of the three strategies read as follows: d𝑥_𝑖d𝜏=𝑥_𝑖(𝐀𝐱)_𝑖−𝐱^𝑇𝐱𝐀𝐱,𝑖=1,2,3,_𝑇𝐀𝐱=³_𝑗=1𝑥_𝑗(𝐀𝐱)_𝑗, (5) where 𝐱=(𝑥₁,𝑥₂,𝑥₃). Obviously, if a strategy brings a better payoff than the average payoff of the population, the proportion of individuals playing this strategy is increasing and conversely is decreasing. In (5), 𝜏 is the fast time scale. Thus, we assume that the R-P-S game is in fast time scale compared to the processes that we will now consider in the model, such as predator death, prey growth, and capturing of preys by predators.

2.2. Dynamics of Predator-Prey at the Slow Time Scale

A classical Lotka-Volterra model with self-limiting prey [31] is assumed to describe the total prey density as follows: d𝑛𝑛d𝑡=𝑟𝑛1−𝐾−𝑎𝑛𝑝, (6) where 𝑡 corresponds to the slow time scale. Usually, we have the relationship 𝑡=𝜖𝜏 for two time scales. From (6) we know, in absence of predators, that the growth of prey population depends on a logistic equation with an intrinsic growth rate 𝑟 towards a carrying capacity 𝐾. Besides, we assume a Lotka-Volterra functional response of type 𝐼; the parameter 𝑎 is positive, which represents predation force and is the same as the one in the functional response.

Predators are assumed to have a constant natural mortality rate 𝜇 (𝜇>0) no matter what kind of strategies they use., the growth rate of each subgroup with the same strategy is proportional to the average payoff obtained by an individual using the tactic on the occasion of each type of encounter. For example, a rock player can encounter either a rock player in proportion of 𝑝_𝑟/𝑝 and gets the gain (𝐺−𝐶)/2, or a paper player in proportion of 𝑝_𝑝/𝑝 and gets zero, or as a proportion of 𝑝_𝑠/𝑝 to encounter a scissor player and gets 𝐺. Consequently, the growth of the rock predator subpopulation is described as follows: d𝑝_𝑟d𝑡=−𝜇𝑝_𝑟+𝑝_𝑟𝛼𝐺−𝐶2𝑝_𝑟𝑝𝑝+𝛼𝐺^𝑠𝑝, (7) where 𝛼 is a conversion positive coefficient of gain and cost into biomass of predators. Analogously, it also holds for the paper subpopulation of predators as follows: d𝑝_𝑝d𝑡=−𝜇𝑝_𝑝+𝑝_𝑝𝑝𝛼𝐺^𝑟𝑝+𝛼𝐺−𝐶2𝑝_𝑝𝑝. (8) For the scissors subpopulation of predators, we obtain d𝑝_𝑠d𝑡=−𝜇𝑝_𝑠+𝑝_𝑠𝑝𝛼𝐺^𝑝𝑝+𝛼𝐺−𝐶2𝑝_𝑠𝑝. (9)

2.3. The Complete Slow-Fast Predator-Prey Model

By combing both fast and slow processes in a similar way as in previous papers [6, 32–35], the complete model reads as follows: 𝜖d𝑛𝑛d𝑡=𝜖𝑟𝑛1−𝐾.𝜖−𝑎𝑛𝑝d𝑝_𝑟d𝑡=𝜖−𝜇𝑝_𝑟+𝑝_𝑟𝛼𝐺−𝐶2𝑝_𝑟𝑝𝑝+𝛼𝐺^𝑠𝑝+𝑝𝑥₁(𝐀𝐱)₁−𝐱^𝑇.𝜖𝐀𝐱d𝑝_𝑝d𝑡=𝜖−𝜇𝑝_𝑝+𝑝_𝑝𝑝𝛼𝐺^𝑟𝑝+𝛼𝐺−𝐶2𝑝_𝑝𝑝+𝑝𝑥₂(𝐀𝐱)₂−𝐱^𝑇.𝜖𝐀𝐱d𝑝_𝑠d𝑡=𝜖−𝜇𝑝_𝑠+𝑝_𝑠𝑝𝛼𝐺^𝑝𝑝+𝛼𝐺−𝐶2𝑝_𝑠𝑝+𝑝𝑥₃(𝐀𝐱)₃−𝐱^𝑇𝐀𝐱, (10) where 𝜖≪1 is a small parameter. By using the fast time scale 𝜏, the system changes to d𝑛𝑛d𝜏=𝜖𝑟𝑛1−𝐾−𝑎𝑛𝑝.d𝑝_𝑟d𝜏=𝜖−𝜇𝑝_𝑟+𝑝_𝑟𝛼𝐺−𝐶2𝑝_𝑟𝑝𝑝+𝛼𝐺^𝑠𝑝+𝑝𝑥₁(𝐀𝐱)₁−𝐱^𝑇𝐀𝐱.d𝑝_𝑝d𝜏=𝜖−𝜇𝑝_𝑝+𝑝_𝑝𝑝𝛼𝐺^𝑟𝑝+𝛼𝐺−𝐶2𝑝_𝑝𝑝+𝑝𝑥₂(𝐀𝐱)₂−𝐱^𝑇𝐀𝐱.d𝑝_𝑠d𝜏=𝜖−𝜇𝑝_𝑠+𝑝_𝑠𝑝𝛼𝐺^𝑝𝑝+𝛼𝐺−𝐶2𝑝_𝑠𝑝+𝑝𝑥₃(𝐀𝐱)₃−𝐱^𝑇𝐀𝐱. (11) Obviously, the game dynamics correspond to the fast time scale while the small terms of the order of 𝜖 correspond to the slow time scale. Next, aggregation method will be carried to reduce the dimension of the system.

2.4. The Aggregated Predator-Prey Model

2.4.1. Derivation of the Aggregated Model

Since 𝜖≪1, we neglect the small terms of the order of 𝜖 and do stability analysis for the fast part of the system which relates to the game dynamics as follows: d𝑝_𝑟d𝜏=𝑝𝑥₁(𝐀𝐱)₁−𝐱^𝑇𝐀𝐱,d𝑝_𝑝d𝜏=𝑝𝑥₂(𝐀𝐱)₂−𝐱^𝑇𝐀𝐱,d𝑝_𝑠d𝜏=𝑝𝑥₃(𝐀𝐱)₃−𝐱^𝑇𝐀𝐱. (12) Since we are only interested in positive solutions, by realizing the fact that 𝑥₁+𝑥₂+𝑥₃=1 always holds at any time 𝑡, we can limit our analysis on the region called a three-dimensional simplex, 𝑆³={(𝑥₁,𝑥₂,𝑥₃)∈𝑅³₊∶𝑥₁+𝑥₂+𝑥₃=1}. By letting the right-hand side of (12) equal zero, the equilibria of the system can be derived, which depends on the parameters 𝐺 and 𝐶.

The system has an interior stable focus (1/3,1/3,1/3) no matter what the parameters are. There exist three fixed points (1,0,0), (0,1,0), and (0,0,1) located on the vertex of 𝑆³; they are all saddle points given 𝐺>𝐶, unstable equilibrium points for the inverse condition of 𝐺<𝐶. Besides, three boundary equilibria ((𝐺+𝐶)/2𝐶,0,(−𝐺+𝐶)/2𝐶)), ((−𝐺+𝐶)/2𝐶),(𝐺+𝐶)/2𝐶,0), and (0,(−𝐺+𝐶)/2𝐶),(𝐺+𝐶)/2𝐶) can also exist as saddle points provided that 𝐺<𝐶. Therefore, 𝑥^∗=(1/3,1/3,1/3) is the only nontrivial interior asymptotically stable point without depending on initial conditions. In this case, the predator population is a polymorphic state with equal probability of three tactics of predators.

In order to aggregate, we make the assumption that the fast process is at fast equilibrium. Then we come back to the complete initial system (10), substitute the previous fast equilibrium and add the three predator equations together. It is necessary to replace the fast variables in terms of the fast equilibrium which are 𝑝_𝑟=𝑝/3, 𝑝_𝑝=𝑝/3, and 𝑝_𝑠=𝑝/3.

After some algebra, one obtains the following system of two equations governing the total prey and predator densities at the slow time scale, that we call the aggregated model (AM) as follows: d𝑛𝑛d𝑡=𝑟𝑛1−𝐾−𝑎𝑛𝑝+𝑜(𝜖).d𝑝d𝑡=−𝜇𝑝+𝛼3𝐺−𝐶2𝑝+𝑜(𝜖). (13) Because 𝜖 is small enough and the system is structurally stable [36]; it is a good approximation of the real dynamics. Then by neglecting these terms of higher order of 𝜖 and realizing the fact that 𝐺(𝑛)=𝑎𝑛, the aggregated model simply reads as follows: d𝑛=𝑛d𝑡𝑟−𝑟𝐾−𝑎𝑝𝑛d𝑝=d𝑡−𝜇+3𝛼𝑎𝑛2−𝛼𝐶2𝑝. (14)

2.4.2. Analyses of the Dynamics of the Aggregated Model

The phase plane {𝑛,𝑝} is separated by the vertical nullcline 𝐿 denoted as 𝑟𝐾−𝑟𝑛−𝑎𝑝𝐾=0 and the horizontal nullcline 𝑀 denoted as −2𝜇+3𝛼𝑎𝑛−𝛼𝐶=0 into different regions, both the signs of d𝑛/d𝑡 and d𝑝/d𝑡 do not change in every part. Next, the dynamics of the system will be studied through the analysis of the intersection of the lines 𝐿 and 𝑀. (1) 𝐾<(2𝑢+𝛼𝐶)/(3𝛼𝑎), 𝐿 and 𝑀 are not intersected in phase plane.

In this case, the system has two equilibria, where (0,0) is a saddle point while (𝐾,0) is an asymptotically stable point (see Appendix A.1). There does not exist any limit cycle. We know that closed trajectory should enclose equilibrium point, but there is no equilibrium point in the interior part of the phase plane. Therefore, the trajectories that originated from the vertical axis representing 𝑛(𝑡)=0 will approach to the point (0,0); otherwise, all trajectories will end at the point (𝐾,0). That is to say, the dynamic result will die out of predators no matter what the initial densities of prey and predators are. Once the initial density of prey is not zero, it will reach to a stable density with a value of 𝐾 in the end. (2) 𝐾>(2𝑢+𝛼𝐶)/(3𝛼𝑎), 𝐿 and 𝑀 intersect in the phase plane.

The system includes three equilibria, (0,0) and (𝐾,0) are saddle points, the interior point 𝑥=(𝑛,𝑝)=((2𝑢+𝛼𝐶)/(3𝛼𝑎),𝑟(3𝛼𝑎𝐾−2𝑢−𝛼𝐶)/(3𝛼𝑎²𝐾)) standing for the intersection point of 𝐿 and 𝑀 is an asymptotically stable point (see Appendix A.1). In addition, we can testify that there is no limit cycle, the details are written in Appendix A.3. Therefore, once the initial densities of prey and predators are all larger than zero, their density will finally reach the point 𝑥.

In reality, 𝜎=(2𝑢+𝛼𝐶)/(3𝛼𝑎) is the least food requirement for viability of predators, and 𝐾 is the carrying capacity. 𝜎>𝐾 means that the least food requirement for predators is larger than the reproduction and fecundity of prey, so predators will extinct without any doubt, or else, the predators and preys will coexist together.

3. Results

A classical LVM will be studied in this part in order to give a comparison with our AM. Suppose that there is no difference between the strategy adopted by predators, the predator-prey interaction can be described as follows: d𝑛𝑛d𝑡=𝑟𝑛1−𝐾−𝑎𝑛𝑝d𝑝d𝑡=−𝜇𝑝+𝛼𝑎𝑛𝑝. (15) By using a similar method as in the last model, the dynamics can also be recognized to be two cases. When 𝐾<𝑢/(𝛼𝑎), the system has two equilibria, (0,0) is a saddle point, (𝐾,0) is an asymptotically stable point. While for 𝐾>𝑢/(𝛼𝑎), the system has three equilibria, both (0,0) and (𝐾,0) are saddle points, 𝑥=(𝑛,𝑝)=(𝑢/(𝛼𝑎),𝑟(𝛼𝑎𝐾−𝑢)/(𝛼𝑎²𝐾)) is an asymptotically stable interior point. The details are given in Appendix A.2. In addition, we can verify there is no limit cycle in the phase plane by using the same method as in Appendix A.3. In summary, all results are concluded in Table 1.

Table 1 

The evolutionary stable state in LVM and AM.

Therefore, we find that the stable state of the predator-prey system will be changed somehow considering that predator population uses different strategies when disputing a prey, which can be recognized as four cases depending on different parameter settings. The stable state can change from only preys with a value at carrying capacity in LVM to coexistence between preys and predators in AM, Figures 1(a) and 1(b) give an example. It also shows from the coexistence of prey and predators in LVM to only prey with a value at carrying capacity in AM (see Figures 2(a) and 2(b)). In addition, the case that the prey population with the density of carrying capacity but predators dying out can be held. Coexistence of the predator and prey populations but with a variation of population densities in the two models can also be kept. When 𝐶<𝑢/𝛼, for the parameter region (2𝑢+𝛼𝐶)/(3𝛼𝑎)<𝐾<𝑢/(𝛼𝑎), the density of prey in the LVM will be higher than that in the AM at the stable state, but the result is inverse for the predators, as shown in Figures 1(c) and 1(d); while for 𝐶>𝑢/𝛼, when the parameters satisfy 𝑢/(𝛼𝑎)<𝐾<(2𝑢+𝛼𝐶)/(3𝛼𝑎), the density of prey in the LVM will be lower than that in the AM, but it is inverse for the predator populations, Figures 2(c) and 2(d) are an example. In conclusion, we can note that fighting between predators could promote the growth of the predator population when the coming cost is smaller than the division of natural mortality of prey by conversion coefficient, so a moderate degree of fighting will also benefit the evolution of the predator population. In other words, a slight extent of competition is helpful for the survival of the population. On the other side, fighting will naturally cause a negative effect on the survival of predator population when the cost is relatively large, thus making an acceleration of the prey population, this phenomenon confirms the reality.

(a)

(b)

(c)

(d)

Figure 1 

The dynamics of LVM and AM given the condition 𝐶<𝑢/𝛼. Parameters are 𝑟=5, 𝐾=23/12, 𝑎=1, 𝛼=1, 𝐶=3/2, and 𝑢=2 in the first two figures, the carrying capacity is set to be 𝐾=3 in the other two.

(a)

(b)

(c)

(d)

Figure 2 

The dynamics of LVM and AM given the condition 𝐶>𝑢/𝛼. Parameters are 𝑟=5, 𝐾=13/12, 𝑎=1, 𝛼=1, 𝐶=3/2, and 𝑢=2 in the first two figures, all parameter values are the same in the other two except 𝐾=2.

4. Conclusions and Discussion

A vulnerable subject in population ecology is understanding that how diversity and complexity of natural populations contribute to the overall stability and persistence of ecological systems. The idea that diversity promotes stability is intuitively appealing. In this work, the stability of a predator-prey interaction system with predators using diversification strategies is investigated to study the effect of different tactics taken by predators when disputing a prey, thus providing us a better understanding of the important effect of biodiversity on the stability of the ecosystem.

In this study, a dynamic modeling method [37] is employed to explore the dynamics of a predator-prey system by connecting individual behavior with population dynamics. Particularly, by resorting to game theory, replicator equations are constructed to derive ESS of the fast dynamics describing the variation of frequencies of different tactics. Here, by taking the R-P-S strategy to investigate the effect of predator behavior on the stability of predator-prey system, we found that the ESS in fast time scale of disputing a prey is independent of parameters and is very symmetric, which locates on the center of the strategy space. This result is dramatically different from that by Auger et al. [6], where the ESS in the fast time scale will stay on different equilibrium points with a variation of the density of preys in their system. Therefore, the main conclusion of their studies is that there exists a relationship between prey density and the strategy adopted by predators. Aggressive behavior is connected to high prey and low predator densities, whereas a polymorphism dove-hawk is found at low prey and high predator densities. This result is consistent well with the previous studies of a domestic cat population [33, 35]. In rural areas, low density, most of the individuals are aggressive while in urban areas, high density, most of them are dove individuals.

In addition, we found that the equilibrium state of the system was determined by parameters 𝐾, 𝑢, 𝛼, 𝑎, 𝐶, 𝑟. More importantly, the change of the stable state by considering predator behavior with R-P-S strategy through comparing with the common Lotka-Volterra system without considering it can be observed, which also depends on 𝐾, 𝑢, 𝛼, 𝑎, 𝐶, but not relates to 𝑟 (the birth rate of prey). Explicitly, it can be firstly recognized as two classes according to the numerical relationship of 𝐶 and 𝑢/𝑎. Furthermore, in every class, constrained by the relationship of 𝐾 with 𝑢/𝛼𝑎 and (2𝑢+𝛼𝐶)/(3𝛼𝑎), it appears that either a change from the extinction of predator populations in LVM to the coexistence state with predators and preys in AM, or from the coexistence in LVM to the extinction of predator population in AM, or coexistence state will be kept in both situations but with a change of the density of the population.

Especially, it is worth noting that given the conditions of 𝐶<𝑢/𝑎 and 𝐾>𝑢/(𝛼𝑎), that is, when carrying capacity of predators is large enough and fighting cost of predators is small, coexistence state will be kept in both Lotka-Volterra and Aggregated system, it shows an increase of the density of predator population and a decrease of the density of prey population. This verified the nature rule that an immediate competition will benefit the survival of the population. We also concluded that a diversification behavior of the predator population would maintain the stability of the ecosystem for certain parameters; that is, 𝐶<𝑢/𝑎 and (2𝑢+𝛼𝐶)/(3𝛼𝑎)<𝐾<𝑢/(𝛼𝑎). Predator population will not go extinct by taking different strategies when disputing a prey.

In future, some more interesting researches are expected to expand by considering other intricate strategies taken by predators which would raise diverse ESS situations depending on different levels of prey densities, thus leading to diversification of the system. In addition, motivated by Huffaker’s experiments [38], it has been documented that the distribution of prey and predators in the spatial configuration as well as the dispersal behaviors are all important factors to influence the dynamics of the predator-prey system; these factors play an important role in understanding the emergence of biodiversity and the stability of ecosystems (i.e., [39–42]). From this point of view, combining spatial effect into our model also deserves research, further studies are expected along this direction.

Appendix

A. Local Stability of the Equilibria 

A.1.  Case of Model (14)

The stability analysis of the aggregated model (14), the Jacobian matrix is 𝑟𝐉=1−2𝑛𝐾−𝑎𝑝−𝑎𝑛3𝛼𝑎𝑝2−𝜇+3𝛼𝑎𝑛2−𝛼𝐶2. (A.1) At the origin point (0,0), we have 𝐉|_(𝟎,𝟎)=𝑟00−𝜇−𝛼𝐶2. (A.2) It has two real eigenvalues with opposite signs, the origin point is a saddle point.

For the equilibrium point (𝐾,0), the Jacobian matrix reads as follows: 𝐉|_(𝐊,𝟎)=−𝑟−𝑎𝐾0−𝜇+3𝛼𝑎𝐾2−𝛼𝐶2, (A.3) if 𝐾<(2𝑢+𝛼𝐶)/3𝛼𝑎, (𝐾,0) is an asymptotically stable point, while for 𝐾>(2𝑢+𝛼𝐶)/3𝛼𝑎, a saddle point.

For the positive equilibrium point (𝑛,𝑝)=((2𝑢+𝛼𝐶)/(3𝛼𝑎),𝑟(3𝛼𝑎𝐾−2𝑢−𝛼𝐶)/(3𝛼𝑎²𝐾)), the Jacobian matrix reads as follows: 𝐉|_(𝐧,𝐩)=−𝑟𝐾𝑛−𝑎𝑛3𝛼𝑎𝑝20. (A.4) So it is straightforward that when the initial state of the population belongs to the positive quadrant, the trace of the Jacobian matrix is negative and the determinant is positive which implies stability.

A.2. Case of Model (15) 

The stability analysis of the Lotka-Volterra model (15): the Jacobian matrix is 𝑟𝐉=1−2𝑛𝐾.−𝑎𝑝−𝑎𝑛𝛼𝑎𝑝−𝜇+𝛼𝑎𝑛 (A.5) At the origin point, we have 𝐉|_(𝟎,𝟎)=𝑟00−𝜇 (A.6) with two real eigenvalues with opposite signs. The origin point is a saddle point.

For the equilibrium point (𝐾, 0), the Jacobian matrix reads as follows: 𝐉|_(𝐊,𝟎)=.−𝑟−𝑎𝐾0−𝜇+𝛼𝑎𝐾 (A.7) if 𝐾>𝜇/(𝛼𝑎), (𝐾, 0) is a saddle point, while for 𝐾<𝜇/(𝛼𝑎), (𝐾, 0) is a stable node.

For the positive equilibrium point 𝑥=(𝑛,𝑝)=(𝑢/𝛼𝑎, (𝑟/𝑎)(1−𝑢/𝛼𝑎𝐾)), the Jacobian matrix reads as follows: 𝐉|_(𝐧,𝐩)=−𝑟𝐾𝑛,𝑛−𝑎𝛼𝑎𝑝0 (A.8) where it is obvious that when the point belongs to the positive quadrant, the trace of the Jacobian matrix is negative; consequently, the determinant is positive so the equilibrium point is stable.

A.3.  Limit Cycle Analysis to System (14)

We assume that there is a closed orbit Γ₀ around the interior point 𝑥, the equation is 𝑛=𝑛(𝑡),𝑝=𝑝(𝑡),−∞<𝑡<+∞. (A.9) Then there must be a 𝑇>0, where 𝑛(𝑇)=0, 𝑝(𝑇)=0. Next, we calculate the index 𝛾₀ of Γ₀, 𝛾₀=1𝑇^𝑇₀𝑃^′_𝑛𝑛(𝑡),𝑝(𝑡)+𝑄^′_𝑝𝑛(𝑡),𝑝(𝑡)𝑑𝑡, (A.10) where 𝑃(𝑛, 𝑝)=𝑛(𝑟−𝑟(𝑛/𝐾)−𝑎𝑝) and 𝑄(𝑛, 𝑝)=𝑝(−𝜇+3𝛼𝑎𝑛/2−𝛼𝐶/2). Thus, (15) can be rewritten as 𝑛𝑟−𝑟𝐾−𝑎𝑝=d𝑛𝑛−𝜇+3𝛼𝑎𝑛2−𝛼𝐶2=d𝑝𝑝. (A.11) Calculating the equations on two sides, we can get 1𝑇^𝑇₀𝑛(𝑡)𝑟−𝑟𝐾−𝑎𝑝(𝑡)𝑑𝑡=_𝑛(0)^𝑛(𝑇)d𝑛𝑛=ln|𝑛||_𝑛(0)^𝑛(𝑇)1𝑇^𝑇₀(𝑡)−𝜇+3𝛼𝑎𝑛2−𝛼𝐶2𝑑𝑡=_𝑝(0)^𝑝(𝑇)d𝑝𝑝=ln𝑝_𝑝(0)^𝑝(𝑇). (A.12) It shows that (1/𝑇)∫^𝑇₀𝑛(𝑡)𝑑𝑡 and (1/𝑇)∫^𝑇₀𝑝(𝑡)𝑑𝑡 also satisfy the following algebra equation: 𝑛𝑟−𝑟𝐾−𝑎𝑝=0−𝜇+3𝛼𝑎𝑛2−𝛼𝐶2=0. (A.13) Therefore, 1𝑇^𝑇₀𝑛(𝑡)𝑑𝑡=𝑛1𝑇^𝑇₀𝑝(𝑡)𝑑𝑡=𝑝. (A.14) Then the index of Γ₀ is 𝛾₀=1𝑇^𝑇₀𝑛𝑟−2𝑟𝐾−𝑎𝑝+−𝜇+3𝛼𝑎𝑛2−𝛼𝐶2=1𝑑𝑡𝑇^𝑇₀𝑛−𝑟𝐾𝑑𝑡=−𝑟𝑛𝐾<0. (A.15) Therefore, the closed orbit should be stable, it contradicts that 𝑥 is an asymptotically stable equilibrium, so the assumption is wrong; there are not any closed orbits.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Hui Zhang gratefully acknowledges the support from the National Natural Science Foundation of China (no. 31300310, no. 71171163, and no. 71271171), the Foundation Sciences Northwestern Polytechnical University (no. GDKY1007), the Projects of International Cooperation and Exchanges NSFC (no. 71311120091), and the Shaanxi Province Science and Technology Research and Development Programm (no. 2014KW03-01). Zhihui Ma gratefully acknowledges the support of the National Natural Science Foundation of China (no. 11301238).