The Dynamical Behavior of a Certain Predator-Prey System with Holling Type II Functional Response

In this paper we analytically and numerically consider the dynamical behavior of a certain predator-prey system with Holling type II functional response, including local and global stability analysis, existence of limit cycles, transcritical and Hopf bifurcations. Mathematical theory derivation mainly focuses on the existence and stability of equilibrium point as well as threshold conditions for transcritical and Hopf bifurcation, which can in turn provide a theoretical support for numerical simulation. Numerical analysis indicates that theoretical derivation results are correct and feasible. In addition, it is successful to show that the dynamical behavior of this predator-prey system mainly depends on some critical parameters and mathematical relationships. All these results are expected to be meaningful in the study of the dynamic complexity of predatory ecosystem.


Introduction
In 1965, C.S. Holling proposed three kinds of functional responses for different kinds of species to model the phenomena of predation, which made the traditional non-autonomous Lotka-Volterra predator-prey system more realistic [1] [2] [3] [4]. These functional responses describe how predators transform the harvested prey into the growth of itself and were discussed by numbers of re-searchers, including the stability of equilibrium points, existence of Hopf bifurcation, limit cycles, homoclinic loops, and even catastrophe [5].
For the Rosenzweig-MacArthur Model (R-M model) or the predator-prey model with Holling type II functional response, [6] studied stability of the R-M model by using graphical method, [7] studied global stability of the R-M model.
In [8], conditions for an interior equilibrium are given, and the stability of this equilibrium is analyzed. Certain critical cases, some of which cannot occur in the usual model are also discussed. (If m = 1, q(y) = 0, the above system reduces to the "usual" system.) [9] investigated a predator-prey system for the global stability and existence of limit cycles; they proved that there exist at least two limit cycles by using qualitative analysis and the idea of Poincare-Bendixson theory.
In recent years, for a more complicated system with Holling type II functional response, Liu et al. [10] investigated a predator-prey model with Holling type II functional response incorporating a constant prey refuge. The authors studied the instability and global stability properties of the equilibrium points and the existence and uniqueness of limit cycle. Lv et al. [11] considered a model to describe the harvesting for the phytoplankton and zooplankton based on plausible toxic-phytoplankton-zooplankton systems. Shanbing Li et al. [12] studied a spatially heterogeneous predator-prey model where the interaction is governed by Holling II functional response. They showed that the degeneracy for the prey and predator has distinctly different effects on the coexistence states of the two species when intrinsic growth rate of the prey is above a certain critical value. P. D. N. Srinivasu et al. [13] considered a prey-predator model with Holling type II of predation and independent harvesting in either species. Their study showed that using the harvesting efforts as controls can break the cyclic behavior of the system and drive it to a required state. By introducing impulsive control strategy, Yongzhen Pei et al. [2] investigated the dynamics behaviors of one-prey multi-predator model with Holling type II functional response with the help of Floquet theorem and small amplitude perturbation method. They showed that multi-predator impulsive control strategy is more effective than the classical and single one. Motivated by the above works, especially the references [8] and [9], we investigate a certain predator-prey system with density-dependent predator specific death rate and predator mutual interference incorporating a square term, which is described by the following nonlinear ordinary differential equations Here, x and y are the prey and predator densities at time t, respectively. All the above positive constants have biological considerations. Parameter 1 r denotes the intrinsic growth rate of the prey; 1 half-saturation constant; α is the search efficiency of predator for prey; 1 m and 2 m are the mortality rate of the prey and predator species, respectively; e is the biomass conversion; d is the intra-specific competition coefficient. The term xy a x α + named Holling type II functional response describes the functional response of the predator. The specific growth term 1 governs the increase of the prey in the lack of predator. The square term 2 dy denotes intrinsic decrease of the predator. Such ordinary differential system of predator-prey populations is familiar to the Lotka-Volterra system in which populations have the addition of damping terms (or self inhabit) [9].
The main purpose and conclusions of this paper are to analyze equilibrium stability and bifurcations of the above system by using qualitative analysis and bifurcation theories. This paper is organized as follows. Preliminaries, such as non-negativity, boundness, permanence and invariant sets are given in Section 2.
In Section 3, we give existence and sufficient conditions for the stability analysis of equilibrium points. In critical cases of interior equilibrium, we describe sufficient criteria to illustrate focus-center problem and problem of higher-order equilibrium. Existence and non-existence conditions of limit cycle(s) are also presented in this section. In Section 4, analyses of transcritical and Hopf bifurcations are presented. Here we choose d as a Hopf bifurcation parameter to consider limit cycle(s) further and point out that the control parameter d how to affect the complexity of this system. In Section 5, we carry out numerical simulations and conclusions of our system.

