Financing of the global shift to renewable energy and energy efficiency

1.   Introduction

The interaction between two or more species has been an important and interesting issue in biology and ecology since the famous Lotka-Volterra model ^[1,2] was proposed. The interaction between different species will generate rich interesting dynamics of biological species, and exhibit the complexity and diversity ^[3,4]. Amensalism, as a typical type of interaction between the species, has been intensively considered in the last decades. Amensalism describes a basic biological interaction in nature, where one species inflicts harm on another not affected by the former, which means that it does not receive any costs or benefits to itself. The first pioneer work for the investigations of amensalism model is due to Sun ^[5], who, in 2003, proposed the following two-species amensalism model:

where $x = x(t)$ and $y = y(t)$ represent the population densities of two species at time $t$, respectively; $r_1$ and $k_1$ describe the intrinsic growth rate and the carrying capacity of the first species, respectively; $r_2$ and $k_2$ describe the intrinsic growth rate and the carrying capacity of the second species, respectively; $\alpha$ describes the effect of the second species on the first species. They explored the stability properties of all possible equilibria in system (1.1).

Setting $a_{11} = \frac{r_1}{k_1}$, $a_{12} = \frac{r_1\alpha}{k_1}$ and $a_{22} = \frac{r_2}{k_2}$, we can see that model (1.1) can be rewritten as

where all parameters $r_1$, $r_2$, $a_{11}$, $a_{12}$, and $a_{22}$ are positive real numbers.

Ever since the first amensalism model was presented, the complicated dynamics of amensalism models have been studied extensively (see ^{[6,7,8,9,10,11,12,13,14,15]} and the references cited therein). For example, in ^[6], Guan and Chen considered a two-species amensalism model with Beddington-DeAngelis functional response and gave some comprehensive bifurcation and global dynamics of this system. Recently, Luo et al. ^[12] proposed an amensalism model with Holling-Ⅱ functional response and weak Allee effect for the first species, and discussed its local dynamics and global structure.

In the real world, from the point of view of human needs, the management of renewable resources, the exploitation of biological resources, and the harvesting of populations are commonly practiced in forestry, fishery, and wildlife management ^[16,17,18]. Hence, it is necessary to introduce and to consider the harvesting of species in population models. During the last decade, population models with harvesting and the role of harvesting in the management of renewable resources have received significant attention from the researchers ^{[19,20,21,22,23,24]}. Generally speaking, there are three types of harvesting: 1) constant harvesting ^[25]; 2) linear harvesting ^[25]; 3) non-linear harvesting ^[26,27]. As is well known, non-linear harvesting is more realistic from biological and economical points of view ^[28]. In ^[26], Clark first proposed the non-linear harvesting term $h(E, x) = \frac{qEx}{cE+lx}$, which is called Michaelis-Menten type functional form of catch rate, here $q$ is the catchability coefficient, $E$ is the external effort devoted to harvesting, $c$ and $l$ are constants. After that, many scholars started to focus attention on the influence of the Michaelis-Menten type harvesting on the population systems ^{[28,29,30,31,32,33,34]}. For example, a reaction-diffusion predator-prey model with non-local delay and Michaelis-Menten type prey-harvesting was investigated by Zhang et al. ^[32], they obtained that the discrete and non-local delays are responsible for a stability switch in this system, and a Hopf bifurcation occurs as the delays pass through a critical value. In ^[33], Hu and Cao discussed a predator-prey system with the nonlinear Michaelis-Menten type predator harvesting and gave a detailed analysis of the stability and bifurcation for this system.

These studies have shown that harvesting activity has important influence to dynamics of population systems and the harvested models can exhibit richer dynamics compared to the model with no harvesting. However, seldom did scholars consider the dynamical behaviors of an amensalism model with a harvesting term. Accordingly, inspired by the previous works, based on the system (1.2), we now introduce the Michaelis-Menten type harvesting to the second species, then we propose the following amensalism model:

where the parameters $r_1$, $r_2$, $a_{11},$ $a_{12}$, and $a_{22}$ are positive constants, and this term $\frac{qEy}{cE+ly}$ represents Michaelis-Menten harvesting.

Setting $\overline{t} = r_1t, \ \overline{x} = \frac{a_{11}}{r_1}x, \ \overline{y} = \frac{a_{12}}{r_1}y, \ \alpha = \frac{a_{12}qE}{lr_1^2}, \ \delta = \frac{r_2}{r_1}, \ \beta = \frac{a_{22}}{a_{12}}, \ \gamma = \frac{a_{12}cE}{lr_1},$ and dropping the bars, then system (1.3) is transformed into

The layout of this paper is arranged as follows. In the next section, the global dynamics including positivity and boundedness of solutions, number of equilibria, local asymptotical stability, codimension one bifurcations, dynamical behaviors near infinity, closed orbits analysis and global phase portraits are shown for system (1.4). Further, in Section 3, we also obtain the global dynamics of system (1.1) without harvesting term. Numerical simulations and discussions are displayed by Section 4 for demonstrating the theoretical results and the impact of harvesting term. Finally, a brief conclusion is presented in Section 5.

2.   Global dynamics of system (1.4)

In this section, we mainly investigate the global dynamics of system (1.4), which include the positivity and boundedness of solutions, the existence and local stability analysis of equilibria, all possible bifurcation behaviors and the global structure of system (1.4).

2.1.   Positivity and boundedness of solutions

Here, we give the positivity and boundedness of the solutions of model (1.4) in the region $\mathbb{R}_{+}^2 = \left\lbrace (x, y): x\geq 0, y\geq 0\right\rbrace$.

2.1.1.   Positivity

Lemma 2.1. All solutions of model $(1.4)$ with positive initial value are positive for all $t\geq0$.

Proof. Solving model (1.4) with positive initial condition $(x(0), y(0))$ gives the result:

Therefore, we can easily see that each solution of model (1.4) starting from the positive initial condition $(x(0), y(0))$ still remains in the first quadrant.

2.1.2.   Boundedness

Lemma 2.2. All solutions of system $(1.4)$ with positive initial value are bounded in region $\varOmega = \left\lbrace(x(t), y(t)): 0\leq x(t) < 1 \ \text{and}\ 0 \leq y(t) < \frac{\delta}{\beta}\right\rbrace$ for the large enough time $t$.

Proof. If the initial value is selected in $\mathbb{R}_{+}^2$, each solution of model (1.4) remains positive from the above result. One can see

from the first equation of model (1.4). And a standard comparison theorem shows that

Also, the following inequality

is obtained from the second equation of system (1.4). It yields

