Complexity Analysis of a Modified Predator-Prey System with Beddington–DeAngelis Functional Response and Allee-Like Effect on Predator

In this paper, complex dynamical behaviors of a predator-prey system with the Beddington–DeAngelis functional response and the Allee-like eﬀect on predator were studied by qualitative analysis and numerical simulations. Theoretical derivations have given some suﬃcient and threshold conditions to guarantee the occurrence of transcritical, saddle-node, pitchfork, and nondegenerate Hopf bifurcations. Computer simulations have veriﬁed the feasibility and eﬀectiveness of the theoretical results. In short, we hope that these works could provide a theoretical basis for future research of complexity in more predator-prey ecosystems.


Introduction
In reference [1], the authors simply considered a predatorprey system with Holling type II functional response and Allee-like effect on predator, which is described by the following nonlinear ordinary differential equations (ODEs): _ y � r 2 y 1 − y K 2 y y + e + e 1 qxy a + x − m 2 y, subject to initial conditions x(0), y(0) ≥ 0. Here we replace the Holling type II functional response (qx)/(a + x) with a functional response (qx)/(a + bx + cy) and denote the parameter q as q 1 for later use, and thus above system (1a) and system (1b) have a modified version: _ y � r 2 y 1 − y K 2 y y + e + e 1 q 1 xy a + bx + cy − m 2 y, where functions x � x(t) and y � y(t) are the densities of prey and predator at time t, respectively. In terms of biology, all above positive constants have practical considerations. Parameters r 1 and r 2 denote the intrinsic growth rate of the prey and predator, respectively; K 1 and K 2 represent the carrying capacity of the environment; a is the half-saturation constant; q 1 is the search efficiency of predator for prey; m 1 and m 2 are the mortality rate of the prey and predator species, respectively; e 1 is the biomass conversion and we denote e 1 q 1 as q 2 for convenience; d is the intraspecific competition coefficient of the prey; e is the Allee effect constant. e specific growth term r 1 x(1 − (x/K 1 )) governs the increase of prey in the lack of predator, while the specific growth term r 2 y(1 − (y/K 2 )) governs the increase of predator. e square term dx 2 is an intrinsic decrease term on prey. e term y/(y + e) on the specific growth term of predator with multiply form is called the Allee-like effect [2,3] and is different from predator-prey systems in [4][5][6][7][8] with Allee effect on prey. e coupled term q 1 x/(a + bx + cy) is named Beddington-DeAngelis (B-D) functional response (named after Beddington and DeAngelis et al.) [9,10]. It is similar to the Holling type II functional response incorporating an extra term cy in denominator, which describes mutual interference among predators [11,12]. is functional response has some of the same qualitative features as the ratio-dependent form but avoids some of singular behaviors of ratio-dependent models at low densities [11]. Obviously, when b � 1 and c � 0, system (2a) and system (2b) reduce to original system (1a) and system (1b). When parameters d � m 1 � e � 0 (without Allee effect) and r 2 − m 2 < 0, system (2a) and system (2b) reduce to System (2.1) in [13] and the author particularly conducted stability (local and global) and bifurcation (saddle-node, transcritical, Hopf-Andronov, and Bogdanov-Takens) analysis with a detailed mathematical analysis. When d � m 1 � r 2 � 0, system (2a) and system (2b) become model system (3a) and model system (3b) in [14], which was also independently and originally proposed in [9][10][11], while in [14], the authors discussed local and global asymptotic stability behavior of various equilibria and Hopf bifurcation occurs when parameter m corresponding to reserved region crosses some critical values. To mimic the real-world scenario, they solved the inverse problem of estimation of system parameter m by using the sampled data. System (1.3) in [15] is similar to above system in [14] except the constant rate harvesting term. In this reference, the authors showed that it undergoes several kinds of bifurcations, such as the saddle-node bifurcation, the subcritical and supercritical Hopf bifurcation, and Bogdanov-Takens bifurcation by choosing the death rate of the predator and the harvesting rate of the prey as bifurcation parameters.
Motivated by previous progress of predator-prey systems with B-D functional response or Allee-like effect, this paper mainly concentrates on dynamical analysis of system (2a) and system (2b). e rest of this paper is structured as follows. Preliminaries, such as boundedness, permanence, and existence of trivial equilibria, are given in Section 2. e existence of interior equilibrium is presented in Section 3 by virtues of the cobweb model and polynomial equations, respectively. In Section 4, we give stability analysis of equilibria and nonexistence of limit cycles. In Section 5, local codimension one bifurcations are analyzed, especially the Hopf bifurcation incorporating numerical simulations and Hopf bifurcation curves. In Section 6, we carry out short conclusions for our system.

