Bifurcation analysis in a singular Beddington-DeAngelis predator-prey model with two delays and nonlinear predator harvesting

Abstract: In this paper, a differential algebraic predator-prey model including two delays, Beddington-DeAngelis functional response and nonlinear predator harvesting is proposed. Without considering time delay, the existence of singularity induced bifurcation is analyzed by regarding economic interest as bifurcation parameter. In order to remove singularity induced bifurcation and stabilize the proposed system, state feedback controllers are designed in the case of zero and positive economic interest respectively. By the corresponding characteristic transcendental equation, the local stability of interior equilibrium and existence of Hopf bifurcation are discussed in the different case of two delays. By using normal form theory and center manifold theorem, properties of Hopf bifurcation are investigated. Numerical simulations are given to demonstrate our theoretical results.

The Beddington et al. [18,19] gave the follows functional response Here, l 1 represents that per predator population per time can eat the maximum number of prey population, l is a positive constant and denotes the effect of handling time for predators, and l 2 is positive and measures the magnitude of interference among predators, d represents the prey density where the attack rate is half-saturated.It is obvious that two cases are possible as following.One case is that, if l = 1, l 2 = 0 and d > 0, then it reduces to a Holling II functional response (or Michaelis-Menten functional response) [20].The other case is that, if l = 0, l 2 = 0 and d > 0, then it reduces to a linear mass-action functional response (or Holling I functional response) [21].Conforto et al. [16] considered a three-dimensional reaction-diffusion system incorporating the dynamics of handling and searching predators, and showed that its solutions converge when a small parameter tends to 0 towards the solutions of a reaction-cross diffusion system of predator-prey type involving a Holling-type II or Beddington-DeAngelis functional response.Employing the upper and lower solution method and comparison theory, Li et al. [12] got the sufficient conditions of the upper ultimate boundedness and permanence of this system which implies that impulse always changes the situation of survival for species, and obtained the conditions for the existence of unique globally stable positive periodic solution.
Generally speaking, the introduction of time lag in the mathematical models [22][23][24][25][26][27][28][29][30][31][32] tends to reflect the interaction and coexistence mechanism of population in the past.System with two delays is discussed by using the equivalent system with a single time delay in reference [30].Liu et al. [33] used a Markovian switching process to model the telephone noise in the environment, proposed a stochastic regime-switching predator-prey model with harvesting and distributed delays, obtained the sufficient and necessary conditions for the existence of an optimal harvesting policy, and gave the explicit forms of the optimal harvesting effort and the maximum of sustainable yield.
Biological resources in the predator-prey system tend to be harvested and sold in order to obtain economic profit.In general, three types of harvesting function have been studied in the literatures [32,[34][35][36][37]. (I) Constant harvesting function is where C is a suitable constant.
(II) Proportionate harvesting function is where q is the catchability coefficient.
(III) Nonlinear harvesting function is where m 1 , m 2 are suitable positive constants.
It is easy to find that there are several unrealistic features in the proportional harvesting, such as stochastic search for prey and unbounded linear harvesting.However, the nonlinear harvesting function (III) eliminates the above unrealistic features and satisfies that is to say, the nonlinear harvesting function shows saturation effects in terms of harvest effort and inventory abundance [38].
As we known that the harvested prey-predator systems have been focused by both theoretical and mathematical biologists [15,39,40].Chakraborty et al. [40] studied the global dynamics and bifurcation of the predator-prey with constant harvesting.Liu et al. [15] investigated Hopf bifurcation and center stability for a predator-prey biological economic model with linear prey harvesting.Liu et al. [39] discussed bifurcation in a prey-predator model with nonlinear predator harvesting, and obtained the conditions for Turing and Hopf bifurcation.However, little work has been done on dynamic effects of economic interest on prey-predator system with nonlinear predator harvesting and Beddington-DeAngelis functional response.
In 1954, the common-property resource economic theory was proposed by Gordon in [41], which investigated the dynamic effects of harvesting efforts on the ecosystem from an economic point of view.Thus, in order to study the economic interest of commercial harvesting, an equation is proposed: Recently, a number of delay differential algebraic systems were proposed to investigate the impact of commercial harvesting on prey-predator model [42,43].Liu et al. [44] investigated a delayed differential-algebraic system with double time delays, Holling type II functional response and linear commercial harvesting effort on predator population.By jointly using the normal form of differential algebraic system and the bifurcation theory, Li et al. [45] discussed the stability and bifurcations, and obtained richer dynamics of the bioeconomic differential algebraic predator-prey model with nonlinear prey harvesting.Due to the above analysis and discussion, we will extend the work in [44,45] by considering nonlinear predator harvesting, double delays and Beddington-DeAngelis functional response function into our bioeconomic differential algebraic predator-prey system in this paper.The aim of our work is to reveal the dynamical behavior of Beddington-DeAngelis predator-prey model with two delays and nonlinear predator harvesting, to obtain a reasonable profit and provide some guidance for the harvest of biological economy system.
The rest of the paper is organized as follows.In Section 2, a differential algebraic predator-prey model including two delays, Beddington-DeAngelis functional response and nonlinear predator harvesting is proposed.In Section 3, positive solutions are analyzed when m = 0 and m > 0, respectively.In Section 4, in the absence of time delay, the singularity induced bifurcation is discussed by using differential-algebraic system theory, what's more, state feedback controllers are designed to remove singular induced bifurcation.In the presence of time delay, by analyzing the corresponding characteristic transcendental equation, the local stability around the interior equilibrium are studied.Furthermore, the existence of Hopf bifurcation is investigated.Directions of Hopf bifurcation and stability of the bifurcating periodic solutions are also discussed by using normal form theory and center manifold theorem.Numerical simulations illustrate the effectiveness of the mathematical conclusions in Section 5. Discussions and conclusions are included in the last section.

Model formation
Leslie and Gower [46] proposed the classical predator-prey system, which is as follows ). (2.1) Here, X(T ), Y(T ) indicate prey and predator population at time T respectively, s 1 , s 2 and K are intrinsic growth rate of prey, intrinsic growth rate of predator and carrying capacity.l 1 is the maximum number of prey each predator can eat each time and l 2 is a semi saturation constant, and the predation rate is The carrying capacity of predator is proportional to the size of prey population.l 3 is the amount of prey needed to support a predator population in equilibrium.
Based on the fact that Beddington-DeAngelis functional response is more authentic [47], we introduce it in system (2.1).Thus, the new system is ( Considering nonlinear predator harvesting, system (2.2) is constructed as follows .
(2.3) By using following transformations where a, b, α, β, r, q, p, c, m 1 and m 2 are all positive constants.Simultaneously, an algebraic equation is also included due to the economic profit of harvesting.Based on Eq (1.2), we have here p is the price per unit harvested biomass, c is the cost per unit harvest, while m means the net economic revenue.
In general, the delay differential equation model can produce more efficient and accurate dynamics than the ordinary differential equation model as capturing oscillation dynamics.
Therefore, by combining the above biological economic algebraic equation and time delay, a differential-algebraic system with two time delays is given as follows where τ 1 is maturation time for prey population, τ 2 is gestation delay for predator population.The initial conditions of system (2.6) take the following form: System (2.6) can be transformed into matrix form as follows where x(t−τ 2 ) .
Remark 1.Compared with system proposed in [44], nonlinear predator harvesting and Beddington-DeAngelis functional response are considered in system (2.6).
Remark 2. Compared with system proposed in [45], system (2.6) contains two delays and Beddington-DeAngelis functional response, and focuses on economic interest of commercial harvest effort on predator.
In the case of m > 0, we suppose (H4): pq(1 If one of the assumptions (H 4 ) and (H 5 ) holds, then the interior equilibrium P * = (x * , y * , E * ) exists, here 1+bx * −b and x * satisfies the following equation here

System without delays
Here, we are interested in the stability of system (2.6) at the interior equilibrium P 0 .
When τ 1 = 0, τ 2 = 0, system (2.6) takes the following form Lemma 4.1.(Singularity induced bifurcation theorem [49,50]).If the differential-algebraic equations satisfy the following conditions at the singular equilibrium: (I) D y f 3 has a simple zero eigenvalue and trace is also nonsingular, here ∆ = D E f 3 , then the singularity induced bifurcation occurs at the singular equilibrium.
Theorem 4.2.A singularity induced bifurcation takes place at the interior equilibrium P 0 of the differential algebraic system (4.1).When the bifurcation parameter m increases through zero, system (4.1) is unstable at P 0 .
Proof.We can obtain that the Jacobin matrix of system (4.1)around P 0 is Let m be bifurcation parameter, D be differential operator.We can obtain the following results.
(II) If J 11 0 holds, then it can be obtained that det Hence a singularity induced bifurcation occurs around P 0 of system (4.1) and the bifurcation value is m = 0. On the other hand, where According to Lemma 4.1 and Theorem 3 in [49], when m increases through 0, one eigenvalue of system (4.1)moves from C − to C + along real axis by diverging through ∞.Hence, the Theorem 4.2 is proved.
By using the similar proof of Theorem 4.2, it is easy to show system (4.1) is unstable around P * , what's more, state feedback controllers are designed to remove singularity induced bifurcation and stabilize system (4.1)around P 0 and P * , respectively.
In the case of m = 0, according to the leading matrix Ā in (4.1) and J P 0 , we can obtain rank (J P 0 , ĀJ P 0 , Ā2 J P 0 ) = 3.It is not hard to find that system (4.1) is locally controllable around P 0 based on Theorem 2-2.1 in [48].Therefore, a feedback controller can be used to stabilize system (4.1)around P 0 .By using Theorem 3-1.2 in [48], a feedback controller is as follows where k is a feedback gain.
Applying the controller into system (4.1),we can obtain that a controlled system is Proof.At first, the Jacobin matrix of system (4.2) around P 0 is The characteristic equation of system (4.2) around P 0 is det λ Ā − ĴP 0 = 0 based on Ā in system (4.2) and ĴP 0 , which can be written as follows where Similarly, in the case of m > 0, a controlled system is where k is a feedback gain.By using the similar analysis in Theorem 4.3, we can obtain the following result.
Remark 3.According to design of feedback controller, we can make interior equilibria be stable, which shows that prey-predator ecosystem can be kept sustainable and economic interest can be kept ideal by controlling commercial harvest effort on predator.

Case II:
When τ 1 = 0 and τ 2 > 0, Eq (4.4) can be written the following form Similarly, it shows that λ = 0 is not a root of Eq (4.9).We suppose that λ = iβ 2 ( β 2 > 0) is a root of Eq (4.9).We can obtain the following two transcendental equations which gives that By computing Eq (4.10), we can obtain τ 2k is 12) The proof of Theorem 4.6 is similar to that of Theorem 4.5.So, we omit it.

Case III:
In this part, τ 1 is considered as a parameter while τ 2 is regarded as a fixed value τ2 ∈ (0, τ 20 ), which is a stable interval calculated in Subsection 4.2.2.Here Ω 6 is defined as follows Let λ = iα 1 ( α 1 is a positive real number) represent a purely imaginary root of Eq (4.4).By separating real and imaginary parts, we have the following two transcendental equations Based on Eq (4.13), it derives that where Due to its complicated form, it is not easy for us to analyze properties of roots of transcendental equation (4.15).Based on dynamical system theory [1], we know that if and only if every eigenvalue has negative real part, system (2.6) is locally asymptotically stable around P * .What's more, by analyzing existence of Hopf bifurcation around the corresponding interior equilibrium P * , the periodic oscillation of system (2.6) is investigated.Hale [51] proposed that when the corresponding eigenvalue has a pair of purely imaginary roots, system usually exhibits Hopf bifurcation.Obviously, if Eq (4.15) has finite positive and simple roots 0 < α 10 < α 11 < • • • < α 1n , Eq (4.4) has a pair of purely imaginary roots.Without loss of generality, we denote α 1c = max {α 1k } , k = 0, 1, 2, • • •n and regard τ 1 as the bifurcation parameter, while we have the following corresponding critical value τ 1c here, ω 1c ∈ [0, 2π] satisfies the following equations cos ω 1c = l 10 (α 1c , τ2 ) l 12 (α 1c , τ2 ) , sin ω 1c = l 11 (α 1c , τ2 ) l 12 (α 1c , τ2 ) .(4.17) Furthermore, it is important to check transversal condition.The necessary condition for the existence of Hopf branches is that the eigenvalues passes through the imaginary axis with non-zero speed [51].
By differentiating λ with respect to τ 1 in Eq (4.4), it derives that where .
By virtue of Eq (4.18), it can be obtained that Θ = sign Re( dλ dτ 1 ) −1 where It is obvious to show Θ > 0 when B 11 B 13 + B 12 B 14 > 0. Therefore, we have the following results on stability and bifurcation in system (2.6).

Properties of Hopf bifurcation
In this part, we shall study the direction of the Hopf bifurcation and the stability of bifurcating periodic solution of system (2.6) when τ 1 is regarded as a parameter, τ 2 = τ2 is a fixed value.Similarly, we can discuss other cases.The approach employed here is the normal form method and center manifold theorem introduced by Hassard et al. [52] and Guckenheimer et al. [53].It follows from implicit function theorem [53] and the third equation of system (2.6) that E(t) = mm 2 y(t) q(py(t)−c)−mm 1 .Hence, system (2.6) can be transformed as follows Some transformations associated with P * (x * , y * ) are given as follows: Here, we define τ1 = µ + τ 1c , then µ = 0 is the Hopf bifurcation value of system (2.6).By simplifying, system (2.6) can be transformed to the following functional differential equation that is in the Banach space of continuous functions mapping C = C([−τ, 0], R 2 ) (τ = max {τ 1c , τ2 }), here τ 1c , τ2 are defined in (4.16) and (4.21), respectively, here, Based to Riesz representation theorem [54], there is a 2 × 2 matrix function η(θ, µ), here, θ ∈ [−τ, 0] and where here, δ denotes the Dirac delta function.
where η(θ) = η(θ, 0).From the above analysis, we know that A and A * are adjoint operators.According to discussion in Subsection 4.2, ±iα 1c are eigenvalues of both A and A * .
Based on q(θ) = (1, χ) T e iα 1c τ 1c θ , we have where ) is a constant vector.Similarly, based on (4.36) and (4.39), we have here, 2 ) is a constant vector.By above definitions and conditions, we have By above analyses and the results of Kuang [54], the following results can be given: By computing Eq (4.40), we have the following theorem.
(iii) The sign of T 2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T 2 > 0(T 2 < 0).

Conclusions
It is well known that commercial harvesting and economic benefits have a strong impact on the dynamical behavior.Liu et al. [44] investigated a differential-algebraic prey-predator system with linear harvesting on predator and Holling-II.Li et al. [45] analyzed a differential-algebraic prey-predator system without time delay.However, nonlinear harvesting and Beddington-DeAngelis functional response more realistic.Therefore, this paper proposed a singular Beddington-DeAngelis predator-prey model with two delays and nonlinear predator harvesting by extending the work of references [44] and [45], and obtained some results.The existence of equilibria was analyzed.Without considering time delay, the existence of singularity induced bifurcation by regarding economic interest as bifurcation parameter was discussed.In order to remove singularity induced bifurcation and stabilize system (4.1),state feedback controllers u(t) = (E(t) − 0.008) (see Figure 1b) and ũ(t) = (E(t) − 0.01) ( see Figure 2b) were designed, which shows that the system can be kept in a stable state with benefits by capturing predators.While considering time delay, stability of system were discussed by analyzing the corresponding characteristic transcendental equation.When τ 2 = 0, the critical value of time delay is τ 10 = 2.4; when τ 1 = 0, the critical value of time delay is τ 20 = 1.65; when τ 1 is regarded as a parameter and τ 2 as a fixed value, the critical value of time delay is    τ 1c = 2.12.At the same time, when τ 2 is regarded as a parameter and τ 1 as a fixed value, the critical value of time delay is τ 2c = 81.It was obvious that system lost local stability around it's corresponding interior equilibrium when time delays crossed corresponding the critical values, and Hopf bifurcations occurred (see Figures 4,5b,7,and 9).Finally, by using normal form theory and center manifold theorem, Hopf bifurcation is subcritical and the bifurcating periodic solutions exist for τ 1 < τ 1c , bifurcating periodic solutions are unstable and the period decreases because of ε 2 < 0, T 2 < 0, which could be found in Theorem 4.9.
In fact, the prey and predator may be captured simultaneously in real world.Thus, in order to make this model more practical, nonlinear predator harvesting and nonlinear prey harvesting can be    introduced into our predator-prey system, which is