Hence, there exists a sufficiently large time $T$, such that $0\leq x(t) < 1$ and $0\leq y(t) < \frac{\delta}{\beta}$ for $t > T$. In summary, all the solutions of system (1.4) are always bounded. Thus we complete the proof.

2.2.   Existence and stability of equilibria

Theorem 2.1. System $(1.4)$ always has two boundary equilibria $E_0(0, 0)$ and $E_1(1, 0)$ for all positive parameters. For the existence of other equilibria, we have

$1)$ For the possible positive equilibria:

$i)$ if $\alpha > \alpha^*$, then system $(1.4)$ has no positive equilibria;

$ii)$ if $\alpha = \alpha^*$ and $0 < \delta-\beta\gamma < 2\beta$, then there is a unique positive equilibrium $E_{33}(1-y_{3}^*, y_{3}^*)$;

$iii)$ if $\alpha^{**} < \alpha < \alpha^*$, when $\alpha < \delta+\delta\gamma-\beta-\beta\gamma$ and $\delta-\beta\gamma > 0$, or $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, there is one positive equilibrium $E_{31}(1-y_{1}^*, y_{1}^*)$; when $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, there are two distinct positive equilibria $E_{31}(1-y_{1}^*, y_{1}^*)$ and $E_{32}(1-y_{2}^*, y_{2}^*)$;

$iv)$ if $\alpha = \alpha^{**}$ and $0 < \delta-\beta\gamma < \beta$, then $E_{31}$ coincides with $E_1$ and there is one positive equilibrium $E_{32}(1-y_{2}^*, y_{2}^*)$, where $y_{2}^* = \frac{\delta-\beta\gamma}{\beta}$;

$v)$ if $0 < \alpha < \alpha^{**}$ and $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$, then there is one positive equilibrium $E_{32}(1-y_{2}^*, y_{2}^*)$.

$2)$ For the other possible boundary equilibria:

$i)$ if $\alpha > \alpha^*$, then system $(1.4)$ has no other boundary equilibria; there is one boundary equilibrium $E_{22}(0, y_{2}^*)$.

$ii)$ if $\alpha = \alpha^*$ and $\delta-\beta\gamma > 0$, then there is a boundary equilibrium $E_{23}(0, y_{3}^*)$;

$iii)$ if $\alpha^{**} < \alpha < \alpha^*$ and $\delta-\beta\gamma > 0$, then there are two distinct boundary equilibria $E_{21}(0, y_{1}^*)$ and $E_{22}(0, y_{2}^*)$;

$iv)$ if $\alpha = \alpha^{**}$ and $\delta-\beta\gamma > 0$, then $E_{21}$ coincides with $E_0$ and there is one boundary equilibrium $E_{22}(0, y_{2}^*)$, where $y_{2}^* = \frac{\delta-\beta\gamma}{\beta}$;

$v)$ if $0 < \alpha < \alpha^{**}$, then

Where $\alpha^* = \frac{(\delta+\beta\gamma)^2}{4\beta}$, $\alpha^{**} = \delta\gamma$, $y_1^* = \frac{\delta-\beta\gamma- \sqrt{(\delta+\beta\gamma)^2-4\alpha\beta}}{2\beta}$, $y_2^* = \frac{\delta-\beta\gamma+\sqrt{(\delta+\beta\gamma)^2-4\alpha\beta}}{2\beta},$ and $y_3^* = \frac{\delta-\beta\gamma}{2\beta}.$

Proof. It is easy to see that system (1.4) always possesses two boundary equilibria given by $E_0(0, 0)$ and $E_1(1, 0)$ for all positive parameters.

1) If $x\neq 0$ and $y\neq 0$, system (1.4) has a positive equilibrium which satisfies

For the positive equilibria, $y$ must satisfy $0 < y < 1$. Let $\Delta$ denote the discriminant of the second equation of (2.1), namely,

and let

Obviously, when $0 < \alpha < \alpha^*$, the second equation of (2.1) has two roots

When $\alpha = \alpha^*$, the second equation of (2.1) has one root

Let

and denote by

Then we have

H1) if $\alpha > \alpha^{**}$, then $k(0) > 0$;

H2) if $\alpha = \alpha^{**}$, then $k(0) = 0$;

H3) if $\alpha < \alpha^{**}$, then $k(0) < 0$.

In addition, $0 < \alpha^{**}\leq \alpha^{*}$, and $\alpha^{**} = \alpha^{*}$ if and only if $\delta = \beta\gamma$.

Based on the above analysis, combining with H1), H2) and H3), we can derive that:

ⅰ) If $\alpha > \alpha^*$, then $\Delta < 0$ which implies $k(y) = 0$ has no real roots. So system (1.4) has no positive equilibria.

ⅱ) If $\alpha = \alpha^*$, then $\Delta = 0$, so $k(y) = 0$ has a unique real roots $y_3^* = \frac{\delta-\beta\gamma}{2\beta}$, combing with the condition $0 < y_3^* < 1$, we get $0 < \delta-\beta\gamma < 2\beta$. In this situation, there exists a unique positive equilibrium $E_{33}(1-y_{3}^*, y_{3}^*)$.

ⅲ) If $\alpha^{**} < \alpha < \alpha^*$, then $\Delta > 0$ and $k(0) > 0$, implying $k(y) = 0$ two positive real roots $y_1^*$ and $y_2^*$ if $\delta-\beta\gamma > 0$. Moreover, when $k(1) < 0$, namely, $\alpha < \delta+\delta\gamma-\beta-\beta\gamma$, we have $0 < y_1^* < 1 < y_2^*$, that is, there is one positive point $E_{31}(1-y_{1}^*, y_{1}^*)$; when $k(1) = 0$ and the symmetry axis of $k(y)$ is $y = \frac{\delta-\beta\gamma}{2\beta} < 1$, namely, $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, we have $0 < y_1^* < y_2^* = 1$, which means there is one positive point $E_{31}(1-y_{1}^*, y_{1}^*)$; when $k(1) > 0$ and the symmetry axis of $k(y)$ is $y = \frac{\delta-\beta\gamma}{2\beta} < 1$, namely, $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, we have $0 < y_1^* < y_2^* < 1$, that is there are two different positive equilibria $E_{31}(1-y_{1}^*, y_{1}^*)$ and $E_{32}(1-y_{2}^*, y_{2}^*)$.

ⅳ) If $\alpha = \alpha^{**}$, then $\Delta > 0$ and $k(0) = 0$, hence $k(y) = 0$ has two real roots $y_1^* = 0$ and $0 < y_2^* = \frac{\delta-\beta\gamma}{\beta} < 1$ if $0 < \delta-\beta\gamma < \beta$, thus $E_{31}$ coincides with $E_{1}$ and there is one positive equilibrium $E_{32}(1-y_{2}^*, y_{2}^*)$.