Preliminaries
In this section, we devote to give some priori foundations. Before presenting the main results, we denote the first quadrant as R +2 , and its closure is R 2 + � R +2 . For biological consideration, system (2a) and system (2b) are defined on the domain R 2 + and all the solutions are non-negative with initial conditions x(0), y(0) ≥ 0, i.e., R 2 + is an invariant set and any orbits starting from it cannot cross the coordinate axes. Furthermore, all solutions are uniformly bounded. But now we only need to prove following theorems.

Boundedness and Permanence
Theorem 1 (uniform boundedness). Suppose that a nonnegative function φ(x, y) and its partial derivatives φ x and φ y are continuous in R 2 + , then the system subject to x(0), y(0) ≥ 0 is uniformly bounded.
is eorem 1 holds obviously after introducing an auxiliary function z � e 1 x + y [1,16]. e following theorem with the help of comparison principle in ODEs is about permanence of system (2a) and system (2b).
Proof. From above eorem 1, we have a positive upper bound ξ 1 such that en, we obtain M 2 > 0 and sufficiently large T 1 , such that By using the lemmas in [17], we admit liminf t⟶∞ x ≥ ω 1 . at is to say, there exist sufficiently large T 2 , such that en, we consider equation (2b) again. It yields _ y ≥ f(y)y; here the function f(y) incorporating its Taylor expression is Here and below, we denote the interior equilibrium E * of system (2a) and system (2b) as (x * , y * ) or (s 1 , s 2 ). is equilibrium must satisfy the following coupled algebraic equations: and we denote the positive root of a quadratic equation y 2 + 2ey − K 2 e � 0 as y s 1 � − e + ������� e 2 + K 2 e for later use.

Cobweb Model.
Based on above approximate analysis and the cobweb model, some cases about the existence of the interior equilibrium E * will be illustrated when isoclines from equations (10a) and (10b) all fall in R 2 + . Case 1. If parameters satisfy then an interior equilibrium exists. Here, we take r 1 � 1, then an interior equilibrium exists. Here, we take r 1 � 1, then an interior equilibrium exists. Here, we take then an interior equilibrium exists. Here, we take then an interior equilibrium exists (see the second equilibrium E 7 in Case 2). Case 6. If y k does not exist and parameters satisfy then an interior equilibrium exists. Here, we take then an interior equilibrium exists. Here, we take and d � 0.11; an interior equilibrium is E 9 ≈ (4.178700, 2.09452).

Polynomial
Equations of x * and y * . From equation (10a), an expression of y is Substituting it into the equation (10b), a quintic algebraic equation of x * can be derived in the form of p(x) ≔ 5 k�0 a k x k � 0 (see coefficients a k in Appendix A.1). anks to Niels Henrik Abel and Evariste Gallois's ingenious works, quintic equations usually have no analytical form solutions. But we can make some special efforts to numerically derive positive roots for this equation where x 1 > 0 and a 0 < 0, then there is a positive root x * such that x * < x 1 .

Case 1. If parameters satisfy
then an interior equilibrium exists. Here, we take 491956, 0.016748), and the following lemma is verified.
a k x k is a polynomial with real coefficients, a n ≠ 0, n > 1. If a n a 0 < 0, then the equation f(x) � 0 has a positive root.

Proof
We only consider the special case a n > 0, and thus a 0 < 0. It is clear that the polynomial f(x) has a decomposition If we take a sufficiently large positive number X such that X > max 0≤k≤n− 1 (n|a k |/a n ) 1/(n− k) , then f(X) > 0, and thus we complete the proof. □ Lemma 2. Suppose that f(x) � n k�0 a k x k is a real polynomial, a n ≠ 0, n > 1. If there is a positive number x 0 > 0 such that a n f(x 0 ) < 0, then the equation f(x) � 0 has a positive root.
On the other hand, from equation (10b), we obtain an expression And a quintic algebraic equation of y * can be written in the form of q(y) ≔ 5 k�0 b k y k � 0 (see coefficients b k in Appendix A.2). Suppose that the denominator and numerator of expression (23) are all positive for some y * > 0; then, x * > 0. Notice that if b 5 > 0 and b 0 < 0, there must exist a positive root for the quintic equation q(y) � 0 (see Lemma 1).
Discrete Dynamics in Nature and Society is a root of the quadric equation in the denominator of (23), then an interior equilibrium exists. Here, we take Case 3. If parameters satisfy b 0 < 0, y 2 > 0, q(y 2 ) > 0, y u 1 > 0, and q(y u 1 ) > 0, then an interior equilibrium exists. Here, we take r 1 � 1, 12 , and an interior equilibrium is E 12 ≈ (0.067460, 0.695847).

