DYNAMICS OF A PREDATOR-PREY MODEL WITH STATE-DEPENDENT CARRYING CAPACITY

Vegetation and plateau pika are two key species in alpine meadow ecosystems on the Tibetan Plateau. It is frequently observed on the field that plateau pika reduces the carrying capacity of vegetation and the mortality of plateau pika increases along with the increasing height of vegetation. This motivates us to propose and study a predator-prey model with state-dependent carrying capacity. Theoretical analysis and numerical simulations show that the model exhibits complex dynamics including the occurrence of saddle-node bifurcation, transcritical bifurcation and Hopf bifurcation, and the coexistence of two stable equilibria.

1. Introduction.The alpine meadow is an important ecosystem in China that provides important ecological services such as purification of air and water, maintenance of biodiversity, decomposition of wastes, soil and vegetation generation and renewal, etc.Over the past few decades, it has been observed that many alpine meadows have been degraded resulting in severe consequences in the ecology, environment and economy [19,25,26].For example, the degradation of alpine meadow can reduce species diversity and vegetation cover resulting in hydrological disturbances, dust storms, commodity scarcity.To effectively restore a degraded alpine meadow, it is vital to identify the underlying causes of alpine meadow degradation and evaluate the efficiency of possible restoration strategies.This requires a better understanding of the relationships among the interacting species in the alpine meadow ecosystem.
The small mammal plateau pika (Ochotona curzoniae) is a keystone species in the alpine meadow ecosystem.In a degraded alpine meadow, plateau pika is often in high density and is labeled as a pest subjected to control measures [6,15].There are more complex relationships between plateau pika and meadow vegetation beyond that plateau pika feeds on vegetation.A common fact is that when a plateau pika digs a hole, a lot of soil is thrown out, piling up a mound.The mound covers the existing vegetation and consequently there will be no vegetation on the mound for a long time due to the harsh climate.Hence the holes and mounds reduce the area where vegetation can grow and thus reduce the carrying capacity of vegetation [20].It has also been noticed that the height of vegetation has a negative impact on the population of plateau pika: higher vegetation prevents a plateau pika from noticing its predators such as owls and has less influence on its predators, so more plateau pika are caught or move to other areas to avoid being caught [14].In general, higher vegetation is accompanied with fewer plateau pika and shorter vegetation is accompanied with more plateau pika.
Though the relationship between plateau pika and meadow vegetation falls into the category of predator-prey relationship, the classical Rosenzweig-MacArthur model [3,12] fails to take the above mentioned facts into account.Motivated by the special features in the interaction between plateau pike and vegetation, we propose the following predator-prey model with state-dependent carrying capacity Here x(t) is the dry weight of vegetation at time t and y(t) denotes the amount/density of plateau pika at time t.The functional response of Holling-II type [4] is qx α+x .Parameter q is the maximal per-capita feeding rate and α is the half-saturation constant, p denotes the maximum per-capita growth rate of plateau pika.Parameter r 2 is the per-capita natural mortality rate of plateau pika and µ is the per-capita death rate related to the height of vegetation which is assumed to be proportional to x(t).Parameter r 1 is the intrinsic per-capita growth rate of vegetation, and K is the carrying capacity of vegetation in the absence of plateau pika.The parameter λ accounts for the reduction in carrying capacity due to holes and mounds produced by plateau pika.All parameters are positive constants.In the light of biological meaningfulness, Model (1) will only be studied in the following domain For convenience, we denote In the rest of this paper, we first present some theoretical results on the dynamics of Model (1) in Section 2 with proofs given in the Appendix.In Section 3, we present some numerical simulations for Model (1) to illustrate our theoretical results.We then discuss some ecological implications of our results in Section 4 and give a brief summary in the last section.
2. Mathematical results.Our first result presented below concerns the existence of possible equilibria.
About the stability of the equilibria, we have the following result.
Further, we have the following result on the existence of a limit cycle.
Remark 1.The above results show that, for Model (1), as the model parameters vary, the number and stability of equilibria change, and a limit cycle may appear.In other words, bifurcations occur at some critical values.Using p as the bifurcation parameter, we have the following possible bifurcation scenarios.
3. Numerical simulations.To illustrate our main results and the above bifurcation scenarios, in this section, we numerically sketch two representative bifurcation diagrams.To this end, we take parameter values comparable to values in the literature.For simulation purposes, all parameters in Model ( 1) are chosen on the basis of 1 ha of alpine meadow and we use year as the time unit.The annual per-capita death rate of plateau pika r 2 has been reported as 0.9476 [18] or 0.9674 [17].Here we take r 2 = 0.95year −1 for our simulations.Parameter r 1 is the intrinsic per-capita growth rate of vegetation.According to [7], the average value of the growth rate of fenced non-degraded, lightly degraded, moderately degraded and heavily degraded alpine meadows is 0.227year −1 .The value of K is the carrying capacity of 1 ha of alpine meadow, that is, the maximum possible yield.Reported yields are about 3536kg [8], 3486kg [7], 2200 ∼ 3000 kg [1], and approximately 4000 kg [21].Here we take K = 4000 kg in simulations.According to the data from non-degraded, lightly degraded, and moderately degraded alpine meadows in the study of Liu et al. [9], the area of one hole and one mound is 0.18m 2 and the hole coefficient (i.e. the amount of hole the amount of plateau pika ) is 4.48.Hence, one plateau pika makes barren 0.8064m 2 of alpine meadow, and the corresponding reduced yield λ is 4000×0.806410000 = 0.32256kg • head −1 .On average, the per-capita feeding rate q is 0.36kg Parameter p is the maximum per-capita growth rate of plateau pika in one year, which is equal to the product of litter size, number of reproduction in one year, the proportion of breeding females among all females, and the proportion of females in all plateau pika.According to [11,13,16,23], the value of p would be in the range of 1 ∼ 12year −1 .As no measured values for α and µ are available in the literature, we will vary the values of α (in the same order as K) and µ to demonstrate the possible occurrence of those bifurcation scenarios mentioned in Remark 1.The default values of parameters are summarized in Table 1.
Here we only present two representative bifurcation diagrams associated with S-I and S-III, respectively.In these figures, equilibria or limit cycles represented by solid curves are stable, while dotted curves represent unstable equilibria.
4. Ecological implications.In the alpine meadow ecosystem, both vegetation and plateau pika serve important ecological functions and they are supposed to coexist.In general, the plateau pika-vegetation system is stabilized at a stable positive equilibrium.When bifurcations occur, three phenomena may happen: (1) The positive equilibrium loses its stability, and the boundary equilibrium E 0 becomes stable.As a consequence, the plateau pika becomes extinct and vegetation approaches its  1) with r 1 = 0.227, α = 1000, µ = 0.0001, r 2 = 0.95, K = 4000, λ = 0.3225, q = 0.36 and p varying from 1.2 to 2.2.The solid curves represent stable equilibria or limit cycle, while the dotted curves represent unstable equilibria.
carrying capacity; (2) The positive equilibrium becomes unstable and a limit cycle appears.Thus the amount of vegetation and plateau pika vary periodically; (3) Both the positive equilibrium and the boundary equilibrium are stable.
In the case that E 2 is unstable, a limit cycle exists and both populations of vegetation and plateau pika undergo sustained oscillations.In the phase plane, the solution trajectories rotate around E 2 counterclockwise (see Figure 3).Depending on the initial conditions, the following four phases would appear in turn periodically: vegetation decreasing-plateau pika increasing; vegetation decreasingplateau pika decreasing; vegetation increasing-plateau pika decreasing; and vegetation increasing-plateau pika increasing.In this case, when either the vegetation population or the plateau pika population is at a minimum, the populations might be prone to extinction from catastrophic events [5].
In a healthy alpine meadow ecosystem, the amount of vegetation and the plateau pika population density should maintain at a constant level, that is, the vegetation and the plateau pika populations should be near an equilibrium rather than a limit cycle.In contrast to a healthy alpine meadow, on a degraded alpine meadow there is less vegetation and there are more plateau pika [24].
Many attempts have been made to identify the causes of alpine meadow degradation and to design an effective restoration strategy.The parameters in Model (1) were formed in long-term biological evolution.Anthropogenic activity or environment changes, however, could cause shifts in the parameter values, possibly leading to alpine meadow degradation.In the case that the equilibrium E 2 = (x 2 , y 2 ) is stable, direct calculations (taking partial derivatives) show that as p increases, or  1) are converging to a stable limit cycle.Parameter values used ere r 1 = 0.227, µ = 0.0001, r 2 = 0.95, K = 4000, λ = 0.3225, q = 0.36, α = 1000 and p = 2; Two sets of initial conditions are (3500, 1000) and (1500, 1000).
µ decrease, the value of x 2 decreases, while the value of y 2 increases.This results in more plateau pika and less vegetation.Plateau pika live in an extreme environment and have an r-life history strategy, and so global warming gives plateau pika an advantage.Although there is not much research on how climate influences plateau pika behavior, it is reasonable to assume that under global warming, the plateau pika death rate will decrease and birth rate will increase, that is, r 2 will decrease and p will increase.Therefore, global warming may be a reason that causes alpine meadow degradation.Higher vegetation would increase the chance of plateau pika to be caught by its predator.Decreased predation corresponds to a decrease in the parameter µ.In recent decades, poaching has caused a dramatic decrease of wild animals, including raptors and carnivorous mammals that are predators of plateau pika.Thus poaching reduces the value of µ, further contributing to the degradation of alpine meadows.
On a degraded alpine meadow, decreasing p, increasing r 2 , or increasing µ would increase vegetation, decrease plateau pika and hence may help restore the degraded alpine meadow.In practice, controlling plateau pika with sterilants can decrease p, the utilization of rodenticides can increase r 2 .Protecting natural enemies of plateau pika can enlarge µ, which can be achieved by prohibiting poaching, creating suitable habitats for these predators, or releasing man-raised natural enemies.5. Summary and discussion.In this paper, we have proposed a predator-prey model to understand the interplay between vegetation and plateau pika in an alpine meadow ecosystem.The dynamics of the model, including the existence and stability of equilibria, bifurcation of equilibria, and existence of limit cycles, have been analyzed.Based on our theoretical analysis, we can identify possible causes of alpine meadow degradation, and provide some insights on designing effective restoration strategies.Our model shows that global warming and poaching could explain observed alpine meadow degradation.Controlling plateau pika or protecting its natural enemy may be effective strategies for the restoration of degraded alpine meadow.
In contrast to the classic predator-prey model, our model incorporated the destruction function of plateau pika in alpine meadows (decreasing the area of vegetation growth) and the negative effect of the height of vegetation on the growth of plateau pika.Model (1) exhibits rich and complex dynamical behavior including the coexistence of two stable equilibria.
Although the causes of alpine meadow degradation and effective restoration strategies have been briefly discussed on the basis of our modeling exercise, more efforts are needed for the development of practical and feasible management policies.For an alpine meadow ecosystem, in addition to the vegetation and plateau pika, livestock, insects, raptors, and carnivorous mammals are all very important components to be included.Moreover, soil, climate, human activity, and social development can also greatly influence the stability of an alpine meadow ecosystem.
Acknowledgments.The authors would like to thank the editor and two anonymous reviewers for their very constructive comments, which greatly helped the improvement of the presentation of this work.1) is a solution to the following system in the domain G.
The existence of O and E 0 follow immediately.If y = 0, then it follows from the second equation of (3) that µx 2 − Σx + r 2 α = 0, which admits two solutions provided that ∆ > 0. If Σ > 0, then both x 1 and x 2 are larger than 0, if Σ < 0, then both x 1 and x 2 are smaller than 0. Therefore, when Σ > 0 and ∆ > 0, that is Σ > 2 √ r 2 αµ (or √ p > √ r 2 + √ αµ), both x 1 and x 2 exist and are positive.It then follows from the second equation of (3) that where x * is x 1 or x 2 .Note that the discriminant of the above quadratic equation is Equation ( 4) has two roots y ± = r1(α+x * ) 2q , which is not feasible, thus we only consider the root y − , which is positive if and only if K .Notice also that Kµ + r2α K ≥ 2 √ r 2 αµ and 2Kµ cannot lie between 2 √ r 2 αµ and Kµ + r2α K , otherwise, K 2 µ < r 2 α and K 2 µ > r 2 α would hold simultaneously.Thus the existence of E 1 and E 2 follows.In a similar way, we can obtain the existence of E 3 .
A.2. Proof of Theorem 2.2.At any equilibrium (x, y), the associated Jacobian matrix is given by At O, the Jacobian matrix is , whose two eigenvalues are −r 2 < 0 and r 1 > 0. Thus O is a saddle.At E 0 , the Jacobian matrix becomes , which has two eigenvalues −r 1 < 0 and −r 2 + pK α+K −Kµ.
Therefore, if Σ < Kµ + r2α K , both eigenvalues are negative, and E 0 is a stable node and when Σ > Kµ + r2α K , one eigenvalue is positive and the other one is negative, thus E 0 is a saddle.
At E 1 and E 2 , the Jacobian matrix reads as where x * and y * are the coordinates of So both eigenvalues of J * having negative real parts is equivalent to A =: pα (α+x * ) 2 − µ > 0 and B =: A contradiction and hence A is not larger than 0 for x 1 .This means that the equilibrium E 1 is unstable.
The second inequality in (b) implies Therefore, (a) and (b) imply that A > 0 always holds for x 2 .
A.3.Proof of Theorem 2.3.Under the stated conditions, the equilibrium E 1 does not exist and the equilibrium E 2 = (x 2 , y 2 ) with x 2 = Σ− √ ∆ 2µ , y 2 = h(x 2 ) is unstable.To use the Poincare-Bendixson theorem to prove the existence of a limit cycle, we construct a trapping region with its outer boundary labeled as OE 1 F GHO in Figure 5.In Figure 5, the line segment E 0 S is given by x = K and the line segment ST is given by y = K λ .Points O and E 0 are singular points and it is easy to know that OE 0 (0 < x < K) and T O(0 < y < K/λ) are piece of solution trajectories.On the line segment E 0 S, we have 0 < y < K λ , x = K and thus x < 0. This implies that any trajectory starting from a point on the line segment E 0 S would go leftward.We pick a point F on the line segment E 0 S with its y−coordinate in (y 2 , K λ ).Draw a horizontal line segment F I that intersects the line x = x 2 at I. On the line segment F I, we have y > 0, x 2 < x ≤ K.Note that in this case x 1 > K. Thus y = −µy α+x (x − x 2 )(x − x 1 ) > 0. This shows that the trajectory starting from the point F enters the rectangular domain F SJI and cannot leave this domain from the border F S or F I.Moreover, the trajectory starting from F cannot leave the domain F SJI from the border SJ because dx dy is not defined on SJ.Thus, the trajectory starting from F leaves the domain F SJI at some point G on the line segment IJ. Draw a horizontal line segment GH that intersects the y−axis at H. On the line segment GH, we have y > 0, x < x 2 < x 1 and hence y = −µy α+x (x − x 2 )(x − x 1 ) < 0, thus the trajectory starting from any point on GH would move downward.Therefore, if the stated conditions are satisfied, the region D 1 enclosed by the curve L consisting of the line segments GH, HO, OE 0 , E 0 F and the solution trajectory F G would be the required trapping region.Note that inside D 1 , there is a unique equilibrium E 2 , which is unstable.The well-known

Figure 5 .
Figure 5.The outer closed curve L: OE 0 F GHO around E 2 .

Table 1 .
The default values of parameters in simulations.