Stability and bifurcation analysis of two-species competitive model with Michaelis–Menten type harvesting in the first species

In this paper, a two-species competitive model with Michaelis–Menten type harvesting in the first species is studied. We have made a detailed mathematical analysis of the model to describe some important results that may be produced by the interaction of biological resources. The permanence, stability, and bifurcation (saddle-node bifurcation and transcritical bifurcation) of the model are investigated. The results show that with the change of parameters, two species could eventually coexist, become extinct or one species will be driven to extinction and the other species will coexist. Moreover, by constructing the Lyapunov function, sufficient conditions to ensure the global asymptotic stability of the positive equilibrium are given. Our study shows that compared with linear harvesting, nonlinear harvesting can exhibit more complex dynamic behavior. Numerical simulations are presented to illustrate the theoretical results.


Introduction
During the past decade, a competitive system has been extensively investigated by many scholars , many excellent results concerned with the permanence, extinction, and global attractivity of the competition system have been obtained.
Traditional two-species Lotka-Volterra competitive system is as follows: where x 1 and x 2 denote the population density of the first and second species at time t, respectively. b i , a ij , i, j = 1, 2, are positive constants. System (1.1) has been investigated in mathematical biology books [25]. Depending on the relationship of the coefficients of system (1.1), it has the following dynamic behaviors.
(1) If holds, system (1.1) has a unique positive equilibrium, which is globally attractive.
holds, then the second species will be driven to extinction, while the first species will approach b 1 a 11 . On the other hand, the study of resource-management, including fisheries, forestry, and wildlife management, has great importance as human activity is the main cause of the extinction of endangered species. It is necessary to harvest the population, but harvesting should be regulated so that both the ecological sustainability and conservation of the species can be implemented in a long run. Many scholars are interested in establishing appropriate biological models to further understand the scientific management of renewable resources. For example, based on model (1.1), Sharma and Samanta [26] further considered the harvesting and obtained the following model: dy(t; p) dt = y(t) (r 2l ) (1-p) (r 2u ) p -(b 22l ) (1-p) (b 22u ) p y(t) where p ∈ [0, 1]. Due to the difficulty in the estimation of the model parameters, the authors argued that taking account of the imprecise of biological parameter values makes some situations more realistic. They developed a method to handle imprecise parameters. In addition, they discussed the existence and stability of the system equilibria, as well as the bionomic equilibrium and optimal harvesting policy. The equilibrium solution of the control problem was given and the harvest strategy was optimized dynamically. Ecosystem with harvesting has been extensively investigated by many scholars [27][28][29][30][31][32][33][34][35][36][37]. May et al. [31] proposed two types of harvesting: (1) constant harvest and (2) linear harvest. For the first case, it is impossible to harvest a certain number of species every year. Although this type of harvest may be relatively easy to study, it is not a reality. For the second case, the harvesting term takes the form h(x) = qEx, obviously, when x or E tends to infinity, h(x) tends to infinity. Clearly this contradicts the facts, because in reality there is limited harvesting capacity or number of species, so the amount of species that can be harvested is limited. To overcome the drawback of the two kinds of harvesting, Clark [32] proposed the Michaelis-Menten type harvesting h(x) = qEx mE+nx when x or E tends to infinity, h(x) tends to qE n or qx m . In this type of harvesting, if the number of species tends to infinity, the final harvest depends on the harvesting capacity, and if the harvesting capacity tends to infinity, the final harvest depends on the number of species, which is in accordance with the actual situation. Since such kind of harvesting is more suitable, it brings many scholars to do works in this direction (see [33][34][35][36][37] and the references cited therein). For example, Yu, Chen, and Lai [35] introduced Michaelis-Menten type harvesting into a May type cooperative system and discussed the extinction of the first species and the global attraction of the unique positive equilibrium. Chen [36] studied the Lotka-Volterra commensal symbiosis model with Michaelis-Menten type harvesting; the modified model takes the following: where x and y denote the population density of the two species at time t, respectively.
q denotes the fishing coefficient of the first species and E denotes the fishing effort. r 1 , The results show that the system has a globally asymptotically stable positive equilibrium. In addition, the two species can coexist stably when α and K 2 are large enough. On the contrary, the first species will be driven to extinction.
After that, Liu et al. [37] considered the two-species amensalism model with Michaelis-Menten type harvesting and a cover for the first species; the modified model takes the following: (1.4) Here, x and y denote the population density of the two species at time t, and k is the refuge (0 < k < 1). When the parameters satisfy certain conditions, saddle-node bifurcation and transcritical bifurcation will occur in the system, and the maximum threshold of species where r 1 , r 2 , k 1 , k 2 , α 1 , α 2 , q 1 , m 1 , h 1 , and E are all positive. For simplicity, we make the following nondimensional scheme: Dropping the bars, we have the following system: , and the initial conditions Biologically, we consider system (1.6) is defined on the set The organization of this paper is as follows. The basic properties of the model are discussed in the next section. We analyze the existence of the equilibria of the system in Sect. 3. The local stability of the equilibria of the system are investigated in Sect. 4. We consider the global stability of the positive equilibrium in Sect. 5. The possible bifurcation of the system is studied in Sect. 6. Numerical simulation is presented to show the feasibility of theoretical results in Sect. 7, and a brief discussion of our results is given in the last section. Proof Since

