Dynamics of a Beddington–deangelis-type Predator–prey System with Constant Rate Harvesting

In the paper, a predator–prey system with Beddington–DeAngelis functional response and constant rate harvesting is considered. Various dynamical behaviors of the system including saddle–node points and a cusp of codimension 2 are investigated by using the analysis of qualitative method and bifurcation theory. Also, it is shown that the system undergoes several kinds of bifurcation such as the saddle–node bifur-cation, the subcritical and supercritical Hopf bifurcation, Bogdanov–Takens bifurcation by choosing the death rate of the predator and the harvesting rate of the prey as the bifurcation parameters. Some numerical examples are illustrated in order to substantiate our theoretical results. These results unveil far richer dynamics compared to the system without harvesting.


Introduction
In population dynamics, a lot of mathematical models to investigate the relationship between preys and predators have been suggested by ecologists and mathematicians [2,16].In fact, there are well-known predator-prey models such as Lotka-Volterra models [12,13], Hollingtype models [7,17], Michaelis-Menten-type models (so called ratio-dependent predator-prey models) [1,5], Ivlev-type models [10] and so on.Especially, many researchers are interested in the following Beddington-DeAngelis predator-prey model ) − γxy βy + x + α , y (t) = −dy + δxy βy + x + α , (x(0 + ), y(0 + )) = (x 0 , y 0 ), Corresponding author.Email: hkbaek@cu.ac.krwhere x, y represent the population density of the prey and the predator at time t, respectively.Usually, K is called the carrying capacity of the prey.The constant r is called the intrinsic growth rate of the prey.The constants δ, d are the conversion rate and the death rate of the predator, respectively.The term βy measures the mutual interference between predators.The reason for the model is that if one take the parameter β close to 0 then the model can be regarded as a Holling-type II predator-prey model and also one take the parameter α close to 0 then the model can be though of a ratio-dependent predator-prey model.The dynamics of system (1.1) has been studied extensively in [3,4,6,8,9].Research on this system (1.1) revealed various dynamics such as deterministic extinction, existence of multiple attractors and stable limit cycles, etc.
For mathematical simplicity, we need to rewrite system (1.1) as a dimensionless system with following scaling Then we can obtain the following dimensionless system where a = K α , b = γ rβ , D = rd and e = δK rα are positive constants.From the point of view of human needs, the exploitation of biological resources and the harvesting of populations are commonly practiced in fishery, forestry, and wildlife management.Concerning the conservation for the long-term benefits of humanity, there is a wide range of interest in the use of bioeconomic models to gain insight into the scientific management of renewable resources like fisheries and forestries [18].For the reason, it is important to deal with mathematical models with the harvesting of populations.Thus, by taking into account the above reasons, in the paper, we will focus on the Beddington-DeAngelis-type predator-prey system with constant rate harvesting as follows; = −Dy + exy y + ax + 1 . (1. 3) The organization of this paper is as follows.In next section, we study the existence of the equilibria and their stability in the small neighborhood for system (1.3).In Section 3 we show that system (1.3) can display bifurcation phenomena such as a saddle-node bifurcation, supercritical and subcritical Hopf bifurcations and Bogdanov-Takens bifurcation for a = 3, b = 1.Finally, we give a conclusion in Section 5.