Invariant Sets
Denotes the first quadrant as 2 R + , and its closure is 2 2 R R + + = . For practical biological consideration, system (1) is defined on the domain 2 R + and all the solutions are positive. It is easy to prove lemma 1. Thus, any trajectory starting from 2 R + cannot cross the coordinate axes. Lemma 1. All the solutions of system (1) are non-negative with the non-negative initial value We will give some cases about invariant sets for system (1) in supplementary to lemma 1.

Boundness
Theorem 1 All the solutions of system (1) are bounded with non-negative ini- Similarly, this implies that ( ) ( ) Proof. Taking an auxiliary function z ex y = + , it is obvious to see that ( ) By applying the theory of differential inequality [15], we obtain .
Thus D V + is negative definite. This completes the proof.

Permanence
Theorem 3. If parameters of system (1) satisfy ≥ . This completes the proof and we get the permanence of the system.

Equilibria
In this section we will discuss equilibria of the system (1) with their sufficient existence conditions and stability analysis.
It is obvious to see that our system has following trivial equilibria:

Equilibria Type
Here we will use the Routh-Hurwitz criteria and the Perron's theorems to analyze local stability and type of these equilibria by the nature of eigenvalues of Jacobian matrices around them, respectively. We observed that the Jacobian matrix of the system (1) is  , 2 E is an asymptotically stable node; b) If For the critical value is also a characteristic direction.
Taking transformations = in the right hand side are analytical functions, with terms starting from second order. From the center manifold Combining with 2 m = (even number), the critical point * E is a saddle node and the parabolic sector is on the left half (u, v)-plane.
For the interior equilibrium * E , from its Jacobian matrix ( ) * J E and existence conditions, denoting discriminant we have:

The Special Case:
, where higher order terms ( ) , P u v and ( ) , Q u v in the right hand side of the above system are all real analytical functions of u and v, with terms starting from second order.
Suppose that the above system has a first integral in the form of ( ) ( ) ( ) . Computing its derivative along the above system, we have its total derivative . In order to get Fourier series for further analysis, we firstly take the polar-coordinate-transformation cos , sin u r v r θ θ = = and assume ( ) , then following recurrence relations are derived by comparing coefficients: :

The Special Case:
In this section, we consider a critical case: 2 E is a critical point of higher order and its Jacobian matrix with zero eigenvalue can be rewritten as ( ) where positive parameter ( ) . By using its Jacobian matrix and the polar-coordinate-transformation, we derive function and its characteristic equation It is clear to see that the real roots 1,2 i.e.
If 3 0 a > , * E is an unstable node; If 3 0 a < , * E is a saddle. When 3 0 a = , continuing the above procedure to obtain 4 a and using the same criteria of 2 a . Journal of Applied Mathematics and Physics

The Special Case: 1 2 = 0 A A =
In this section, we consider a critical case: 1 2 0 A A = = , i.e. * E is a critical point of higher order and its Jacobian matrix with two zero eigenvalues can be rewritten as From the above section we know that its characteristic equation is