ⅴ) If $\alpha < \alpha^{**}$, then $\Delta > 0$ and $k(0) < 0$, so there are $y_1^* < 0$ and $y_2^* > 0$. When $k(1) > 0$, i.e., $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$, we get one positive root $0 < y_2^* < 1$, thus there exits one positive equilibrium $E_{32}(1-y_{2}^*, y_{2}^*)$.

The proof of the corresponding boundary equilibria in 2) is similar to that of 1), and hence omitted here. This completes the proof of Theorem 2.1.

In the following, we focus on the local stability of each equilibrium of system (1.4).

Theorem 2.2. For the boundary equilibria $E_0$ and $E_1$, the following statements are true.

$1)$ If $0 < \alpha < \alpha^{**}$, then $E_0$ is a hyperbolic unstable node and $E_1$ is a hyperbolic saddle.

$2)$ If $\alpha > \alpha^{**}$, then $E_0$ is a hyperbolic saddle and $E_1$ is a hyperbolic stable node.

$3)$ If $\alpha = \alpha^{**}$, then

$i)$ when $\delta\neq \beta\gamma$, $E_0$ and $E_1$ are both saddle-nodes;

$ii)$ when $\delta = \beta\gamma$, $E_0$ is a non-hyperbolic saddle and $E_1$ is a non-hyperbolic stable node.

Proof. The Jacobian matrix of system $(1.4)$ evaluated at any equilibrium is

where $H(x, y) = \delta-2\beta y-\frac{\alpha\gamma}{(\gamma+y)^2}$.

For the equilibrium $E_0$, the Jacobian matrix is

and the two eigenvalues of $J(E_0)$ are $\lambda_1(E_0) = 1 > 0$ and $\lambda_2(E_0) = \delta-\frac{\alpha}{\gamma}$. Obviously, if $0 < \alpha < \alpha^{**}$, then $\lambda_2(E_0) > 0$, so $E_0$ is a hyperbolic unstable node. If $\alpha > \alpha^{**}$, then $\lambda_2(E_0) < 0$, thus $E_0$ is a hyperbolic saddle. If $\alpha = \alpha^{**}$, then $\lambda_2(E_0) = 0$, this means that the equilibrium $E_0$ is non-hyperbolic, so it is hard to directly judge its type from their eigenvalues. We further discuss its stability properties by applying Theorem 7.1 in Chapter 2 in ^[35].

In order to change system (1.4) into a standard form, we expand system (1.4) in power series up to the fourth order around the origin

where $Q_0(y)$ represents a power series with the terms $y^i$ $(i\geq 5)$. From $\frac{dx}{dt} = 0,$ we obtain that there is a unique implicit function $x = \varphi_0(y) = 0$ such that $\varphi_0(y)+P(\varphi_0(y), y) = 0$ and $\varphi_0(0) = \varphi_0'(0) = 0$. Then substituting $x = \varphi_0(y) = 0$ into the second equation of (2.7), we get that

Because $\alpha = \alpha^{**} = \delta\gamma$, the coefficient at $y^2$ is $\frac{\delta}{\gamma}-\beta$. By Theorem 7.1 in ^[35], we have when $\delta\neq \beta\gamma$, the equilibrium $E_0$ is a saddle-node.

When $\delta = \beta\gamma$, (2.8) becomes

Employing the notations of Theorem 7.1 in Chapter 2 in ^[35], we have $m = 3$ and $a_m = -\frac{\beta}{\gamma} < 0$, so the equilibrium $E_0$ is a non-hyperbolic saddle.

For the equilibrium $E_1$, the Jacobian matrix is

The two eigenvalues of the above matrix $J(E_1)$ are $\lambda_1(E_1) = -1 < 0$ and $\lambda_2(E_1) = \delta-\frac{\alpha}{\gamma}$. Clearly, if $\alpha < \alpha^{**}$, then $\lambda_2(E_1) > 0$, so $E_1$ is a hyperbolic saddle. If $\alpha > \alpha^{**}$, then $\lambda_2(E_1) < 0$, then $E_1$ is a hyperbolic stable node. If $\alpha = \alpha^{**}$, then $\lambda_2(E_1) = 0$, the equilibrium $E_1$ is non-hyperbolic and its stability cannot be given directly from the eigenvalues. In order to determinate the stability of $E_1$, we translate $E_1$ to the origin by the translation $(\overline{x}, \overline{y}) = (x-1, y)$ and expand system (1.4) in power series up to the fourth order around the origin, which makes system (1.4) to be the following form:

where $Q_1(\overline{y})$ represents a power series with terms $\overline{y}^i$ $(i\geq 5)$.

To transform the Jacobian matrix into a standard form, we use the invertible translation

then system (2.11) becomes

By introducing a new time variable $\tau = -t$, we get

Based on the implicit function theorem, from $\frac{du}{d\tau} = 0$, we can deduce a unique function

which satisfies $\varphi_1(0) = \varphi_1'(0) = 0$ and $\varphi_1(v) + P(\varphi_1(v), v) = 0$. Then substituting it into the second equation of (2.14), we have

From $\alpha = \delta\gamma$ it follows that the coefficient at $v^2$ is $\beta-\frac{\delta}{\gamma}$. Thus, by using Theorem 7.1 in Chapter 2 in ^[35], we get that when $\delta\neq \beta\gamma$, the equilibrium $E_1$ is a saddle node.

When $\delta = \beta\gamma$, (2.15) can be written as

By Theorem 7.1 in ^[35] again, we obtain $m = 3$ and $a_m = \frac{\beta}{\gamma} > 0$. Then $E_0$ is a non-hyperbolic unstable node. Hence, $E_1$ is a non-hyperbolic stable node due to that we have used the transformation $\tau = -t$. Accordingly, Theorem 2.2 is proved.

Theorem 2.3. For the boundary equilibria $E_{21}$, $E_{22}$ and $E_{23}$, the following statements are true.

$1)$ Assume $\alpha = \alpha^*$ and $\delta-\beta\gamma > 0$, then there exists the boundary equilibrium $E_{23}(0, y_{3}^*)$. Moreover,

$i)$ If $\delta-\beta\gamma > 2\beta$, then $E_{23}$ is a saddle-node, which includes a stable parabolic sector.

$ii)$ If $0 < \delta-\beta\gamma < 2\beta$, then $E_{23}$ is a saddle-node, which includes an unstable parabolic sector.

$iii)$ If $\delta-\beta\gamma = 2\beta$, then $E_{23}$ is non-hyperbolic.