Basic properties of the model
So x(t) > 0, y(t) > 0 for all t ≥ 0 with initial condition (1.7). This completes the proof.
Proof Based on the positivity of variable x, y, from system (1.6), we have This completes the proof.
Proof From (2.2) and (2.4), for ε > 0 small enough, there is T > 0 such that, for t > T, we have Then from the first equation of system (1.6), one could get From Lemma 2.1, when 1 - Let ε → 0 in this inequality, we can obtain In the same way, we can obtain the following results for y(t): where ω 2 = 1a 2 > 0, i.e., 0 < a 2 < 1. So we choose m = min(ω 1 , ω 2 ), M = 1. This completes the proof.

Stability of equilibria
Theorem 4.1 For the boundary equilibria E 0 and E 1 of system (1.6), which always exist, we have: (1) E 0 is always unstable.
(2) The Jacobian matrix of system (1.6) at E 1 is given by It is obvious that the eigenvalues of J(E 1 ) are λ 1 = 1-a 1 -b 1 c 1 and λ 2 = -ρ < 0, so the stability of J(E 1 ) depends on λ 1 , we cannot directly come to the conclusion.
Let us first move E 1 to the origin by transforming (X, Y ) = (x, y -1) and make Taylor's expansion around the origin, under which system (1.6) is as follows: Next, making the following transformation and letting τ = ρt, system (4.4) becomes where And From the implicit function theorem, there exists a function Y 1 = ϕ(X 1 ) that satisfies ϕ(0) = 0 and P(X 1 , ϕ(X 1 )) = 0. By using the second equation of system (4.6), we can obtain Substituting (4.7) into the first equation of system (4.6),we have (1) E 21 is always unstable.
(2) For E 22 , we have the following results: If a 2 (c 1 -1) = -2, E 23 is a saddle node; on the contrary, E 23 consists of a hyperbolic sector and an elliptic sector.
Case 1: If λ 2 = 0, the stability of E 23 can be proved by the same method in Theorem 4.2(2), so E 23 is a saddle node if x 23 = 1 a 2 .
This completes the proof.

Theorem 4.3 When the positive equilibria exist, we have:
(1) E 31 is a stable node.
Proof Notice that when the positive equilibria exist, (4.1) can be simplified as follows: where i = 1, 2, 3. Thus, (1) According to (3.5), we get that   Figure 2 is an enlarged version of Fig. 1 so the eigenvalues of J(E 33 ) are λ 1 < 0 and λ 2 = 0, we use the same method as Theorem 4.2(2) and easily get E 33 is a saddle node.
Proof Through the discussion in Theorem 4.3, we can get a 1 a 2 c 1 ) ρy * .
In order to verify the above results, let a 1 = 0.6, a 2 = 2, b 1 = 0.42, c 1 = 0.4, ρ = 1. By simple computation, we have A > 0, B < 0, C < 0. E * (x * , y * ) and E 21 (0.03542, 0) are saddle, E 1 (0, 1) and E 22 (0.56458, 0) are stable nodes. From Fig. 3, it is easy to get the red line that divides the first quadrant into two parts, recorded as I (left one) and II (right one). Assume that the initial conditions are located in region I, all the solutions tend to E 1 (0, 1) which is a stable manifold, E * (x * , y * ) and E 21 (0.03542, 0) are unstable manifolds. From the biological point of view, when the initial values are located in region I, the first species will be driven to extinction. On the contrary, assume that the initial conditions are located in region II, all the solutions tend to E 22 (0.56458, 0) which is a stable manifold, E * (x * , y * ) and E 21 (0.03542, 0) are unstable manifolds. From the biological point of view, when the initial values are located in region II, the second species will be driven to extinction. This is the bistable phenomenon which is shown in Fig. 3.
This completes the proof.

Global stability of equilibrium
Theorem 5.1 When E * is locally stable, which is globally asymptotically stable.
Proof We will adapt the idea of Yu [40] to prove Theorem 5.1. More precisely, we will consider a Lyapunov function It is easy to see that the function V (x, y) is zero at the equilibrium E * (x * , y * ), which is positive everywhere in the first quadrant except at E * . Then the time derivative of V (x, y) along the trajectories of (1.6) is Note that B > 0, C > 0 when E * exists, and we have b 1 c 1 < 1a 1 < c 1 (1a 1 a 2 ), i.e., b 1 c 1c 1 < 0, so one could obtain D + V (t) < 0 strictly for all x, y > 0 except the equilibrium E * (x * , y * ), where D + V (t) = 0. Hence, V (x, y) satisfies Lyapunov's asymptotic stability theorem, then the equilibrium E * (x * , y * ) of system (1.6) is globally asymptotically stable. From the biological point of view, two species can always coexist (see Fig. 4). This completes the proof.

Bifurcation analysis
In this section, we are interested in studying the various possible bifurcations of system (1.6). From Theorem 3.1, we know that system (1.6) may undergo saddle-node bifurcation at E 23 and E 23 , respectively, and transcritical bifurcation around the equilibria E 0 and E 1 , which is a very interesting phenomenon.

Saddle-node bifurcation
From Theorem 3.1, it is easy to find that when c 1 < b 1 < b * 1 and 0 < c 1 < 1 the system has two different boundary equilibria, which may coincide if b 1 = b * 1 and 0 < c 1 < 1 or which The emergence or appearance of the equilibria is due to the saddle-node bifurcation at E 23 (see Fig. 5 and Fig. 6).  Proof By the proof of Theorem 4.2, we have an eigenvalue of J(E 23 ) that is zero, named λ 1 . Let V 1 and W 1 represent the eigenvectors of λ 1 for the matrices J(E 23 ) and J(E 23 ) T , respectively. After simple calculation, we have  Clearly, we can get that V 1 and W 1 satisfy This implies that when b 1 = b 1SN , the saddle-node bifurcation occurs at E 23 . This completes the proof.
Similarly, the conditions for the existence of positive equilibria of system (1.6) are given in Theorem 3.1, and we could find that when 2 > 0, A < 0, B < 0, C < 0, the system has two different positive equilibria, which may coincide if 2 = 0 and A < 0, B < 0 or which may disappear if 2 < 0. The emergence or appearance of the equilibria is due to the saddlenode bifurcation at E 33 (see Fig. 5 (a) and Fig. 7). Theorem 6.2 System (1.6) undergoes a saddle-node bifurcation around E 33 when b 1 =b 1SN and A < 0, B < 0, where b 1 is the bifurcation parameter.
Proof By the proof of Theorem 4.3, we know that an eigenvalue of J(E 23 ) is zero, named λ 1 . Let V 2 and W 2 represent the eigenvectors of λ 1 for the matrices J(E 33 ) and J(E 33 ) T , respectively. After simple calculation, we have Moreover, So, we can obtain V 2 and W 2 satisfy which means that when b 1 = b 1SN , the saddle-node bifurcation occurs at E 33 . This completes the proof.

Transcritical bifurcation
Through the discussion of Theorem 3.1, we find an interesting phenomenon: when b 1 = c 1 , E 21 will coincide with E 0 if 0 < c 1 < 1; E 22 will coincide with E 0 if c 1 > 1. Therefore, the emergence of this phenomenon is owing to the transcritical bifurcation at E 0 (see Fig. 8).
Then we obtain the following. Proof Here, we use Sotomayor's theorem to verify the transversality conditions for transcritical bifurcation. The Jacobian matrix of system (1.6) evaluated at the point E 0 is given by (3.2). Obviously, the eigenvalue λ 1 = 0 of J(E 0 ) if b 1 = c 1 . Let V 3 and W 3 be the eigenvectors of J(E 0 ) and J(E 0 ) T corresponding to λ 1 = 0, respectively. Then we can obtain Furthermore, Thus, we have So from Sotomayor's theorem system (1.6) undergoes a transcritical bifurcation at E 0 . This completes the proof.
In the same way, for the positive equilibria of system (1.6), when b 1 = c 1 (1a 1 ) and A < 0, E 32 will coincide with E 1 if B > 0; E 31 will coincide with E 1 if B < 0. Hence, the appearance of this phenomenon is owing to the transcritical bifurcation at E 1 (see Fig. 9). Then we can get the following.  Proof Here, we use Sotomayor's theorem to verify the transversality conditions for transcritical bifurcation. The Jacobian matrix of system (1.6) evaluated at the point E 1 is given by (3.2). Clearly, the eigenvalue λ 1 = 0 of J(E 1 ) if b 1 = c 1 (1a 1 ). Let V 4 and W 4 be the eigenvectors of J(E 1 ) and J(E 1 ) T corresponding to λ 1 = 0, respectively. Then we can Furthermore, Thus, we have So from Sotomayor's theorem system (1.6) undergoes a transcritical bifurcation at E 1 . This completes the proof.

Numeric simulations
Now, let us give the following examples to illustrate the main results.
(3) For a 2 = 4, we get E 23 consists of a hyperbolic sector and an elliptic sector. Figure 5 shows the above results.  Figure 6 shows the above results.

Conclusion
In this paper, we have considered a two-species competitive system with Michaelis-Menten type harvesting. The model shows rich dynamic behaviors. We have studied the permanence condition of the system, and by analyzing the stability of the system equilibria, we obtained that the system cannot collapse for any parameter value as the origin is always unstable. In addition, from the global asymptotic stability of the positive equilibrium, it can be seen that two species can coexist stably under appropriate conditions. We also get that there are two different cases of bistability in the system: on the one hand, a boundary equilibrium and a positive equilibrium are globally asymptotically stable; on the other hand, the two boundary equilibria are globally asymptotically stable. Qualitative analysis indicates that Michaelis-Menten type harvesting plays an important role in the dynamic behaviors and bifurcations of the system. Firstly, the parameters b 1 and c 1 of the Michaelis-Menten type harvesting term will affect the number and stability of the system equilibria, compared with system (1.1), the boundary equilibria and positive equilibria of system (1.6) are both increased. Secondly, from the Michaelis-Menten type harvesting term h(x) = b 1 x c 1 +x , we know that system (1.6) is supplementary to system (1.1). Because system (1.6) becomes unharvested situation if b 1 = 0, system (1.6) becomes constant harvested situation if c 1 = 0. Thirdly, we give a strict proof of the bifurcation of system (1.6) by Sotomayor's theorem, which has important ecological significance. Through saddle-node bifurcation and transcritical bifurcation, one could obtain the maximum threshold without extinction risk of species in continuous harvest. This provides important reference for decision makers to make reasonable strategies to ensure the sustainable development of ecosystem and maximize economic benefits.