Global Stability
In this section, we will give following theorems to illustrate the global stability by constructing proper Lyapunov functions.
and computing its Dini derivative along orbits of Equations (1), we have From the condition we know that D V + is negative definite. By using the Krasovskii's theorem, the proof is completed.

Closed Orbits and Limit Cycles
For notation simplicity, we still retain the symbols x, y, t. Define a Dulac function ( ) 1 By using the Bendixson-Dulac criteria, the proof is completed.
Remark Proof. We can take another type of Dulac function e , 0 ux u > for our discussion. Obviously, along system (37), we have where these coefficients 30 20 10 , , , a a a  in the summation are  , then there are no closed orbits. This procedure can also derive our theorem while it is uncomplicated but tedious. Corollary 1. If * E is asymptotically stable, then combining with conditions in theorem 8 or 9, * E is globally asymptotical stable.
Remark. See reference [18] for the above corollary.
Theorem 10. If 2 x a > and the equilibrium point * E is unstable, then system (1) has limit cycle(s).
Proof. We see that  ) ( ) For this system, we can construct a Bendixson ring including unstable singular point * E . By the Poincare-Bendixson theorem, system (37) has at least one limit cycle in the first quadrant. This completes the proof.

Transcritical Bifurcation
Without loss of generality, denoting the system (1) as ( ) , .   Theorem 13. Assume that the parameters satisfy conditions for the existence of equilibrium * E , the condition (43) and the transversality condition (44), then the system (1) undergoes a Hopf bifurcation around equilibrium * E as parameter d passes through the bifurcation value [ ] H d .
We will calculate the first Lyapunov number σ at the point * E of the system which is used to determine the stability of limit cycles around Hopf bifurcation. Therefore we first translate the point * E to the origin by a non-singular linear transformation ( ) ( )   2  2  3  2  2  3  10  01  20  11  02  30  21  12  03  1   2  2  3  2  2  3  10  01  20  11  02  30  21  12  03  2 , , , , x a x a y a x a xy a y a x a x y a xy a y F x y   y b x b y b x b xy b y b x b x y b   . From [19] [21], we know that the first Lyapunov number for a planar system is given by By computing the first Lyapunov number at the point * E , we will know that the Hopf bifurcation is supercritical and limit cycle(s) is(are) stable if 0 σ < ; the Hopf bifurcation is subcritical and limit cycle(s) is(are) unstable if 0 σ > . Hence we have to give the following numerical simulations in the next section, since the above expression is much too complicated.

Numerical Simulations and Conclusions
In this section we numerically give simulations of dynamical behavior between the prey and predator. In order to verify the feasibility and correctness of the theoretical derivation results, we will give some numerical simulations with a control parameter d and take 1 0.6 , the equilibrium point * E is (2.898588, 2.753889), then we can verify that the values of these parameters can meet theorem 7, which means that the equilibrium point * E is globally asymptotical stable. The time series diagram and phase diagram of the system (1) can be seen in Figure 1. It is easy to see from Figure 1 that the population x and y are permanent and can approach a fixed value respectively. These results demonstrate that the system (1) is permanent and the equilibrium point * E is globally asymptotical stable.
To investigate in detail how the control parameter affects the dynamical behavior of the system (1), Figure 2 and Figure 3 depict the time series diagram and   equilibrium point * E , which means that the equilibrium point * E is asymptotical stable; this result can be seen in Figure 2 with 0.04915 d = .
When the value of the control parameter d is less than 0.04910, the other solutions of the system can gradually converge to a limit cycle, which implies that the Hopf bifurcation occurs in the system (1), and this result can be seen in Figure 3 with Based on mathematical theory derivation and numerical simulation analysis, it is successful to show that dynamical behavior of this certain predator-prey system mainly depends on some critical parameters and relationships. Then, it is particularly significant to point out how the control parameter d affects the complexity of the system (1), which can lead this predator-prey system to emerge Hopf bifurcation. Moreover, these results also show that some control parameters can directly or indirectly affect the dynamic complexity of our predator-prey system.