$2)$ Assume $\alpha^{**} < \alpha < \alpha^*$ and $\delta-\beta\gamma > 0$, then there are two boundary equilibria $E_{21}(0, y_{1}^*)$ and $E_{22}(0, y_{2}^*)$. Moreover,

$i)$ If $\alpha < \delta+\delta\gamma-\beta-\beta\gamma$, then $E_{21}$ is a hyperbolic unstable node and $E_{22}$ is a hyperbolic stable node.

$ii)$ If $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $\delta-\beta\gamma > 2\beta$, then $E_{21}$ is a hyperbolic saddle and $E_{22}$ is a hyperbolic stable node.

$iii)$ If $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, then $E_{21}$ is a hyperbolic unstable node and $E_{22}$ is a hyperbolic saddle.

$iv)$ If $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$ and $\delta-\beta\gamma > 2\beta$, then $E_{21}$ a saddle-node, which includes an unstable parabolic sector and $E_{22}$ is a hyperbolic stable node.

$v)$ If $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, then $E_{21}$ a hyperbolic unstable node and $E_{22}$ is a saddle-node, which includes a stable parabolic sector.

$3)$ Assume $\alpha = \alpha^{**}$ and $\delta-\beta\gamma > 0$, then there is one boundary equilibrium $E_{22}(0, y_{2}^*)$, where $y_{2}^* = \frac{\delta-\beta\gamma}{\beta}$. Moreover,

$i)$ If $0 < \delta-\beta\gamma < \beta$, then $E_{22}$ is a hyperbolic saddle.

$ii)$ If $\delta-\beta\gamma = \beta$, then $E_{22}$ is a saddle-node, which includes a stable parabolic sector.

$iii)$ If $\delta-\beta\gamma > \beta$, then $E_{22}$ is a hyperbolic stable node.

$4)$ Assume $0 < \alpha < \alpha^{**}$, then there is one boundary equilibrium $E_{22}(0, y_{2}^*)$.

$i)$ If $\alpha < \delta+\delta\gamma-\beta-\beta\gamma$, then $E_{22}$ is a hyperbolic stable node.

$ii)$ If $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$, then $E_{22}$ is a saddle-node, which includes a stable parabolic sector.

$iii)$ If $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$, then $E_{22}$ is a hyperbolic saddle.

Proof. For the equilibrium $E_{2i}$($i$ = 1, 2, 3), the Jacobian matrix is

where $H(x_i^*, y_i^*) = \delta-2\beta y_i^*-\frac{\alpha\gamma}{(\gamma+y_i^*)^2}$. The two eigenvalues of the above matrix $J(E_{2i})$ are $\lambda_1(E_{2i}) = 1-y_i^*$ and $\lambda_2(E_{2i}) = H(x_i^*, y_i^*)$.

Because $(x_i^*, y_i^*)$ satisfies $\delta-\beta y_i^*-\frac{\alpha }{\gamma+y_i^*} = 0$, we can derive

Then

and

In addition, for $\lambda_1(E_{2i}) = 1-y_i^*$, we discuss it in the following four cases.

Case 1. $\alpha = \alpha^{*}$ and $\delta-\beta\gamma > 0$.

If $\delta-\beta\gamma > 2\beta$, then $y_3^* > 1$, which implies $\lambda_1(E_{23}) = 1-y_1^* < 0$, so $E_{23}$ is non-hyperbolic. In order to determinate the stability of the equilibrium $E_{23}$, we translate $E_{23}$ to the origin by letting $(\overline{x}, \overline{y}) = (x, y-y_3^*)$, and expand the system in power series up to the third order around the origin, under which system (1.4) can be transformed into

where $c_0 = \delta-2\beta y_3^* -\frac{\alpha\gamma}{(\gamma+y_3^*)^2} = 0$, $c_1 = \frac{\alpha\gamma}{(\gamma+y_3^*)^3}-\beta$, $c_2 = -\frac{\alpha\gamma}{(\gamma+y_3^*)^4}$, and $Q_2(y)$ represents for a power series with terms $\overline{y}^i$ $(i\geq 4)$.

Then introducing a new time variable $\tau = (1-y_3^*)t$, we get

We see that the coefficient at $\overline{y}^2$ is $\frac{c_1}{1-y_3^*} > 0.$ Hence, from Theorem 7.1 in Chapter 2 of ^[35], we have $E_{23}$ is a saddle-node, which includes a stable parabolic sector and this parabolic sector is on the upper half-plane.

If $\delta-\beta\gamma < 2\beta$, then $y_3^* < 1$, which means $\lambda_1(E_{23}) = 1-y_1^* > 0$. Same analysis as the above we can get $E_{23}$ is also a saddle-node, which includes an unstable parabolic sector and the parabolic sector is on the lower half-plane.

If $\delta-\beta\gamma = 2\beta$, then $y_3^* = 1$, which means $\lambda_1(E_{23}) = 1-y_1^* = 0$, so $E_{23}$ is a non-hyperbolic critical point with two zero eigenvalues.

Case 2. $\alpha^{**} < \alpha < \alpha^{*}$ and $\delta-\beta\gamma > 0$.

1) If $\alpha < \delta+\delta\gamma-\beta-\beta\gamma$, then $y_1^* < 1 < y_2^*$, that implies $\lambda_1(E_{21}) = 1-y_1^* > 0$ and $\lambda_1(E_{22}) = 1-y_2^* < 0$, so $E_{21}$ is a hyperbolic unstable node and $E_{22}$ is a hyperbolic stable node;

2) If $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $\delta-\beta\gamma > 2\beta$, then $1 < y_1^* < y_2^*$, implying $\lambda_1(E_{21}) = 1-y_1^* < 0$ and $\lambda_1(E_{22}) = 1-y_2^* < 0$, thus $E_{21}$ is hyperbolic saddle and $E_{22}$ is hyperbolic stable node;

3) If $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $\delta-\beta\gamma < 2\beta$, then $y_1^* < y_2^* < 1$, that means $\lambda_1(E_{21}) = 1-y_1^* > 0$, $\lambda_1(E_{22}) = 1-y_2^* > 0$, hence $E_{21}$ is a hyperbolic unstable node and $E_{22}$ is a hyperbolic saddle;

4) If $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$ and $\delta-\beta\gamma > 2\beta$, then $1 = y_1^* < y_2^*$, which means $\lambda_1(E_{21}) = 1-y_1^* = 0$ and $\lambda_1(E_{22}) = 1-y_2^* < 0$, so $E_{21}$ is non-hyperbolic and $E_{22}$ is a hyperbolic stable node.

In order to determinate the stability of the equilibrium $E_{21}$, we translate the equilibrium $E_{21}$ to the origin by the translation $(\overline{x}, \overline{y}) = (x, y-1)$, and expand system (1.4) in power series up to the third order around the origin, which makes system (1.4) to be the following form:

