The Dynamical Behavior of a Certain Predator-Prey System with Holling Type II Functional Response ()
1. 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 researchers, 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
(1a)
(1b)
subject to initial conditions
. Here, x and y are the prey and predator densities at time t, respectively. All the above positive constants have biological considerations. Parameter
denotes the intrinsic growth rate of the prey;
represents the carrying capacity of the environment; a is the half-saturation constant;
is the search efficiency of predator for prey;
and
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
named Holling type II functional response describes the functional response of the predator. The specific growth term
governs the increase
of the prey in the lack of predator. The square term
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.
2. Preliminaries
2.1. Invariant Sets
Denotes the first quadrant as
, and its closure is
. For practical biological consideration, system (1) is defined on the domain
and all the solutions are positive. It is easy to prove lemma 1. Thus, any trajectory starting from
cannot cross the coordinate axes.
Lemma 1. All the solutions of system (1) are non-negative with the non-negative initial value
, i.e.
is an invariant set.
We will give some cases about invariant sets for system (1) in supplementary to lemma 1.
Remark. a) If
and
, domain
is an invariant set; b) Domain
and
is not an invariant set, if
; c) If
or
, domain
(
) is not an invariant set.
2.2. Boundness
Theorem 1 All the solutions of system (1) are bounded with non-negative initial conditions
.
Proof. From the equation (1a) we have inequality
. With the help of [14], this implies that
. We denote the positive upper bound as
. A similar argument, from the Equation (1b) we have inequality
(2)
Similarly, this implies that
. This completes the proof and the system under consideration is dissipative.
Lemma 2. All the solutions of system (1) satisfy
(3)
Proof. Taking an auxiliary function
, it is obvious to see that
(4)
where
. By applying the theory of differential inequality [15], we obtain
(5)
which, upon letting
, yields
. This completes the proof.
Remark. It is clear to see that all the solutions of system (0) are uniformly bounded.
Theorem 2. If
, then the system (1) is Lagrange stable.
Proof. Taking the positive definite Lyapunov function
with infinitesimally small upper and large lower bounds, computing its Dini derivative along the trajectories of Equations (1), we have
(6)
Thus
is negative definite. This completes the proof.
2.3. Permanence
Theorem 3. If parameters of system (1) satisfy
(7a)
(7b)
where
, then system (1) is permanent.
Proof. By theorem (2) we have a positive upper bound
such that
(8)
From Equation (1a) we have
here
is defined in theorem 2. With the help of [14], this implies that
. For
, there must exit T, such that for any
we have
, thus
Similarly, this implies that
. This completes the proof and we get the permanence of the system.
3. 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:
,
,
, where
. For practical or biological considerations, we omit the singular point
. The point
is a desired equilibrium only if
. When
, the critical point
becomes
.
Then we make a special effort to derive the existence conditions of an interior equilibrium denoted by
, where it can be given from the following algebraic equations
(9a)
(9b)
From Equation (9a) we know
as
;
as
; from Equation (9b) we know
as
;
as
,
. Thus, if the following conditions
,
satisfy, an interior equilibrium
exists and meanwhile we have
,
, where
.
3.1. 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
(10)
For zero point
, from its Jacobian matrix
we have: a) If
,
is an asymptotically stable node; b) If
,
is a saddle.
For the critical value
,
is a higher order singular point. According to Frommer’s method, by using its Jacobian matrix and polar-coordinate-transformation:
, we derive
and characteristic equation
. For the first simple real root
of the characteristic equation, from the calculation
,
,
and
, we know that
is a characteristic direction. While for the second simple root
, it is clear to see an expansion
, i.e.
(odd number) and
has the opposite sign as
. Thus
is also a characteristic direction. By the sign of
, we know that there is only one orbit tending to the critical point
along the direction
in the normal regions of second type.
Taking the transformation
, we have following normal form (retain the symbol “
” for simplicity)
(11)
where
,
in the right hand side are analytical functions, with terms starting from second order. From center manifold
,
the series of function
with respect to y, and
(even number), the critical point
is a saddle node and the parabolic sector is on the right half (x,y)-plane.
For point
, from its Jacobian matrix
and existence condition we have: a) If
,
is an asymptotically stable node; b) If
,
is a saddle.
For the critical value
,
is a higher order singular point. By
using its Jacobian matrix and polar-coordinate-transformation, we derive
and characteristic equation
, where
The real simple root
is in
, it is clear to see an expansion
and
, i.e.
(odd number) and
has contrary sign with
. Thus
is indeed a characteristic direction. While
changes its sign in the small neighbourhood of
, thus we know that there is only one orbit tending to the critical point
along the direction
in the normal regions of second type. Another real simple root
is in
, it is clear to see that
and
i.e.
, thus
is also a characteristic direction.
Taking transformations
;
,
and
, we have likewise real Jordan normal form
(12)
in here
,
in the right hand side are analytical functions, with terms starting from second order. From the center manifold
, the series of function
with respect to v is derived, where
(13)
Combining with
(even number), the critical point
is a saddle node and the parabolic sector is on the left half (u, v)-plane.
For the interior equilibrium
, from its Jacobian matrix
and existence conditions, denoting discriminant
with notations
(14a)
(14b)
we have:
(a) If
and
(a1)
, then
is an asymptotically stable node; (a2)
, then
is a saddle;
(b) If
and
(b1)
, then
is a center or a focus; (b2)
, then
is a saddle;
(c) If
, then
is unstable and
(c1)
, then
is a node; (c2)
, then
is a focus; (c3)
and
, then
is a node; (c4)
and
, then
is a saddle.
Remark. If
,
is asymptotically stable.
3.2. The Special Case:
If the interior equilibrium
exists and
, i.e. the Jacobian matrix
has a pair of purely imaginary eigenvalues
, in here
, or
is a center of local linear approximation of system (1). We consider this critical case or the center and focus problem by the H. Poincare’s formal series method. With the application of linear transforms
and
(or
), the above system (1) transforms into following normal form
(15)
where higher order terms
and
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
, where
is a homogeneous polynomial with order k (
). Computing its derivative along the above system, we have its total derivative
(16)
where
and
are homogeneous polynomials in u, v of degree k (
). In order to get Fourier series for further analysis, we firstly take the polar-coordinate-transformation
and assume
(
), then following recurrence relations are derived by comparing coefficients:
From a complicated expression of function
or the integral
, we know that
is precisely a periodic function of period
. Thus the focus value
of
reads
(17)
where
,
and
for later convenience. When
, the point
is a fine focus of order 1. If
,
is an unstable fine focus; If
,
is a stable fine focus. When
, continuing this procedure to obtain
.
3.3. The Special Case:
In this section, we consider a critical case:
, i.e.
is a critical point of higher order and its Jacobian matrix with zero eigenvalue can be rewritten as
(18)
where positive parameter
. By using its Jacobian matrix and the polar-coordinate-transformation, we derive function
(19)
and its characteristic equation
(20)
It is clear to see that the real roots
in
respectively satisfy
Case I:
In this case, for
,
if and only if
. With some calculation, we derive
and series of function
at
:
, where
(21)
i.e.
(odd number),
and
have opposite sign. Thus
is a characteristic direction. While
changes its sign in the small neighbourhood of
, thus we know that there is only one orbit tending to the critical point
along the direction
in the normal regions of second type.
If
, then
, while
Thus
is also a characteristic direction.
Case II:
In this case, for
, we have
,
and expression of function
at
:
, where
(22)
If
, then
and
is not a special direction; If
, then
is a characteristic direction.
From now on, making a change of variables
,
, and taking transformations
,
and
. Substituting them into the system (1), we have following normal form
(23)
where
,
in the right hand side are analytical functions, with terms starting from second order. From center manifold
and series of function
with respect to v, where
(24)
one can judge the type of critical point
in this special case. When
, i.e.
, combining with
(even number), the critical point
is a saddle node. Furthermore, if the coefficient
(or < 0), then the parabolic sector is on the right (or left) half (u,v)-plane. If
, continuing the above series to obtain
(25)
If
,
is an unstable node; If
,
is a saddle. When
, continuing the above procedure to obtain
and using the same criteria of
.
3.4. The Special Case:
In this section, we consider a critical case:
, i.e.
is a critical point of higher order and its Jacobian matrix with two zero eigenvalues can be rewritten as
(26)
From the above section we know that its characteristic equation is
, and function
Similarly,
is a special direction, where
.
Making a change of variables and taking transformations
;
, we have following normal form
(27)
From center manifold
, we have series of function
:
and series of following function:
where coefficients
(28)
If
, i.e.
(even number),
, or
, then
is a degenerate singular point. If
, but
(29)
i.e.
(odd number),
, and combining a discriminant
, where
(30)
then
is a complicated singular point and its neighbourhood
consists of one hyperbolic sector and one elliptic sector, if
and
, or just
, then
is a center or a focus.
3.5. Global Stability
In this section, we will give following theorems to illustrate the global stability by constructing proper Lyapunov functions.
Theorem 4. If
, then the zero point
is globally asymptotical stable.
Proof. Taking an unbounded positive definite Lyapunov function in lemma 2, we know that
is negative definite, similarly. This completes the proof.
Remark. Taking an unbounded positive definite Lyapunov function
, where positive constant
, and computing its Dini derivative along orbits of Equations (1), we have
(31)
where auxiliary function
On account of the discriminant
of a quadratic function, it is obvious that
or
is negative definite. This also completes the proof.
Theorem 5. Assume that the equilibrium
exists and
, then
is globally asymptotical stable.
Proof. Taking an unbounded positive definite Lyapunov function
(32)
and computing its Dini derivative along orbits of Equations (1), we have
(33)
This completes the proof.
Theorem 6. Assume that the equilibrium
exists and
, then
is globally asymptotical stable.
Proof. Taking an unbounded positive definite Lyapunov function
with positive constant
(34)
and computing its Dini derivative along orbits of Equations (1), we have
,
, where auxiliary functions
Conditions in this theorem imply that discriminant of a quadratic function in
is
thus
, i.e.
. By using the Krasovskii’s theorem, the proof is completed.
Theorem 7. Assume that the equilibrium
exists and
, then
is globally asymptotical stable.
Proof. Taking an unbounded positive definite Lyapunov function
(35)
and computing its Dini derivative along orbits of Equations (1), we have
(36)
From the condition we know that
is negative definite. By using the Krasovskii’s theorem, the proof is completed.
3.6. Closed Orbits and Limit Cycles
As we see, if
, i.e. the point
does not exist, there are no closed orbits in
. The trace of Jacobian matrix
is non-positive if
whether the existence condition
holds or not. By using the Bendixson criteria, there are no closed orbits in
. Such closed orbits and limit cycles surrounding the existed point
cannot cross the x and y coordinate axes and are confined in an invariant domain
, where the upper bounds A, B are all sufficiently large. The point
is not a saddle point for index
on this occasion. In this section, we will construct proper Dulac functions to further illustrate the existence problem of closed orbits and limit cycles.
Theorem 8. If
or
, then for system (1), there are no closed orbits in
.
Proof. Taking the diffeomorphism
which preserving the orientation of time, the system (1) is topologically equivalent to following system [16] [17]
(37)
For notation simplicity, we still retain the symbols x, y, t. Define a Dulac function
, along the above system, we have
where the numerator
By using the Bendixson-Dulac criteria, the proof is completed.
Remark. (a) If
, taking Dulac function
, then there are no closed orbits; (b) If
, taking Dulac function
, then there are also no closed orbits.
Theorem 9. If
, and
, then for system (1), there are no closed orbits in
.
Proof. We can take another type of Dulac function
for our discussion. Obviously, along system (37), we have
(38)
where these coefficients
in the summation are
For a quadratic function
of variable x, its discriminant implies that it is non-positive. Since other coefficients preserved are all non-positive, this completes the proof.
Remark. For the system (0) and the above form Dulac function
, if following conditions hold:
, and
, then there are no closed orbits. This procedure can also derive our theorem while it is uncomplicated but tedious.
Corollary 1. If
is asymptotically stable, then combining with conditions in theorem 8 or 9,
is globally asymptotical stable.
Remark. See reference [18] for the above corollary.
Theorem 10. If
and the equilibrium point
is unstable, then system (1) has limit cycle(s).
Proof. We see that
with
. Thus the straight line
is an untangent line of the system (37). Taking the Dulac function
, where
is sufficiently large, and computing
along the orbits of Equations (37), we derive
where
. For this system, we can construct a Bendixson ring including unstable singular point
. By the Poincare-Bendixson theorem, system (37) has at least one limit cycle in the first quadrant. This completes the proof.
4. Bifurcation Analysis
4.1. Transcritical Bifurcation
Without loss of generality, denoting the system (1) as
(39)
Firstly, by fixing values of all the parameters and varying the bifurcation parameter
, we observed that the two equilibrium points
and
collide with each other as
cross the critical value
. Thus, there is a chance of bifurcation around this point
.
Theorem 11. The system (39) undergoes a transcritical bifurcation around
when the parameters satisfy
.
Proof. For the Jacobian matrix at
, we have eigenvectors v and w corresponding to matrix
and its transpose, respectively:
. Then the following transversality conditions are hold:
(40)
Clearly, one of the eigenvalues of the Jacobian matrix
is zero and the other is negative. From the Sotomayor’s theorem( see [19] [20]), the system (39) undergoes a transcritical bifurcation around
at
. This completes the proof.
Here we take a as bifurcation parameter for the consideration of
.
Theorem 12. Assume that the parameters satisfy conditions for the existence of
and
, then the system (39) undergoes a transcritical bifurcation around
when the parameters satisfy:
(41)
Proof. For the Jacobian matrix at
, we have eigenvectors v and w corresponding to matrix
and its transpose, respectively:
,
. Then the following transversality conditions are hold:
(42)
It is clear to see that one of the eigenvalues of the Jacobian matrix
is zero and the other is negative. From the Sotomayor’s theorem, the system (39) undergoes a transcritical bifurcation around
at
. This completes the proof.
4.2. Hopf Bifurcation
We choose parameter d as a bifurcation parameter. Assume that the parameters satisfy conditions for the existence of
and
(43)
Define
as a pair of purely imaginary eigenvalues of matrix
, where
. It is clear to see that, if following conditions hold at
:
,
, and transversality condition
(44)
for instance,
. Then
changes its stability through Hopf bifurcation threshold the critical value
.
Theorem 13. Assume that the parameters satisfy conditions for the existence of equilibrium
, the condition (43) and the transversality condition (44), then the system (1) undergoes a Hopf bifurcation around equilibrium
as parameter d passes through the bifurcation value
.
We will calculate the first Lyapunov number
at the point
of the system which is used to determine the stability of limit cycles around Hopf bifurcation. Therefore we first translate the point
to the origin by a non-singular linear transformation
,
, then the system (39) in power series around the origin(drop the hats for the sake of convenience as usual) is:
(45)
with coefficients
(46)
and smooth functions
and
. From [19] [21], we know that the first Lyapunov number for a planar system is given by
(47)
By computing the first Lyapunov number at the point
, we will know that the Hopf bifurcation is supercritical and limit cycle(s) is(are) stable if
; the Hopf bifurcation is subcritical and limit cycle(s) is(are) unstable if
. Hence we have to give the following numerical simulations in the next section, since the above expression is much too complicated.
5. 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
,
,
,
,
,
and
. Thus, when
, the equilibrium point
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
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
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
Figure 1. (a) Time series diagram of population x; (b) Time series diagram of population y; (c) Phase diagram of population x and y.
Figure 2. (a) Time series diagram of population x; (b) Time series diagram of population y; (c) Phase diagram of population x and y; (d) Enlarged phase diagram of population x and y.
phase diagram of the system (1). On the basis of theoretical derivation and numerical calculus, the critical threshold of the control parameter d is approximately 0.04910. When the value of the control parameter d is greater than 0.04910, the other solutions of the system (1) can gradually converge to the
Figure 3. (a) Time series diagram of population x; (b) Time series diagram of population y; (c) Phase diagram of population x and y; (d) Enlarged phase diagram of population x and y.
equilibrium point
, which means that the equilibrium point
is asymptotical stable; this result can be seen in Figure 2 with
. 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
. At the bifurcation parameter
, the first Lyapunov number
, thus the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable.
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.
Acknowledgements
We thank the Editor and the referee for their comments. This research is funded by the National Natural Science Foundation of China (Grant nos.31570364) and the National Key Research and Development Program of China (Grant No.2018YFE0103700). This support is greatly appreciated.