Stability Analysis of System (2a) and System (2b)
In this section, we use the Routh-Hurwitz criterion and Perron's theorems to analyze local stability of above equilibria in their existence interval, respectively. A theorem about global asymptotic stability and a theorem about nonexistence of limit cycles are also considered.

Local Stability Analysis.
e Jacobian matrix of system (2a) and system (2b) takes the following form J � (J ij ) 2×2 , where four components are For the trivial equilibrium E 0 with r 1 ≠ m 1 and the axial equilibrium E 1 with (q 2 x 1 /(bx 1 + a)) ≠ m 2 , we omit their stability [1]. When r 1 � m 1 , the transformation τ � − m 2 t and equations (2a) and (2b) yield _ x � (((r 1 /K 1 ) + d)/m 2 )x 2 (we still use symbol t), and thus E 0 is a saddle node and the parabolic sector is on the right half plane.
In the case that two axial equilibria E (k) 2 (k � 1, 2) exist, the Jacobian matrices are Since 2 ) > 0, and we have the following: 2 is an unstable higherorder singular point if r 1 − m 1 � ((q 1 y 2 )/(a + cy 2 )). In the special case that two axial equilibria E (k) 2 (k � 1, 2) collide with each other, the Jacobian matrix is and thus E 2 is a higher-order singular point. If Discrete Dynamics in Nature and Society For the interior equilibrium E * , its Jacobian matrix is where (32) E * is a node; (c2) Δ * < 0, then E * is a focus; (c3) Δ * > 0 and A 2 > 0, then E * is a node; (c4) Δ * > 0 and A 2 < 0, then E * is a saddle point.
When A 2 � 0 but A 1 ≠ 0, E * is a stable (unstable) node if A 1 < 0 (A 1 > 0) (see eorem 7.1 in Zhifen Zhang's book [18] for more details). It is probable that E * has a cusp case of codimension at least 2 which ensures potential Bogdanov-Takens bifurcation when A 1 � A 2 � 0.

Global Asymptotic Stability.
Combining the stability conclusions of the point E 0 in above section, the positive definite Lyapunov function V � e 1 x + y ensures that E 0 is globally asymptotically stable if one of the following conditions holds: Furthermore, conditions r 1 < m 1 , r 2 + (q 2 /b) − m 2 ≤ 0 and eorem 1 deduce global asymptotical stability of E 0 pronto. If 2 , E 2 , and E * do not exist, and E 0 is unstable, then eorem 1 ensures that E 1 is globally asymptotically stable. For the further consideration of point E * , the following theorem explains its global asymptotic stability.
Theorem 3 (global asymptotic stability of E * ). If a unique interior equilibrium E * exists and parameters satisfy then E * is globally asymptotically stable.
Proof. Here we take an unbounded positive definite Lyapunov function with A � ((q 1 (a + bx * ))/(q 2 (a + cy * ))). Introducing new variables x � x − x * and y � y − y * , computing derivative along orbits of system (2a) and system (2b), we have It is obvious that (dV/dt)| (2) is negative definite. Consequently, the Lyapunov function V satisfies the asymptotic stability theorem in [19]. us, we complete the proof.
(41) us, θ � θ 2 actually shows a fact that trajectories enter the stable node along this direction.

Theorem 4 (nonexistence of limit cycles). For system (42a) and system (42b), if parameters satisfy
then there are no closed orbits and limit cycles in R +2 .
Proof. Here we take a Dulac function B(x, y) � 1/xy and calculate following partial derivative: where coefficients a 02 � − 2K 1 c m 2 − r 2 K 2 + ar 2 , and other unlisted coefficients are all nonpositive. us, we complete the proof.
Here we take over values of parameters from Section 4.2 but m 1 � m 2 � 2, and system (2a) and system (2b) do not have meaningful equilibria except a globally asymptotically stable node E 0 (the origin O) with the characteristic direction θ � 0 (the positive x-axis) since the characteristic equation Θ(θ) � (m 1 − m 2 − r 1 )sin(θ)cos(θ), and conditions in eorem 4 all hold. Hence, there are no closed orbits and limit cycles in this numerical case. In addition, system (2a) and system (2b) merely own a saddle point E 0 and a globally asymptotically stable node E 1 Discrete Dynamics in Nature and Society when we use parameters from Section 4.2 but set a � 2,

Local Bifurcations of System (2a) and System (2b)
In this section, we will give sufficient conditions to show the existence of saddle-node bifurcation, transcritical bifurcation, and Hopf bifurcation of system (2a) and system (2b). Firstly, we denote this system in the following form: for simplicity and convenience.

Saddle-Node
Bifurcation. e two trivial equilibria E (k) 2 (k � 1, 2) collide with each other and system (47) has a unique boundary equilibrium E 2 when y > 0, if r 2 > m 2 and y 1 � y 2 , or en, there is a chance of bifurcation around this higherorder singular point. Here we choose m 2 as a bifurcation parameter and select eigenvector v � 0 1 corresponding to the zero eigenvalue for matrix (30). e eigenvector corresponding to the zero eigenvalue for the transpose of matrix Suppose r 1 − m 1 − ((q 1 y 2 )/(a + cy 2 )) < 0, then the following transversality conditions hold: (50) us, we have following theorem by using Sotomayor's theorem [22,23].
For the special case there is a chance the system (47) undergoes a pitchfork bifurcation. We still use the bifurcation parameter r 1 and eigenvectors v and w; then, the fourth transversality condition is where 1 ≔ m 1 + ((q 1 y 1 )/(a + cy 1 )). us, we have another theorem by using Sotomayor's theorem [22,23].

Hopf Bifurcation.
In this section, we consider the Hopf bifurcation of system (2a) and system (2b). Here the point E * exists and we choose d as bifurcation parameter. Suppose (54) en, system (47) undergoes a Hopf bifurcation around the point E * with respect to the bifurcation parameter d.
We will calculate the first Lyapunov number σ at the point E * , which is used to determine the stability of limit cycles and Hopf bifurcation direction. e method and calculations of the first and second Lyapunov coefficients can be found in [22,24,25]. erefore, translating the point E * to the origin O � (0, 0) by a linear transformation (I): X � x − x * , Y � y − y * , system (47) in power series around the origin is where coefficients are a 02 � q 1 s 1 c bs 1 + a bs 1 + cs 2 + a 3 , a 11 � − a 2 + bs 1 + cs 2 a + 2s 1 s 2 bc q 1 and O(|X, Y| 4 ) stands for some smooth functions. After that, by using a transformation (II): u � Y, v � (− Xb 10 + Ya 10 )/β with β � ������ � A 2 (E * ) > 0, the above system becomes a standard form: (58) be two corresponding eigenvectors of a matrix A such that Aq � iβq, A ⊤ p � − iβp and <p, q > � 1; the operation <x, y > � x † y (x, y ∈ C) with the Hermitian transpose (upper symbol) † represents the usual inner product, and the Jacobian matrix is We should rewrite the functions in the right hand side of system (57) in the form of power series where the components of linear functions B and C are while ‖x‖ is the two-dimensional Euclidean norm of x. Define T c be the largest subspace invariant by the matrix A and the generalized eigen subspace corresponding to the pair of purely imaginary eigenvalues ±iβ, i.e., T c � span q, q . at is to say, for any element z ∈ T c , there must exist a linear expansion z � wq + wq. Now we can construct a twodimensional parameterized center manifold in which h jk ∈ C 2 and its complex conjugate is h jk � h kj . Combining equations (57) and H(w, w), we arrive at a complex equation without consideration: e "chart" w for the central manifold H should be extracted from following differential equation: So, substituting it into equation (63) and comparing coefficients of these w j w k , we recursively derive and an equality from coefficient of the cubic term w|w| 2 : Here letter I represents the unit matrix with rank 2. It is quite apparent that equality (66) admits a solution: At last, the first Lyapunov coefficient is calculated as us, the first Lyapunov number σ for the focus of planar system (57) is given by the formula σ � (6π/β)l 1 . Since above expression is much too complicated, we need to present some numerical simulations and figures around the point E * with computer simulations.
Theorem 8 (nondegenerate Hopf bifurcation). Assume that the equilibrium E * exists and parameters satisfy condition (54) and σ ≠ 0; then, system (47) undergoes a nondegenerate Hopf bifurcation around this equilibrium as parameter δ � d passes through the critical value d [H] . e Hopf bifurcation is supercritical (subcritical), the interior equilibrium E * is a multiple stable (unstable) focus with multiplicity one, and limit cycles are stable (unstable) if σ < 0 (σ > 0).

Remark 2.
e first Lyapunov coefficient l 1 can be defined by (1/2β)Re(G 21 ) at times. On occasion, there are times when our system (55) exists some values of parameters such that σ � 0 or the system may undergoes a degenerate Hopf bifurcation. Accompanied by a proper transformation, for planar system (55), the first Lyapunov number is given by the following formula [22]: Discrete Dynamics in Nature and Society 11 which is an advanced version of σ or l 1 . Finally, for the Hopf bifurcation, we numerically give following examples to simulate how the parameter d controls dynamical behavior of system (2a) and system (2b). Example 1. Consider Case 6 in Section 3.1 and eorem 4 again. Figure 2( is implies that the Hopf bifurcation occurs in system (2a) and system (2b) when d � d [H] 8 . e first Lyapunov coefficient l 1 ≈ − 0.122574, and thus the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable. is interior equilibrium is a multiple stable focus with multiplicity one. Besides, for a perturbed system with sufficiently small parameter vector δ � (δ 1 , δ 2 ) in a neighbourhood of the origin O in the parameter plane and Hopf bifurcation analysis with two bifurcation parameters d and m 2 , we let δ ≠ 0 and suppose that E * � (x * , y * ) is an interior equilibrium of above perturbed system, where y * � y 8 + w, |w| is sufficiently small, and Substituting it into A 1 and A 2 , the solutions δ 1 � δ 1 (w) and δ 2 � δ 2 (w) can be directly solved, which are written as the form up to third order: and hence the slope k � lim w⟶0 ((δ 2 (w))/( δ 1 (w) )) of approximation straight line of the Hopf bifurcation curve Hp (see Figure 2( As a matter of fact, we perceive the phenomenon that eorem 8 describes the special case once δ lies on the horizontal line δ 2 � 0. Similar to above example, system (2a) and system (2b) undergo a Hopf bifurcation when d passes through d [H] 12 . e first Lyapunov coefficient l 1 ≈ − 0.002729 is also found to be negative. us, the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable. e interior equilibrium is also a multiple stable focus with multiplicity one. Figure 6 is the Hopf bifurcation curve corresponding to system (69) with slope k ≈ − 0.094044 and Hp � δ|δ 1 ≈ − 40.550983w + 639.229323w 2 + O w 3 , δ 2 ≈ 3.813575w − 4.617377w 2 + O w 3 . (76)

Conclusions
In summary, we firstly considered stability analysis and bifurcations of system (2a) and system (2b) with B-D functional response and Allee-like effect, which is a modified version of a predator-prey system in [1]. e polynomial's method, derived from eliminants R y (f, g) and R x (f, g), can be extended to more complicated polynomial systems. Some conclusions are the same as reference [1], such as the uniform boundedness ( eorem 1), the existence of equilibria E * , the nonexistence ( eorem 4) of limit cycles, and so on. It is supposed that some methods and conclusions can be available in original system (1a) and system (1b), such as pitchfork bifurcation. Lemmas 1 and 2 are available in more complicated predator-prey systems. Some critical cases, such as r 1 − m 1 � ((q 1 y 1 )/(a + cy 1 )) and A 2 � 0, need to be researched further with the help of topologically equivalent systems or the "blow-up" method (horizontal and vertical blowups). Other parameters can also be considered as a bifurcation parameter δ in Hopf bifurcation, although it is described by Hopf bifurcation curve Hp, for instance δ � x * . Notice that eorem 2 (permanence) can be extended to its reaction-diffusive version [26,27]. Under the conditions or the discussion of eorem 3, the Turing instability of its corresponding reaction-diffusion system subject to the homogeneous Neumann boundary conditions [28]: