Abstract
In this paper, complex dynamical behaviors of a predator-prey system with the Beddington–DeAngelis functional response and the Allee-like effect on predator were studied by qualitative analysis and numerical simulations. Theoretical derivations have given some sufficient and threshold conditions to guarantee the occurrence of transcritical, saddle-node, pitchfork, and nondegenerate Hopf bifurcations. Computer simulations have verified the feasibility and effectiveness 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.
1. Introduction
In reference [1], the authors simply considered a predator-prey system with Holling type II functional response and Allee-like effect on predator, which is described by the following nonlinear ordinary differential equations (ODEs):subject to initial conditions . Here we replace the Holling type II functional response with a functional response and denote the parameter as for later use, and thus above system (1a) and system (1b) have a modified version:where functions and are the densities of prey and predator at time , respectively. In terms of biology, all above positive constants have practical considerations. Parameters and denote the intrinsic growth rate of the prey and predator, respectively; and represent the carrying capacity of the environment; 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; is the biomass conversion and we denote as for convenience; is the intraspecific competition coefficient of the prey; is the Allee effect constant. The specific growth term governs the increase of prey in the lack of predator, while the specific growth term governs the increase of predator. The square term is an intrinsic decrease term on prey. The term 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–8] with Allee effect on prey. The coupled term 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 in denominator, which describes mutual interference among predators [11, 12]. This 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 and , system (2a) and system (2b) reduce to original system (1a) and system (1b).
When parameters (without Allee effect) and , 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 , system (2a) and system (2b) become model system (3a) and model system (3b) in [14], which was also independently and originally proposed in [9–11], while in [14], the authors discussed local and global asymptotic stability behavior of various equilibria and Hopf bifurcation occurs when parameter corresponding to reserved region crosses some critical values. To mimic the real-world scenario, they solved the inverse problem of estimation of system parameter 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). The rest of this paper is structured as follows. Preliminaries, such as boundedness, permanence, and existence of trivial equilibria, are given in Section 2. The 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.
2. Preliminaries
In this section, we devote to give some priori foundations. Before presenting the main results, we denote the first quadrant as , and its closure is . For biological consideration, system (2a) and system (2b) are defined on the domain and all the solutions are non-negative with initial conditions , i.e., 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.
2.1. Boundedness and Permanence
Theorem 1 (uniform boundedness). Suppose that a non-negative function and its partial derivatives and are continuous in , then the systemsubject to is uniformly bounded.
This Theorem 1 holds obviously after introducing an auxiliary function [1, 16]. The following theorem with the help of comparison principle in ODEs is about permanence of system (2a) and system (2b).
Theorem 2 (permanence). If parameters satisfywhere is a positive constant, , and , then system (2a) and system (2b) are permanent.
Proof. From above Theorem 1, we have a positive upper bound such thatThen, we obtain and sufficiently large , such thatBy using the lemmas in [17], we admit . That is to say, there exist sufficiently large , such that , . Then, we consider equation (2b) again. It yields ; here the function incorporating its Taylor expression isin which the second order derivative of isFrom the conditions in this theorem, we know that , when . Thus, , . By using the lemmas in [17] again, the proof is completed.
2.2. Existence of Trivial Equilibria
In this section, we will discuss the existence conditions of trivial equilibria of system (2a) and system (2b). Firstly, from [1], it is obvious that system (2a) and system (2b) have trivial equilibria: , , and , , where and
The point exists when , while the points all exist if and . If but , then the two equilibria and collide with each other and we denote this equilibrium as .
3. Existence of Interior Equilibrium
Here and below, we denote the interior equilibrium of system (2a) and system (2b) as or . This equilibrium must satisfy the following coupled algebraic equations:
From equation (10a), we have (1) or (if ) when ; (2) when . Bearing equation (10b) in mind, we have (1) (if exists) when ; (2) or (if ). The implicit derivative from equation (10b) isand we denote the positive root of a quadratic equation as for later use.
3.1. Cobweb Model
Based on above approximate analysis and the cobweb model, some cases about the existence of the interior equilibrium will be illustrated when isoclines from equations (10a) and (10b) all fall in . Case 1. If parameters satisfy then an interior equilibrium exists. Here, we take , , , , , , , , , , , , and ; the first equilibrium is and the second equilibrium is . Case 2. If parameters satisfy then an interior equilibrium exists. Here, we take , , , , , , , , , , , and ; the first interior equilibrium is , and the second equilibrium is . Case 3. If parameters satisfy then an interior equilibrium exists. Here, we take , , , , , , , , , , , , and ; an interior equilibrium is . Case 4. If parameters satisfy then an interior equilibrium exists. Here, we take , , , , , , , , , , , and ; an interior equilibrium is . Notice that another interior equilibrium verifies Case 1. Case 5. If parameters satisfy then an interior equilibrium exists (see the second equilibrium in Case 2). Case 6. If does not exist and parameters satisfy then an interior equilibrium exists. Here, we take , , , , , , , , , , , , and ; an interior equilibrium is . Case 7. If does not exist and parameters satisfy then an interior equilibrium exists. Here, we take , , , , , , , , , , , and ; an interior equilibrium is .
3.2. Polynomial Equations of and
From equation (10a), an expression of is
Substituting it into the equation (10b), a quintic algebraic equation of can be derived in the form of (see coefficients in Appendix A.1). Thanks 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 . For instance, ifwhere and , then there is a positive root such that . Case 1. If parameters satisfy then an interior equilibrium exists. Here, we take , , , , , , , , , , , , and ; an interior equilibrium is , and the following lemma is verified.
Lemma 1. Suppose that is a polynomial with real coefficients, , . If , then the equation has a positive root.
Proof. We only consider the special case , and thus . It is clear that the polynomial has a decomposition If we take a sufficiently large positive number such that , then , and thus we complete the proof.
Lemma 2. Suppose that is a real polynomial, , . If there is a positive number such that , then the equation has a positive root.
On the other hand, from equation (10b), we obtain an expression And a quintic algebraic equation of can be written in the form of (see coefficients in Appendix A.2). Suppose that the denominator and numerator of expression (23) are all positive for some ; then, . Notice that if and , there must exist a positive root for the quintic equation (see Lemma 1). Case 2. If parameters satisfy , , , , and , where is a root of the quadric equation in the denominator of (23), then an interior equilibrium exists. Here, we take , , , , , , , , , , , , and , and an interior equilibrium is . Case 3. If parameters satisfy , , , , and , then an interior equilibrium exists. Here, we take , , , , , , , , , , , , and , and an interior equilibrium is .
Remark 1. Here we rewrite equations (10a) and (10b) asFurthermore, if we sort them as , , the first eliminantalso yields above quintic equation . Similarly, if we sort them as , , the second eliminantyields the quintic equation .
4. 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.
4.1. Local Stability Analysis
The Jacobian matrix of system (2a) and system (2b) takes the following form , where four components are
For the trivial equilibrium with and the axial equilibrium with , we omit their stability [1]. When , the transformation and equations (2a) and (2b) yield (we still use symbol ), and thus is a saddle node and the parabolic sector is on the right half plane.
In the case that two axial equilibria exist, the Jacobian matrices are
Since , we have the following: (a) is a saddle point if ; (b) is an asymptotically stable node if . For the equilibrium , it is obvious that , and we have the following: (a) is an unstable node if ; (b) is a saddle point if ; (c) is an unstable higher-order singular point if . In the special case that two axial equilibria collide with each other, the Jacobian matrix isand thus is a higher-order singular point. If , is unstable.
For the interior equilibrium , its Jacobian matrix iswhere
Denoting a new discriminant with the trace and the determinant , we have(a)If and (a1) , then is an asymptotically stable node; (a2) , then is an asymptotically stable focus; (a3) , then is a saddle point.(b)If and (b1) , then is a center or a focus; (b2) , then is a saddle point.(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 point.
When but , is a stable (unstable) node if () (see Theorem 7.1 in Zhifen Zhang’s book [18] for more details). It is probable that has a cusp case of codimension at least 2 which ensures potential Bogdanov–Takens bifurcation when .
4.2. Global Asymptotic Stability
Combining the stability conclusions of the point in above section, the positive definite Lyapunov function ensures that is globally asymptotically stable if one of the following conditions holds:
Furthermore, conditions , and Theorem 1 deduce global asymptotical stability of pronto. If , , and , equilibria , , and do not exist, and is unstable, then Theorem 1 ensures that is globally asymptotically stable. For the further consideration of point , the following theorem explains its global asymptotic stability.
Theorem 3 (global asymptotic stability of ). If a unique interior equilibrium exists and parameters satisfythen is globally asymptotically stable.
Proof. Here we take an unbounded positive definite Lyapunov functionwith . Introducing new variables and , computing derivative along orbits of system (2a) and system (2b), we haveIt is obvious that is negative definite. Consequently, the Lyapunov function satisfies the asymptotic stability theorem in [19]. Thus, we complete the proof.
If conditions and hold, then , and is an asymptotically stable node or focus. Additionally, we could solve a potential interior equilibrium with two control variables and from such conditions, whereand parameters and are constrained bywhere coefficients areHere we set some values of parameters as , , , , , , , , , , and , and the unique globally asymptotically stable node with and is depicted in Figure 1 with characteristic direction which is solved from the characteristic functionafter we make the polar-coordinate-transformation , . The residual real simple root of equation in interval is . The power series of above characteristic function up to order two at the point admits indexes and with same negative sign, i.e., (odd number) and , where functionThus, actually shows a fact that trajectories enter the stable node along this direction.
4.3. Nonexistence of Limit Cycles
In this section, we consider nonexistence of closed orbits and limit cycles of system (2a) and system (2b). Firstly, taking a diffeomorphism which preserves the orientation of time, system (2a) and system (2b) are topologically equivalent to following system:since , where matrix
Notice that we still denote , , and as , , and . Above system is a -qualitatively equivalent polynomial form of system (2a) and system (2b), and it is more convenient to consider limit cycles [4, 5, 20, 21].
Theorem 4 (nonexistence of limit cycles). For system (42a) and system (42b), if parameters satisfythen there are no closed orbits and limit cycles in .
Proof. Here we take a Dulac function and calculate following partial derivative:where coefficientsand other unlisted coefficients are all nonpositive. Thus, we complete the proof.
Here we take over values of parameters from Section 4.2 but , and system (2a) and system (2b) do not have meaningful equilibria except a globally asymptotically stable node (the origin ) with the characteristic direction (the positive -axis) since the characteristic equation , and conditions in Theorem 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 and a globally asymptotically stable node when we use parameters from Section 4.2 but set , , , and . Conditions in Theorem 4 are also satisfied.
5. 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.
5.1. Saddle-Node Bifurcation
The two trivial equilibria collide with each other and system (47) has a unique boundary equilibrium when , if and , or
Then, there is a chance of bifurcation around this higher-order singular point. Here we choose as a bifurcation parameter and select eigenvector corresponding to the zero eigenvalue for matrix (30). The eigenvector corresponding to the zero eigenvalue for the transpose of matrix (30) is , where
Suppose , then the following transversality conditions hold:
Thus, we have following theorem by using Sotomayor’s theorem [22, 23].
Theorem 5 (saddle-node bifurcation). Suppose that the point exists; if , then system (47) undergoes a saddle-node bifurcation around point with respect to the bifurcation parameter .
5.2. Transcritical and Pitchfork Bifurcation
The equilibrium changes its stability when crosses the threshold ; in other words, is a higher-order equilibrium when . Let . Thus, we choose the parameter as a bifurcation parameter and an eigenvector corresponding to the zero eigenvalue for the Jacobian matrix when , where . The eigenvector corresponding to the zero eigenvalue for the transpose of matrix is ; then, the transversality conditions are
Thus, we have the following theorem by using Sotomayor’s theorem [22, 23].
Theorem 6 (transcritical bifurcation). Suppose that the two axial equilibria () coexist and ; if or , then system (47) undergoes a transcritical bifurcation around the point with respect to the bifurcation parameter .
For the special casethere is a chance the system (47) undergoes a pitchfork bifurcation. We still use the bifurcation parameter and eigenvectors and ; then, the fourth transversality condition iswhere , . Thus, we have another theorem by using Sotomayor’s theorem [22, 23].
Theorem 7 (pitchfork bifurcation). Suppose that the two axial equilibria () coexist and ; if condition (52) and hold, then system (47) undergoes a pitchfork bifurcation around the point with respect to the bifurcation parameter .
5.3. Hopf Bifurcation
In this section, we consider the Hopf bifurcation of system (2a) and system (2b). Here the point exists and we choose as bifurcation parameter. Suppose that are a pair of conjugate eigenvalues of matrix , where . The critical value satisfies
Then, system (47) undergoes a Hopf bifurcation around the point with respect to the bifurcation parameter .
We will calculate the first Lyapunov number at the point , which is used to determine the stability of limit cycles and Hopf bifurcation direction. The method and calculations of the first and second Lyapunov coefficients can be found in [22, 24, 25]. Therefore, translating the point to the origin by a linear transformation (I): , , system (47) in power series around the origin iswhere coefficients areand stands for some smooth functions. After that, by using a transformation (II): , with , the above system becomes a standard form:where
Let be two corresponding eigenvectors of a matrix such that , and ; the operation () 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 and arewhile is the two-dimensional Euclidean norm of . Define be the largest subspace invariant by the matrix and the generalized eigen subspace corresponding to the pair of purely imaginary eigenvalues , i.e., . That is to say, for any element , there must exist a linear expansion . Now we can construct a two-dimensional parameterized center manifoldin which and its complex conjugate is . Combining equations (57) and , we arrive at a complex equation without consideration:
The “chart” for the central manifold should be extracted from following differential equation:
So, substituting it into equation (63) and comparing coefficients of these , we recursively deriveand an equality from coefficient of the cubic term :
Here letter 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
Thus, the first Lyapunov number for the focus of planar system (57) is given by the formula . Since above expression is much too complicated, we need to present some numerical simulations and figures around the point with computer simulations.
Theorem 8 (nondegenerate Hopf bifurcation). Assume that the equilibrium exists and parameters satisfy condition (54) and ; then, system (47) undergoes a nondegenerate Hopf bifurcation around this equilibrium as parameter passes through the critical value . The Hopf bifurcation is supercritical (subcritical), the interior equilibrium is a multiple stable (unstable) focus with multiplicity one, and limit cycles are stable (unstable) if ().
Remark 2. The first Lyapunov coefficient can be defined by at times. On occasion, there are times when our system (55) exists some values of parameters such that 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]:which is an advanced version of or .
Finally, for the Hopf bifurcation, we numerically give following examples to simulate how the parameter controls dynamical behavior of system (2a) and system (2b).
Example 1. Consider Case 6 in Section 3.1 and Theorem 4 again. Figure 2(a) reveals the existence of this interior equilibrium in Case 6 by using the cobweb method. To investigate how the control parameter affects dynamical behavior of our system (2a) and system (2b), Figures 2(c) and 3 depict phase diagrams corresponding to values and , respectively. When , this interior equilibrium is an asymptotically stable focus since , , and . When , it gets , , and , and there exists a limit cycle (by using the Poincare–Bendixson theorem) around this unstable focus. Here we notice thatThis implies that the Hopf bifurcation occurs in system (2a) and system (2b) when . The first Lyapunov coefficient , and thus the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable. This interior equilibrium is a multiple stable focus with multiplicity one.
Besides, for a perturbed system with sufficiently small parameter vector in a neighbourhood of the origin in the parameter planeand Hopf bifurcation analysis with two bifurcation parameters and , we let and suppose that is an interior equilibrium of above perturbed system, where , is sufficiently small, andSubstituting it into and , the solutions and can be directly solved, which are written as the form up to third order:and hence the slope of approximation straight line of the Hopf bifurcation curve (see Figure 2(b)) in a small neighbourhood of the origin in parameter plane is approximately about . The supercritical Hopf bifurcation curve of system (71a) and system (71b) at the point is numerically defined by solution (73), i.e., , while the variables and both ensure the existence of andAs a matter of fact, we perceive the phenomenon that Theorem 8 describes the special case once lies on the horizontal line .
(a)
(b)
(c)
(a)
(b)
Example 2. Here we set parameters as well as Case 3 in Section 3.2. Figures 4 and 5 depict phase diagrams corresponding to values and , respectively. One should notice thatSimilar to above example, system (2a) and system (2b) undergo a Hopf bifurcation when passes through . The first Lyapunov coefficient is also found to be negative. Thus, the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable. The interior equilibrium is also a multiple stable focus with multiplicity one. Figure 6 is the Hopf bifurcation curve corresponding to system (69) with slope and
(a)
(b)
6. 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]. The polynomial’s method, derived from eliminants and , can be extended to more complicated polynomial systems. Some conclusions are the same as reference [1], such as the uniform boundedness (Theorem 1), the existence of equilibria , the nonexistence (Theorem 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 and , need to be researched further with the help of topologically equivalent systems or the “blow-up” method (horizontal and vertical blow-ups). Other parameters can also be considered as a bifurcation parameter in Hopf bifurcation, although it is described by Hopf bifurcation curve , for instance . Notice that Theorem 2 (permanence) can be extended to its reaction-diffusive version [26, 27]. Under the conditions or the discussion of Theorem 3, the Turing instability of its corresponding reaction-diffusion system subject to the homogeneous Neumann boundary conditions [28]:will not occur, and thus it is available to find potential Hopf bifurcation points and consider transversality conditions if we choose as Hopf bifurcation parameter. Here and are two positive diffusive constants; is the Laplacian operator; and are the densities of prey and predator, respectively; is an one-dimensional bounded domain with smooth boundary ; the symbol is the outer flux, and no flux boundary condition is imposed; thus, the system is closed [29]; the admissible initial functions and are all continuous functions on ; to describe an environment surrounded by dispersal barriers, we take zero flux at [28]. The method of calculating the first Lyapunov coefficient in Section 5.3 is a reference for deducing of Hopf bifurcation direction in reaction-diffusion systems [29, 30].
Though the dynamical behavior of predator-prey systems in single species or multispecies has been researched by many previous literatures, we still need further study in biomathematics, especially the phytoplankton and zooplankton systems. Meanwhile, how to keep ecosystems in balance or coexistence states and avoid harmful effect is our next direction in this area.
Appendix
A. Polynomials and
The polynomial in Section 3.2 is , where coefficients are
The polynomial in Section 3.2 is , where coefficients are
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors thank teachers Min Zhao and Chuanjun Dai. This study was supported by the National Natural Science Foundation of China (grant nos. 31570364 and 61871293) and the National Key Research and Development Program of China (grant no. 2018YFE0103700).