where $d_0 = \delta-2\beta -\frac{\alpha\gamma}{(\gamma+1)^2}$, $d_1 = \frac{\alpha\gamma}{(\gamma+1)^3}-\beta$, $d_2 = -\frac{\alpha\gamma}{(\gamma+1)^4}$, and $Q_3(\overline{y})$ represents for a power series with terms $\overline{y}^i$ satisfying $i\geq 4$.

Then introducing a new time variable $\tau = d_0t$, we get

From $\frac{d\overline{y}}{d\tau} = 0$, we can derive a unique implicit function $\overline{y} = \phi(\overline{x}) = 0,$ which satisfies $\phi(0) = \phi'(0) = 0$ and $\phi(\overline{x}) + P(\overline{x}, \phi(\overline{x})) = 0$. Substituting $\overline{y} = \phi(\overline{x}) = 0$ into the first equation of (2.20), we have that

The coefficient at $\overline{x}^2$ is

By Theorem 7.1 in Chapter 2 of ^[35], we know $E_{21}$ is a saddle-node, which includes an unstable parabolic sector and this parabolic sector is on the left half-plane.

5) If $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$ and $\delta-\beta\gamma < 2\beta$, then $y_1^* < y_2^* = 1$, implying $\lambda_1(E_{21}) = 1-y_1^* > 0$ and $\lambda_1(E_{22}) = 1-y_2^* = 0$, therefore $E_{21}$ a hyperbolic unstable node and $E_{22}$ is non-hyperbolic. Similarly to the proof of $E_{21}$ in 4), we can get that $E_{22}$ is a saddle-node, which includes a stable parabolic sector.

Case 3. $\alpha = \alpha^{**}$ and $\delta-\beta\gamma > 0$.

If $\delta-\beta\gamma > \beta$, then $y_2^* > 1$, that is $\lambda_1(E_{22}) = 1-y_2^* < 0$, so $E_{22}$ is a hyperbolic stable node. If $\delta-\beta\gamma < \beta$, then $y_2^* < 1$, that is $\lambda_1(E_{22}) = 1-y_2^* > 0$, thus $E_{22}$ is a hyperbolic saddle. If $\delta-\beta\gamma = \beta$, then $y_2^* = 1$, that is $\lambda_1(E_{22}) = 1-y_2^* = 0$, thus $E_{22}$ is a non-hyperbolic. Similar to the proof of $E_{21}$ in 4), we can obtain $E_{22}$ is a saddle-node, which includes a stable parabolic sector.

Case 4. $\alpha < \alpha^{**}$.

If $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$, then $y_2^* < 1$, which implies $\lambda_1(E_{22}) = 1-y_2^* > 0$, thus $E_{22}$ is a hyperbolic saddle. If $\alpha < \delta+\delta\gamma-\beta-\beta\gamma$, then $y_2^* > 1$, which means $\lambda_1(E_{22}) = 1-y_2^* < 0$, thus $E_{22}$ is a hyperbolic stable node. If $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$, then $y_2^* = 1$, which implies $\lambda_1(E_{22}) = 1-y_2^* = 0$, thus $E_{22}$ non-hyperbolic. Similar to the proof of 4) in Case 2, we can derive $E_{22}$ is a saddle-node, which includes a stable parabolic sector.

This completes the proof of Theorem 2.3.

Theorem 2.4. For the positive equilibria $E_{31}$, $E_{32}$ and $E_{33}$, the following statements are true.

$1)$ If $\alpha^{**} < \alpha < \alpha^*$, when $\alpha < \delta+\delta\gamma-\beta-\beta\gamma$ and $\delta-\beta\gamma > 0$, or $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, or $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, then $E_{31}$ is a hyperbolic saddle.

$2)$ If $0 < \alpha < \alpha^{**}$ and $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$, or $\alpha = \alpha^{**}$ and $0 < \delta-\beta\gamma < \beta$, or $\alpha^{**} < \alpha < \alpha^*$, $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, then $E_{32}$ is a hyperbolic stable node.

$3)$ If $\alpha = \alpha^*$ and $0 < \delta-\beta\gamma < 2\beta$, then $E_{33}$ is a saddle-node.

Proof. For the equilibrium $E_{3i}$($i$ = 1, 2, 3), the Jacobian matrix is

where $H(x_i^*, y_i^*) = \delta-2\beta y_i^*-\frac{\alpha\gamma}{(\gamma+y_i^*)^2}$. The two eigenvalues of the above matrix $J(E_{3i})$ are $\lambda_1(E_{3i}) = y_i^*-1 < 0$ and $\lambda_2(E_{3i}) = H(x_i^*, y_i^*)$. From the analysis of the proof of Theorem 2.3, it follows that $\lambda_2(E_{31}) > 0$, $\lambda_2(E_{32}) < 0$, and $\lambda_2(E_{33}) = 0$.

Hence, when $\alpha^{**} < \alpha < \alpha^*$, $\alpha < \delta+\delta\gamma-\beta-\beta\gamma$ and $\delta-\beta\gamma > 0$, or $\alpha^{**} < \alpha < \alpha^*$, $\alpha = \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, or $\alpha^{**} < \alpha < \alpha^*$, $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, we have $\lambda_1(E_{31}) < 0$ and $\lambda_2(E_{31}) > 0$, which means $E_{31}$ is a hyperbolic saddle.

When $0 < \alpha < \alpha^{**}$ and $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$, or $\alpha = \alpha^{**}$ and $0 < \delta-\beta\gamma < \beta$, or $\alpha^{**} < \alpha < \alpha^*$, $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$, there exists the positive equilibrium $E_{32}$ with two eigenvalues $\lambda_1(E_{32}) < 0$ and $\lambda_2(E_{32}) < 0$, hence $E_{32}$ is a hyperbolic stable node.

When $\alpha = \alpha^*$ and $0 < \delta-\beta\gamma < 2\beta$, there exists the positive equilibrium $E_{33}$ with two eigenvalues $\lambda_1(E_{33}) < 0$ and $\lambda_2(E_{33}) = 0$. Then we further to study the stability of $E_{33}$ by transforming $E_{33}$ to the origin by the translation $(\overline{x}, \overline{y}) = (x - x_3^*, y-y_3^*)$ and expand system (1.4) in power series up to the third order around the origin, which makes system (1.4) to be the following form:

where $c_0 = \delta y_3^*-\beta (y_3^*)^2-\alpha+\frac{\alpha\gamma}{\gamma+y_3^*} = 0,$ $c_1 = \delta-2\beta y_3^*-\frac{\alpha\gamma}{(\gamma+y_3^*)^2} = 0$, and $Q_4(\overline{y})$ stands for a power series with terms $\overline{y}^i$ $(i\geq 4)$.

In order to transform the Jacobian matrix into a standard form, we use the invertible translation

then system (2.22) can be expressed as

Then introducing a new time variable $\tau = -x_3^*t$, we get

According to the implicit function theorem, from $\frac{du}{d\tau} = 0$, we can derive a unique function

which satisfies $\varphi_4(0) = \varphi_4'(0) = 0$ and $\varphi_4(v) + P(\varphi_4(v), v) = 0$. Substituting it into the second equation of (2.25), we have that

The coefficient at the term $v^2$ is

Therefore, on basis of Theorem 7.1 in ^[35], the equilibrium $E_{33}$ of system (1.4) is a saddle-node.

The proof of Theorem 2.4 is finished.

For readers' convenience, the local dynamical properties of equilibria of system (1.4) are totally summarized in Table 1. It can be clearly seen that the number and stability of equilibria are different according to the value range of the parameters. Therefore, we divide $(\alpha, \beta, \gamma, \delta)\in R_+^{4}$ into the following twelve regions.

For $R_1 = R_{11}\cup R_{12}$, system (1.4) possesses two equilibria $E_0$ and $E_1$. In $R_{11}$, $E_0$ is a saddle and $E_1$ is a stable node. In $R_{12}$, $E_0$ and $E_1$ represent saddle-nodes. And specifically, in $R_1$, $E_1$ is stable in the first quadrant. Therefore, the qualitative properties of these two equilibria can be seen in Figure 1(a).

For $R_2 = R_{21}\cup R_{22}$, there are three equilibria for system (1.4). In $R_{21} = R_{2a}\cup R_{2b}$, system (1.4) admits $E_0$ (unstable node or saddle-node), $E_1$ (saddle or saddle-node) and $E_{22}$ (stable node or saddle-node). In particular, in $R_{21}$, $E_0$ is unstable and $E_{22}$ is stable in the first quadrant. In $R_{22}$, there are $E_0$ (saddle), $E_1$ (stable node) and $E_{23}$ (non-hyperbolic or saddle-node). The qualitative properties of these equilibria in $R_2$ are displayed by Figure 1(b), (c).

For $R_3 = R_{31}\cup R_{32}\cup R_{33}$, there exist four equilibria for system (1.4). Considering $R_{31} = R_{3a}\cup R_{3b}$, there are $E_0$ (unstable node or saddle-node which is unstable in the first quadrant), $E_1$ (saddle or saddle-node), $E_{22}$ (saddle) and $E_{32}$ (stable node). Considering $R_{32}$, system (1.4) admits $E_0$ (saddle), $E_1$ (stable node), $E_{21}$ (saddle or saddle-node) and $E_{22}$ (stable node). Considering $R_{33}$, there exist $E_0$, $E_1$, $E_{23}$ and $E_{33}$, which are saddle, stable node, saddle-node and saddle-node. The detailed image descriptions of the qualitative properties of the equilibria in $R_3$ are illustrated by Figure 1(d)–(f).

For $R_4 = R_{41}\cup R_{42}$, system (1.4) has five equilibria $E_0$ (saddle), $E_1$ (stable node), $E_{21}$ (unstable node), $E_{22}$ (stable node or saddle-node) and $E_{31}$ (saddle). Particularly, $E_{22}$ is stable in the first quadrant. As shown in Figure 1(g), the qualitative properties of these five equilibria are given.

For $R_5$, system (1.4) possesses six equilibria $E_0$, $E_1$, $E_{21}$, $E_{22}$, $E_{31}$ and $E_{32}$, which represent saddle, stable node, unstable node, saddle, saddle and stable node, respectively. Hence, the qualitative properties of these six equilibria are displayed by Figure 1(h).

2.3.   Bifurcation analysis

In this section, all possible bifurcations of system $(1.4)$ will be found out. Moreover, the feasible conditions for the occurence of those bifurcations will be given.

Theorem 2.5. System $(1.4)$ undergoes a saddle-node bifurcation at equilibrium $E_{33}$ when the parameters satisfy the restriction $\alpha = \alpha_{SN} = \frac{(\delta+\beta\gamma)^2}{4\beta}$ along with the conditions $0 < \delta-\beta\gamma < 2\beta$ and $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$ as given in Theorem 2.1.

Proof. Now we will verify the transversality condition for the occurrence of a saddle-node bifurcation at $\alpha = \alpha_{SN}$ by utilizing the Sotomayor's theorem ^[36]. From Section 2.2, we have the two eigenvalues of $J(E_{33})$ are $\lambda_1(E_{33}) < 0$ and $\lambda_2(E_{33}) = 0$. Denote $V$ and $W$ the eigenvectors corresponding to the eigenvalue $\lambda_2(E_{33})$ for the matrices $J (E_{33})$ and $J (E_{33})^T$, respectively. Simple computations yield

Furthermore, we can obtain

and

Obviously, when $0 < \delta-\beta\gamma < 2\beta$, the vectors $V$ and $W$ satisfy the transversality conditions

and

Hence, by Sotomayor's theorem, when $\alpha = \alpha_{SN}$, system (1.4) undergoes a saddle-node bifurcation at non-hyperbolic critical point $E_{33}$. The number of positive equilibria of system (1.4) changes from zero to two as $\alpha$ passes from the right of $\alpha = \alpha_{SN}$ to the left. Thus, the proof of Theorem 2.5 is completed.

Theorem 2.6. System $(1.4)$ undergoes a transcritical bifurcation at the boundary equilibrium $E_{1}$ when the parameters satisfy the conditions $\alpha = \alpha_{TC} = \delta\gamma$ and $0 < \delta-\beta\gamma < 2\beta$.

Proof. Now we also apply Sotomayor's theorem ^[36] to prove the transversality condition for the occurrence of transcritical bifurcation at $E_{1}$ as $\alpha = \alpha_{TC}$. From Section 2.2, we know that $\lambda_1(E_{1}) < 0$ and $\lambda_2(E_{1}) = 0$. Let $V$ and $W$ represent the two eigenvectors corresponding to the zero eigenvalue $\lambda_2(E_{1})$ of the matrices $J (E_{1})$ and $J (E_{1})^T$, respectively, then they are given by

Furthermore, we can obtain

and

Clearly, the vectors $V$ and $W$ satisfy the transversality conditions

and

Therefore, by Sotomayor's theorem, if $\delta\neq \beta\gamma$, system (1.4) experiences a transcritical bifurcation at non-hyperbolic critical point $E_{1}$ as the parameter $\alpha$ varies through the bifurcation value $\alpha = \alpha_{TC}$.