Equilibria and stabilities
In this section, we investigate dynamical behaviors of system (1.3) around equilibria.From biological point of view, we shall assume that the dynamics of system (1.3) is on the positive cone R 2  + = {(x, y) : x ≥ 0, y ≥ 0}.Also, we shall consider the biologically meaningful initial condition x(0) ≥ 0, y(0) ≥ 0. Since the right-hand side two equations of system (1.3) are continuous and Lipschitzian in the close first quadrant R 2  + , it is easy to see that solution of system (1.3) exists and is unique for positive initial conditions.Also, we can see that the x-axis is invariant under the flow.However, this is not the case on the y-axis.In fact, all solutions to the system (1.3) touching the y-axis cross out of the first quadrant due to the harvesting rate h of the prey.Thus the first quadrant is no longer positively invariant under the flow generated by system (1.3), which makes the qualitative analysis of system (1.3) difficult.
The Jacobian matrix of system (1.3) is defined by where x and y are coordinates of equilibria, respectively.Thus elementary calculations yield, for i = 1, 2, and Proof.From Lemma 2.1, system (1.3) has a unique equilibrium (x 0 , y 0 ) if h 0 < h = 1 4 .From (2.5), we get and hence the determinant of the matrix J(x 0 , y 0 ) is zero.Thus the equilibrium (x 0 , y 0 ) is degenerate.In order to determine the dynamics of system (1.3) in the neighborhood of the equilibrium (x 0 , y 0 ), we need to transform the equilibrium (x 0 , y 0 ) of system (1.3) to the origin and expand the right hand side of system (1.3) as a Taylor series.Then system (1.3) can be written as where P 1 (x, y) and Q 1 (x, y) are C ∞ functions of order at least four in (x, y).
In order to show the conclusion (i), suppose that the condition e = D(a + 2) holds.Then both eigenvalues of the matrix J(x 0 , y 0 ) are zero.We use the procedure used in [22] to reduce system via the normal form method. First, let τ = − b 2+a t in system (2.8).Then system (2.8) becomes x (τ) = y + P 2 (x, y), where P 2 (x, y) and Q 2 (x, y) are C ∞ functions of order at least two in (x, y).From y + P 2 (x, y) = 0, we can obtain the implicit function where R 1 (x) is a C ∞ function of order at least four.And we have where R 2 (x) is a C ∞ function of order at least five.Using the notation of Theorem 7.2 of Chapter 2 in [22], we can find k = 2m + 1 = 3, m = 1, and a k = 4e b 2 > 0. Thus the equilibrium (0, 0) of system is a saddle point.Therefore (x 0 , y 0 ) is a saddle if e = D(a + 2).
In order to prove the conclusion (ii), assume that e = D(a + 2) hold.In fact, since the trace of the matrix J(x 0 , y 0 ) is not zero, one of the eigenvalues of the matrix J(x 0 , y 0 ) is zero and the other is nonzero.By letting x = X − b 2+a Y and y = (−D + e 2+a )Y, system (2.8) can be written as where P 3 (X, Y) and Q 3 (X, Y) are C ∞ functions of order at least three in (X, Y).Again, let τ = (−D + e 2+a )t.Then dτ = (−D + e 2+a )dt and hence system (2.12) becomes where P 4 (X, Y) and Q 4 (X, Y) are C ∞ functions of order at least three in (X, Y).Now using Theorem 7.1 of Chapter 2 in [22], we can obtain that (x 0 , y 0 ) is attracting(repelling) if e < D(a + 2) (e > D(a + 2)).
(iii.1) (x 1 , y 1 ) and (x 2 , y 2 ) are hyperbolic saddles if h < h 1 ; (iii.2) (x 1 , y 1 ) is a hyperbolic saddle and (x 2 , y 2 ) is a hyperbolic stable node if h > h 1 and aD < e ≤ (a + 2)D; Proof.It is easy to show the existence of equilibria under given conditions due to Lemma 2.1.Thus we have only to consider their stabilities.It follows from (2.5) that the eigenvalues of the Jacobian matrix J at the equilibria (x 1 , y 1 ) and (x 2 , y 2 ) are 1 − 2x i and −D + ex i ax i +1 , the dynamics of system (1.3) depends on the sign of the eigenvalues −D Thus, by linear stability analysis, (x 1 , y 1 ) is a hyperbolic saddle and (x 2 , y 2 ) is a hyperbolic stable node.Now consider the following cases.
is negative.On the other hand, the sign of the eigenvalue E 2 depends on that of e − (a (ii) Note that the condition be ≥ abD + e implies e ≥ abD.
(ii.1) From cases (I) and (II), we know that the sign of the multiplication of the eigenvalues E 1 and E 2 is always negative if h < h 1 .Thus the equilibria (x 1 , y 1 ) and (x 2 , y 2 ) are hyperbolic saddles.
(ii.3) Similar to (ii.2), by using case (II), we have that (x 1 , y 1 ) is a hyperbolic unstable node and (x 2 , y 2 ) is a hyperbolic saddle if h > h 1 and (a + 2)D < e.
It is easy to see that the determinant of the matrix J(x * 0 , y * 0 ) is zero.Thus the equilibrium (x * 0 , y * 0 ) is degenerate and at least one of the eigenvalues of the matrix J(x * 0 , y * 0 ) is zero.In fact, from elementary calculation, we obtain the eigenvalues are E 1 = 0 and (2.14) In order to determine the dynamics of system (1.3) in the neighborhood of the equilibrium (x * 0 , y * 0 ), we need to transform the equilibrium (x * 0 , y * 0 ) of system (1.3) to the origin and expand the right hand side of system (1.3) as a Taylor series.Then system (1.3) can be written as where g 1 (x, y) and g 2 (x, y) are C ∞ functions of at least the third order with respect to (x, y) and e(−abD − e + be) 2  and (2.17) Now, assume that the condition of (i) holds.Then it follows from (2.14) that the eigenvalue E 2 is nonzero, which means that the trace of Consider C ∞ changes of coordinates in a small neighborhood of (0, 0) as follows; and (2.18) Then system (2.15) can be transformed into where k 1 (x 2 , y 2 ) and k 2 (x 2 , y 2 ) are C ∞ functions with at leat the third order with respect to (x 2 , y 2 ) and , where l 1 (x 2 , y 2 ) and l 2 (x 2 , y 2 ) are C ∞ functions with at leat the third order with respect to (x 2 , y 2 ).
) is a nonzero constant depending on the parameter (a, b, D, e, h) and l 1 (x 2 , 0) is a C ∞ function with at least the third order with respect to x 2 .Thus it follows from Theorem 7.1 of Chapter 2 in [22] that the equilibrium (x * 0 , y * 0 ) is a saddle-node.For the conclusion (ii), assume that a (2.15).Thus system (2.15) can be transformed into where K 3 (X, Y) and K 4 (X, Y) are C ∞ functions with at least the second order with respect to (X, Y).By taking 2 where 2 is C ∞ function with at leat the third order with respect to (X 1 , Y 1 ).Elementary calculations yield that H 2 = −2 and If H 1 = 0, then be = abD + e, which contradicts to the condition be < abD + e. Hence H 1 is a nonzero constant.Using the notation of Theorem 7.3 of Chapter 2 in [22], we can find k = 2m = 2, m = 1, and a k = H 1 = 0 and b 1 = H 2 = 0. Hence the equilibrium (0, 0) of system (1.3) is a cusp of codimension 2. Therefore (x * 0 , y * 0 ) is a cusp of codimension 2.
Theorem 2.5.If max{h 1 , bD e } < h < min{h 0 , 1 4 }, be < abD + e and (e − aD)(abD + e − be) > 2eD, then system (1.3) has four equilibria as shown in Lemma 2.1 and the dynamics of system (1.3) are as follows: (i) (x 1 , y 1 ) is a hyperbolic saddle and (x 2 , y 2 ) is a hyperbolic stable node if aD < e ≤ (a + 2)D; Proof.By Theorem 2.4, (i) and (ii) can be easily obtained.Next, let us think about the sign of the determinant of the matrix

Saddle-node bifurcation
From Lemma 2.1 and Theorem 2.2 and 2.3, we can easily see that is a saddle-node bifurcation surface.In fact, the number of equilibria of system (1.3) changes from zero to two, and the two equilibria which are axial equilibrium points are the hyperbolic saddle and node when the parameters pass through the surface SN.Ecologically speaking, the system collapses for h > 1 4 and the prey population becomes extinct if h = 1 4 , but the prey population does not go to extinction for some initial value if 0 < h < 1 4 .From Theorem 2.4 and 2.5, we can know that the surface is also saddle-node bifurcation surface.This saddle-node bifurcation yields two positive equilibrium points.This means that the predator goes either extinct or out of R 2 + in finite time when h > h 0 , and both the predator and prey populations coexist in the form of a positive interior equilibrium for some initial values when h < h 0 , be < abD + e, (e − aD)(abD

Bogdanov-Takens bifurcation
From Theorem 2.4, we can see that system (1.3) could have a positive equilibrium which is a saddle-node or a cusp.Especially, the part (ii) of Theorem 2.4 shows that the positive equilibrium (x * 0 , y * 0 ) is a cusp of codimension 2 under some conditions, which implies that there may exist a Bogdanov-Takens bifurcation.In order to discuss such bifurcation phenomena of system (1.3) we need to fix some values of parameters for system (1.3) since the condition given in Theorem 2.4 are too complicated to analyze the number of codimension of the unique positive equilibrium.For the reason, in this section, we fix a = 3, b = 1 to investigate a Bogdanov-Takens bifurcation of system (1.3) in the parameter space (D, e, h).Although we will deal with the special case a = 3, b = 1, the procedure of investigating for bifurcations of system (1.3) can be used to discuss bifurcations for general cases.
It follows from Theorem 2.4 that (x * 0 , .Thus there exists a parameter space In order to show that system (1.3) undergoes a Bogdanov-Takens bifurcation the values D and h can be chosen as bifurcation parameters.We need to find the universal unfolding of (x * 0 , y * 0 ).To do this let h = h 0 − µ 1 and D = D 0 − µ 2 in system (1.3),where µ 1 and µ 2 are very small parameters.Then for (D, e, h) ∈ BT we can have the following system where D 0 = e(e−3)+e √ e 2 −78e+225 18(e−3) , h 0 = D 0 (9D 0 +4e) 4e 2   and e > 9D 0 .Translating (x * 0 , y * 0 ) to the origin by letting x 1 = x − x * 0 and y 1 = y − y * 0 in (3.8).Using the Taylor expansion system (3.8) can be written as where p 1 and q 1 are C ∞ functions of x 1 , y 1 of order of at least three and , , and Consider C ∞ changes of coordinates in a small neighborhood of (0, 0) as follows; (3.10) Then system (3.9) can be written as x 3 (t) = y 3 + p 3 (x 3 , y 3 ), where p 3 and q 3 are C ∞ functions of x 3 , y 3 of order at least three and δ , Next, by letting s = t 1−δ 5 x 3 , x 4 = x 3 and y 4 = (1 − δ 5 x 3 )y 3 in (3.11) and rewriting s as t, system (3.11) can be transformed into x 4 (t) = y 4 + p 4 (x 4 , y 4 ), where p 4 and q 4 are C ∞ functions of x 4 , y 4 of order at least three and λ 0 = δ 0 , λ 1 = δ 1 − 2δ 0 δ 5 , Taking the change of variables , y 5 = y 4 + p 4 (x 4 , y 4 ), we obtain the following system by setting t as t; x 6 (t) = y 6 , where q 6 is a C ∞ function of x 6 , y 6 of order at least three and ω Therefore, we have the following theorem.Theorem 3.2.
3) undergoes a Bogdanov-Takens bifurcation.Hence there exist values of the parameters (D, e, h) such that system (1.3) has a unique unstable limit cycle for some parameter values, and system (1.3) has an unstable homoclinic loop for other parameter values.

Numerical examples
In this section, we will illustrate some numerical examples to substantiate the validity of our theoretical results obtained in the previous sections.Then we can easily check that the hypotheses of Theorem 2.5 are satisfied.Moreover, we can obtain (x * 2 , y * 2 ) ≈ (0.4067, 0.7081) and h * ≈ 0.1846.If the value of h is set to be as h * − 0.001 ≈

Conclusion and discussion
In summary, we have analyzed dynamical behaviors of the Beddington-DeAngelis predatorprey system with constant harvesting rate.We have shown that there exist some parameter regions in which both predator and prey species can be existed or extinct simultaneously.Also, it has been obtained some sufficient conditions for the existence of limit cycles of the system via Hopf bifurcation and Bogdanov-Takens bifurcation.Finally, we have displayed some numerical examples to testify our theoretical results.In this paper, following the work of L. Ling and W. Wang in [10], M. Haque in [6], J. Xia, Z. Liu, R. Yuan and S. Ruan in [18], D. Xiao and L. S. Jennings in [19], and D. Xiao, W. Li and M. Han in [20] we continue investing the harvesting effect on the dynamics of a predatorprey system with Beddington-DeAngelis functional response.As a result we figure out that harvesting may play a significant role in dynamics of a predator-prey system.
According to conditions of parameter values, we have investigated the existence of equilibrium points and the stabilities of them.In Lemma 2.1 we have considered parameter values in almost all cases, but some parameter values are not dealt with because they are very minor cases and their classification method is similar to Lemma 2.1 as shown in Lemma 5.1.+ bD e , and assume that be < abD + e and (e − aD)(e + abD − be) > 2eD.Then (i) system (1.3) has a unique equilibrium (x * 0 , Thus using Lemma 5.1 and mathematical methods used in this paper, one can acquire similar results to those obtained in the paper.However we will not mention it since it could be very routine and tedious. As mentioned in the abstract, by letting the parameters α = 0 and β = 0, and adding the harvesting rate on the prey species in system (1.1), we can obtain the following two kinds of predator-prey systems (5.1) and (5.2) which are called a ratio-dependent system with the harvesting rate and a Holling-type II system with the harvesting rate, respectively.(5.2) Also, in [6], in order to account for the influence of intra-species competition among the predator population in system (1.1), the author considers the following predator-prey system.Table 5.1:This table lists the stability and bifurcation of all systems (1.1), (1.3), (5.1), (5.2) and (5.3).
Although system (5.3) is a predator-prey system with competition effect on the predator, we may regard the system as a predator-prey system with harvesting term since the predator population could decrease due to the term hy 2 .Now we will compare the dynamics of these systems around equilibrium points.We will focus on a positive equilibrium state among equilibria since the existence of a positive equilibrium state could give rise to the coexistence of all species, which is a biologically significant fact from the point of view of biodiversity.In fact, all systems which mentioned above have at least one positive equilibrium like (x * i , y * i ), i = 0, 1, 2. Now we summarize the dynamical behavior of systems around positive equilibrium points in Table 5.1 by letting LAS, SNB, HB, TB and BTB stand for locally asymptotically stable, and saddle-node bifurcation, Hopf bifurcation, transcritical bifurcation and Bogdanov-Takens bifurcations, respectively.Comparing the first row to the second and fifth row of Table 5.1, we figure out that predator-prey systems with harvesting term exhibit much richer dynamics than the system without harvesting.From the second, third and fourth rows, we assert that predator-prey systems with harvesting rate on the prey species have similar dynamical behaviors around its positive equilibrium points even if their functional responses are different from each other.From Table 5.1, we conclude that the harvest of any species in a predator-prey system can make various dynamical behaviors regardless of functional responses.