Thus, the proof of Theorem 2.6 is completed.

2.4.   Global structure

In this subsection, we pay some attention to the global structure of the plane nonlinear dynamic system (1.4).

2.4.1.   The dynamics near infinity

For the sake of the research of the global dynamics for system (1.4), we have to investigate the qualitative properties of equilibria at infinity, which is important and useful to study the behavior of orbits when $|x|+|y|$ $\rightarrow$ $\infty$.

Firstly, using Poincaré transformation of Chapter 5 in ^[35]:

system $(1.4)$ can be transformed into

In the first quadrant of $u$-$z$ plane, setting $z = 0$, system $(2.30)$ has one boundary equilibrium $e\left(\frac{1}{\beta-1}, 0\right)$ if $\beta > 1$ and no boundary equilibrium if $\beta \leq 1$ on the nonnegative $u$-axis.

Proceed to the next step, we apply the second type of Poincaré transformation:

system $(1.4)$ becomes

In the nonnegative cone of $v$-$z$ plane, there exist two boundary equilibria $e_1(0, 0)$ and $e_2(\beta-1, 0)$ if $\beta > 1$ and only one $e_1(0, 0)$ if $\beta \leq 1$ for system $(2.31)$.

In this paper, of concern is the global dynamics of system (1.4) when $\beta > 1$. The case $\beta\leq1$ is similar, hence, it is omitted.

Theorem 2.7. In system $(2.30)$, $e\left(\frac{1}{\beta-1}, 0\right)$ is a saddle. In system $(2.31)$, $e_1(0, 0)$ is an unstable node and $e_2(\beta-1, 0)$ is a saddle.

Proof. The Jacobian matrix at $e\left(\frac{1}{\beta-1}, 0\right)$ is evaluated as follows:

obviously, which has two eigenvalues $\lambda_1(e) = -1 < 0$ and $\lambda_2(e) = \frac{\beta}{\beta-1} > 0$. Hence, $e$ is a saddle in system $(2.30)$.

Similarly, the corresponding Jacobian matrices of system $(2.31)$ evaluated at $e_1 (0, 0)$ and $e_2 (\beta-1, 0)$ are

Clearly, $e_1$ is an unstable node and $e_2$ is a saddle in system $(2.31)$. This proof of Theorem 2.7 is completed.

From the point of view of Poincaré transformation, we get

Consequently, $e$ and $e_2$ are equivalent to an infinite singular point $I$ for system $(1.4)$, $e_1$ is equivalent to an infinite singular point $I_1$ for system $(1.4)$. By combining the Poincaré transformation of Chapter 5 in ^[35] and Theorem 2.7, we illustrate the dynamical behaviors near infinity of model $(1.4)$ with Figure 2.

2.4.2.   Nonexistence of closed orbit

We, here, for further research of the global structure of system (1.4), need to discuss the existence of closed orbits in phase space.

Theorem 2.8. No closed orbit exists for model $(1.4)$ in $R_+^{2}$.

Proof. On the contrary, assume that model (1.4) exhibits a closed orbit in $R_+^{2}$. From the discussion in Theorem 2.1, we can deduce that all positive equilibria will lie on invariant lines $y = \frac{\delta-\beta\gamma- \sqrt{\Delta}}{2\beta}$ or $y = \frac{\delta-\beta\gamma+\sqrt{\Delta}}{2\beta}$. On the basis of Theorem 4.6 of ^[35], in the interior of the closed orbit, system (1.4) must admit a positive equilibrium. As such, the closed orbit has two points of intersection with $y = \frac{\delta-\beta\gamma- \sqrt{\Delta}}{2\beta}$ or $y = \frac{\delta-\beta\gamma+\sqrt{\Delta}}{2\beta}$, which is contrary to the existence and uniqueness of solutions. Therefore, there is no closed orbit for model (1.4). We end the proof of Theorem 2.8.

2.4.3.   Global phase portraits

System (1.4), according to the aforementioned Theorem 2.7, admits two equilibria $I$ and $I_1$ at infinity, which stand for a saddle and an unstable node, respectively. A step further, we derive that model (1.4) exists no closed orbit in $R_+^{2}$ from the above Theorem 2.8. It now follows that all possible cases of global phase portraits for model (1.4) can be given.

For the purpose of simplicity, $\xi_I$ and $\omega(\xi_I)$ are utilized to refer the unstable manifold of $I$ and the $\omega$-limit set of $\xi_I$ in the nonnegative cone of $R^{2}$, respectively.

For $\left(\alpha, \beta, \gamma, \delta \right)$ $\in$ $R_1$, it is easy to be checked that $\omega(\xi_I)$ is $E_1$. See Figure 3(a) for the global phase portrait of model (1.4) related to the case $R_1$.

For $\left(\alpha, \beta, \gamma, \delta \right)$ $\in$ $R_2$, we can easily verify that $\omega(\xi_I)$ is $E_{22}$ in $R_{21}$ and $\omega(\xi_I)$ is $E_{23}$ in $R_{22}$. See Figure 3(b), (c) for the global phase portraits of model (1.4) related to the cases $R_{21}$ and $R_{22}$.

For $\left(\alpha, \beta, \gamma, \delta \right)$ $\in$ $R_3$, we analyze the global phase portraits of model (1.4) from the three subregions $R_{31}$, $R_{32}$ and $R_{33}$, respectively. When $\left(\alpha, \beta, \gamma, \delta \right)$ $\in$ $R_{31}$, we can verify that $\omega(\xi_I)$ is $E_{32}$. The global phase portraits of model (1.4) related to $R_{31}$ are given by Figure 3(d), (e). When $\left(\alpha, \beta, \gamma, \delta \right)$ $\in$ $R_{32}$, $\omega(\xi_I)$ is $E_{22}$. The global phase portrait of model (1.4) related to $R_{32}$ is shown by Figure 3(f). When $\left(\alpha, \beta, \gamma, \delta \right)$ $\in$ $R_{33}$, it is clear that $\omega(\xi_I)$ is $E_{33}$ and the global phase portraits of model (1.4) in this region are shown in Figure 3(g), (h).

For $\left(\alpha, \beta, \gamma, \delta \right)$ $\in$ $R_4$, it is clear that $\omega(\xi_I)$ is $E_{22}$. See Figure 3(i) for the global phase portrait of model (1.4) related to the case $R_4$.

For $\left(\alpha, \beta, \gamma, \delta \right)$ $\in$ $R_5$, it is also clear that $\omega(\xi_I)$ is $E_{32}$. And the global phase portraits of model (1.4) in $R_5$ are shown in Figure 3(j), (k).

3.   Global dynamics of system $(1.1)$

Now, we mainly consider the dynamic properties of all possible equilibria and bifurcation behaviors of model $(1.1)$. More, in $R_+^{2}$, we analyze the global structure of this model and illustrate our analysis with the global phase portraits.

To facilitate this, system $(1.1)$ takes the following form

On account of the dynamic properties of system $(3.1)$ corresponding to system $(1.1)$, we study the dynamics of system $(3.1)$ in the following.

3.1.   Existence and stability of equilibria

As it has been shown in ^[5], the existence and stability properties of equilibria of system $(3.1)$ are clearly summarized in Table 2. Next, we divide $(\delta, \beta) \in R_+^{2}$ into two regions as follows.

Consider $(\delta, \beta) \in G_1$, system $(3.1)$ has three equilibria $E_0$, $E_1$ and $E_2$. Specifically, $E_0$ and $E_1$ represent an unstable node and a saddle. $E_2$ is a stable node in $G_{11}$ and non-hyperbolic in $G_{12}$, which is stable in the first quadrant. See Figure 4(a), the qualitative properties of these three equilibria are illustrated.

Consider $(\delta, \beta) \in G_2$, system $(3.1)$ admits four equilibria $E_0$, $E_1$, $E_2$ and $E_3$, which are unstable node, saddle, saddle and stable node, respectively. See Figure 4(b), the qualitative properties of these four equilibria are illustrated.

3.2.   Bifurcation analysis

Here, we will analyze some of all the possible bifurcation scenarios for model $(3.1)$. According to the previous methods applied in Section 2.3, a heuristic discussion similar to the one used by Theorem 2.6 can easily get the conditions for a transcritical bifurcation. Therefore, the following theorem can be directly obtained.

Theorem 3.1 System $(3.1)$ undergoes a transcritical bifurcation at $E_2$ when the parameters satisfy the condition $\delta = \beta$.

3.3.   Global structure

In order to investigate the global structure of system $(3.1)$, the methods proposed by Section 2.4 can be extended. By applying Poincaré transformation and analysis of closed orbit, Figure 5(a) describes the existence and stability of the singular points at infinity, as well as all the boundary equilibria. In addition, Figure 5(b)–(d) describes minutely the global phase portraits in $R_+^{2}$ similarly.

4.   Numerical simulations and discussions

Some numerical examples, in this section, are proposed so as to check the obtained parameter conditions and theoretical results of systems $(1.4)$ and $(3.1)$. Further, we are interested to show the role of harvesting on the local dynamical properties.

Example 4.1. We consider parameter values as below: $\alpha$ = 2, $\beta$ = 1.2, $\gamma$ = 0.8, $\delta$ = 2.5. Simple calculations lead to $\alpha = \alpha^{**}$, $\delta-\beta\gamma \geq \beta$ and $\delta > \beta$, which means both systems (1.4) and $(3.1)$ have three boundary equilibria. The numerical phase portraits related to this case are shown in Figures 6(a) and 7(a). The impact of harvesting is shown by Figure 8.

Example 4.2. Some parameter values are taken $\alpha$ = 1.85, $\beta$ = 0.3, $\gamma$ = 2.5 and $\delta$ = 0.8. In this case, parameters satisfy the conditions $0 < \alpha < \alpha^{**}$ and $\alpha > \delta+\delta\gamma-\beta-\beta\gamma$. Consequently, system (1.4) has a positive equilibrium (global asymptotically stable) and three boundary equilibria. From $\delta > \beta$, we can effortlessly observe there exists no positive equilibrium for model $(3.1)$. See Figures 6(b) and 7(b), the numerical phase portraits corresponding to this situation are presented. The impact of harvesting is illustrated with curves in Figure 9.

Example 4.3. We select parameter values as follows: $\alpha$ = 0.8, $\beta$ = 2.5, $\gamma$ = 0.4, $\delta$ = 2. For these parameters, we derive $\alpha = \alpha^{**}$ and $0 < \delta-\beta\gamma < 2\beta$. Therefore, system (1.4) possesses the same types of equilibria as the above Example 4.2. Whereas, there is a positive equilibrium (global asymptotically stable) in system $(3.1)$ for $\delta < \beta$. See Figures 6(c) and 7(c), the corresponding numerical phase portraits are presented. The impact of harvesting is illustrated with curves in Figure 10.

It can be seen that the harvesting effect has a significant role on the population densities. From a biological point of view, the permanence and extinction of species are influenced to differential extent by the harvesting effect on the second species, for instance, greatly slowing the extinction of the first species (Figure 8) or guaranteeing the permanence of the system by preserving the first species from extinction (Figure 9). In addition, both species need much more time to reach their stable coexistence state due to the harvesting effect, which is shown in Figure 10.

5.   Conclusions

In this article, by combining local and global dynamical analysis, we have investigated the dynamics of the amensalism system with Michaelis-Menten type harvesting for the second species. In order to further explore the influence of harvesting effect, a complete qualitative analysis of the system with no harvesting is also performed. We could show some differences for the model (1.4) with harvesting on the second species and the system (3.1) with no harvesting.

$1)$ When harvesting is present in the system, it has shown that model (1.4) has more different types of equilibria. System (3.1) has at most four equilibria including one interior point in ${R}^{2}_+$ and the unique interior point (if it exists) of system (3.1) is globally stable. However, when the Michaelis-Menten type harvesting is on the second species, we shown that the model (1.4) has at most six equilibria including two interior points in ${R}^{2}_+$ and the interior point $E_{31}$ is a saddle.

$2)$ The complex existence of equilibria leads to more bifurcation behaviors of model (1.4). Two kinds of bifurcation, under some parameter restrictions, saddle-node bifurcation and transcritical bifurcation will occur for system (1.4), while there exhibits only one saddle-node bifurcation for the system with no harvesting.

$3)$ More types of equilibria and more bifurcation behaviors can induce potentially dramatic changes to the dynamics of the system. By the classifications of global phase portraits of model (1.4), we find that more complex dynamics occur for the appearance of harvesting on the second species. Such as, for model (1.4) with harvesting, the first species survives permanently rather than extinction or both species take a lot longer to approach the coexistence state.

These results reveal that the dynamical properties of system (1.4) get more complicated and richer compared to the system with no harvesting. Furthermore, it would also be interesting to study system (3.1) with both the first species and the second species with harvesting terms since we usually harvest, or would like to harvest, both populations.

Acknowledgments

The first author was supported by the National Natural Science Foundation of China (No.12001503) and the third author was supported by the Project of Beijing Municipal Commission of Education (KM 202110015001).

Conflict of interest

The authors declare there is no conflicts of interest.