Qualitative Analysis of a Single-Species Model with Distributed Delay and Nonlinear Harvest

In this paper, a single-species population model with distributed delay and Michaelis-Menten type harvesting is established. Through an appropriate transformation, the mathematical model is converted into a two-dimensional system. Applying qualitative theory of ordinary differential equations, we obtain sufficient conditions for the stability of the equilibria of this system under three cases. The equilibrium A1 of system is globally asymptotically stable when br−c>0 and η<0. Using Poincare-Bendixson theorem, we determine the existence and stability of limit cycle when br−c>0 and η>0. By computing Lyapunov number, we obtain that a supercritical Hopf bifurcation occurs when η passes through 0. High order singularity of the system, such as saddle node, degenerate critical point, unstable node, saddle point, etc, is studied by the theory of ordinary differential equations. Numerical simulations are provided to verify our main results in this paper.


Introduction
Single-species is a unit of the whole ecosphere. Although most biological populations are multi-species, there is no single population in the strict sense, but because of artificial breeding, there are a lot of single population resources in many human-created environments, which provide indispensable roles for human production and life. Due to the economic interests and production needs, people need to develop the single population resources for a long time and continuously, and to give full play to the best use value on the basis of the lowest cost consumption as possible, so the control and prediction of the single population is particularly important.
Single-species models are widely used in the field of mathematical biology, such as pest management, optimal management of biological resources, epidemic prevention and control, cell growth regulation and so on. Many scholars (see [1][2][3]) have obtained some good properties of these models by quantitative analysis, and the results help us forecast and control the actual production. Wang et al. [4] studied a single-species model and found an optimal harvesting strategy which allows the output to reach a maximum and remains constant. At the same time, the population quantity can attain the maximum level in a precise time interval when the population has been harvested. Dou et al. [5] addressed a non-autonomous Logistic single-species model. The analytical expressions about the best harvesting strategy have been obtained by Pontryagin control and optimization principles of impulsive systems. Recently, many articles have also focused on the study of singlespecies population model (see [6][7][8][9]).
We know that both discrete and distributed delays tend to bring more rich and complex dynamics to a population dynamic system. For researchers, although models with delay are more complex than without delay, the research results of these models are closer to actual life. In 1948, Hutchinson [10] analyzed a Logistic delay equation characterizing animal populations, and the form is which could vividly simulate the population size. In (1), the food supply at time t depends on the population number at time t − τ. More generally, many scholars studied some models in which the current population density continuously depends on the population density in the past period. For example, Cushing [11] proposed the following single population model with distributed delay and detailed qualitative results were obtained. Pang et al. [12] considered the following single population model with distributed delay and impulsive state feedback control and sufficient conditions for the existence and stability of the order-1 periodic solution and limit cycles were obtained, where ∆x(t) = x(t + ) − x(t). Ruan and Wolkowicz [13] built a chemostat model about single population with distributed delay, and they obtained that the system existed Hopf bifurcations by taking the average time delay as the bifurcation parameter. Lian et al. [14] considered a predator-prey system with Holling type IV functional response and time delay. They studied the effects of time delay on this system and obtained the conditions of local asymptotic stability of the positive equilibrium and the existence of local Hopf bifurcations by applying the delay as a bifurcation parameter. Yao et al. [15] investigated the global asymptotic stability of fractional-order complex-valued differential equations with distributed delays. Based on the Laplace transform method, a novel necessary and sufficient condition for the stability is established by embedding the characteristic equation into a two-dimensional complex system. The algebraic criterion is expressed by the fractional exponent, coefficients, and the delay. Many single-species population models with time delay were studied and rich results were obtained (see [16][17][18][19][20][21][22]). In addition, many scholars have studied the single population model with different forms of harvest. Clark [23] introduced the following single population model with linear harvest, Clark introduced the single population model (4) with linear harvest, analytical expressions for optimal harvest of a renewable resource stock, which is subject to a stochastic process are found. These expressions give the optimal harvest as an explicit feedback control law. All relations in system (4), including the stochastic process, may be arbitrary functions of the state variable (stock). In 1978, Ludwig et al. [24] put forward the following single population model with nonlinear harvest, Nonlinear harvesting has been proposed and some novel conclusions have been obtained. Tan et al. [25] considered the following single population system with implusive disturbance and nonlinear harvest: They obtained good dynamic properties. Li et al. [26] established and analyzed a singlespecies model with delay weak kernel and constant rate harvesting: The existence of the equilibrium point of the model (7) is obtained, and the properties of the positive equilibrium state of the population are studied.
Population model with distributed delay and harvest plays an important role in ecosystem, which can reflect the growth rule of invertebrate population better and help us get the best policy of harvesting. This paper is mainly aimed at qualitative analysis of a single-species population model with distributed delay and nonlinear harvest. By the establishment of a single-species population model of differential equations, and using stability analysis methods to discuss the existence and stability of equilibria of the system, we can forecast the development of the population.
The rest of this paper is organized as follows. In Section 2, we establish a single-species model with distributed delay and nonlinear harvest and make a statement about the model. Then, we will discuss the existence and stability of the equilibria of the system under three different cases from Sections 3-5. In Section 6, we verify our results by numerical simulations. The conclusions will be given in Section 7.

The Model
We combine distributed delay and nonlinear harvest to consider a new single-species model. From biological and economic points of view, Michaelis-Menten type harvesting is more realistic among the several types of harvesting. Thus, a single-species model with distributed delay and Michaelis-Menten type harvesting is established as follows: where x(t) denotes the density of the population at time t, r denotes the intrinsic rate of increase of the population, ω, α, b and c are positive constants.
x(s)ds represents the effect of distributed delay from the following aspects: The integral for s from −∞ to t represents the continuous influence on current state x(t) from the past time, which shows the influence interval for distributed delay. That is, the entire interval from the past to current time t has impact on current state x(t). The kernel function α exp(−α(t − s))(s ∈ [0, ∞)) represents the influential degree to current state x(t) from time s to t. x(s)(s ∈ (∞, t]) represents the density of population at the past time s. Therefore, the term ω t −∞ α exp(−α(t − s))x(s)ds represents the distributed delay in the model.
b+x(t) is simplified Michaelis-Menten type harvesting (the details of this kind of harvesting can be seen in the Refs. [27][28][29]).
We consider the existence of equilibria of system (9). Let the solutions of the equations are x = y = 0 and we can also write We denote p = ωb − r, q 2 = (r − ωb) 2 + 4ω(br − c) = (r + ωb) 2 − 4ωc and k = r + ωb throughout this paper. In the following discussion, we always assume q ≥ 0. Next, we analyze the existence and stability of the equilibria of system (9) under three cases.

The Existence and Stability of the Equilibria for br − c > 0
Let us consider the existence and stability of the equilibria of system (9) for the case br − c > 0. Obviously, q > 0 and q > p in this case.
Proof. When br − c > 0, the Jacobi matrix of system (9) is At the equilibrium O, the Jacobi matrix of system (9) becomes which allows us to obtain Since br − c > 0, then det(J| O ) < 0. Obviously, O(0, 0) is a saddle point. At the equilibrium A 1 , the Jacobi matrix of system (9) becomes Then we get Since we can also write q 2 = k 2 − 4ωc and p = ωb − r, it is easy for us to simplify det(J| A 1 ) as follows Since br − c > 0, we know that q − p > 0 and q + k > 0, then det(J| A 1 ) > 0. Next, we simplify tr(J| A 1 ) as follows We have tr(J| A 1 ) < 0 for η < 0, tr(J| A 1 ) > 0 as η > 0, tr(J| A 1 ) = 0 when η = 0. We discuss the situation of tr(J| A 1 ) < 0, i.e., η < 0. If br − c > 0, then tr(J| A 1 ) < 0, the positive equilibrium A 1 is locally stable in this case. Next, we are going to prove that there is no closed cycle in Ω. There is a unique positive equilibrium ). Suppose that there is a closed orbit Γ around A 1 , which is, x = x(t), y = y(t), −∞ < t < +∞, so there must exist a periodic T such that x(T) = x(0), y(T) = y(0).
We rewrite system (9) as follows Integrating the two sides along the closed orbit Γ, we get So we know that 1 T T 0 x(t)dt and 1 T T 0 y(t)dt also satisfy the following algebraic equation As The closed orbit Γ is stable as γ 0 < 0, which is contradicted by the fact that the equilibrium A 1 is asymptotically stable. Therefore, there is no closed orbit. Now, we are going to prove that the solution of system (9) is eventually uniformly bounded (see Figure 2). Take a point F 1 on x-axis, draw the trajectory passing through the point F 1 of the following equation and the trajectory intersects isoclinic line r − ωy − c b+x = 0 at F 2 . Comparing the direction of trajectory of system (9) with the tangent direction of trajectory F 1 F 2 of the above system So the trajectory on F 1 F 2 of system (9) goes through into the interior of F 1 F 2 from the outer of F 1 F 2 when the trajectory of system (9) intersects F 1 F 2 . Then draw the trajectory of the following equation passing through the point F 2 , and the trajectory will intersect isoclinic line y = x at F 3 .
Comparing the direction of trajectory of system (9) with the tangent direction of trajectory So the trajectory on F 2 F 3 of system (9) goes through into the interior of F 2 F 3 from the outer of F 2 F 3 when the trajectory of system (9) intersects F 2 F 3 . Then draw the line l from the point F 3 , which is perpendicular to y-axis. The trajectory of system (9) goes through into the interior of l from the outer of l when the trajectory of system (9) intersects line l. Furthermore, y-axis is an isoclinic line, and the trajectory of system (9) goes through into the interior of x-axis from the outer of x-axis when intersects x-axis. According to the above analysis, the solution of system (9) is eventually uniformly bounded in Ω. In addition, there is no closed orbit in Ω. So the unique positive equilibrium Next, we discuss the case of tr(J| A 1 ) = 0, i.e., η = 0. We shift point A 1 (x 1 , y 1 ) of system (9) to the origin and system (9) can be rewritten as where [x] 4 is a C ∞ function of at least fourth order about x. In this case, tr(J| A 1 ) = 0, that is η = 0, then J| A 1 becomes We can obtain n det( So we can know that the eigenvalues of J| A 1 are λ 1 = ni and λ 2 = −ni (i is the complex unit, similarly hereinafter) by simple calculation. Taking non-singular transformations x =x α −ȳ n and y = −ȳ n , and also denotingx,ȳ by x, y, respectively, system (11) becomes where [x, y] 4 is a C ∞ function of at least fourth order about x and y. f (x, y) and g(x, y) are C ∞ functions and satisfy f (0, 0) = g(0, 0) = 0.
Using the formula for the Lyapunov number of Andronov et al. [30], we have the first Lyapunov number α 1 of the equilibrium A 1 of system (9) as follows holds, it follows that α 1 < 0, where n = −α 2 − α p−q 2 . So we obtain the following theorem: If br − c > 0 and condition (12) holds, a supercritical Hopf bifurcation occurs near Then, we discuss the case of tr(J| A 1 ) > 0, i.e., η > 0. When br − c > 0 and η > 0, the positive equilibrium A 1 is an unstable focus or node. Given what has been proved before, it is not difficult to obtain that the solution of system (9) is also uniformly bounded in Ω. Since the unique positive equilibrium A 1 is unstable, we can obtain that there exists at least one limit cycle which is stable in region Ω by Poincaré-Bendixson Theorem.
Particularly, we have known that a Hopf bifurcation occurs when br − c > 0 and condition (12) holds. That is to say, the equilibrium A 1 is not a true center. The derivative of system (9) at the point A 1 is denoted as A(η). We can get the eigenvalues of A(0) are ±ni, i.e., We also define the matrix B(η) which satisfies We obtain that TrB(0) = 0. Using the Friedrich Theorem of Hopf bifurcation, we get the truth that there exists a unique stable limit cycle near A 1 when br − c > 0, η > 0 and condition (12) holds. So we get the following Theorem: If br − c > 0 and η > 0, there exists at least one limit cycle which is stable in Ω. Particularly, there exists a unique stable limit cycle near A 1 when br − c > 0, η > 0 and condition (12) 4. The Existence and Stability of the Equilibria for br − c < 0 Let us consider the existence and stability of the equilibria for the case br − c < 0.
Considering the practical significance, this case will not be discussed (see Figure 3III).
) is a negative equilibrium of system (9), this case will not be discussed (see Figure 3IV).
(6) When br − c < 0 and r − ωb = 0, system (9) has no more equilibrium except O(0, 0) (see Figure 3V).   (9) is globally asymptotically stable. If br − c < 0, r − ωb > 0 and q 2 > 0, the positive equilibrium A 2 of system (9) is a saddle point. If br − c < 0, r − ωb > 0, q 2 > 0 and β < 0, the positive equilibrium A 3 of system (9) is Proof. Firstly, we analyze the stability of the equilibrium O(0, 0) of system (9). Through the preceding discussion, we have When br − c < 0, it is easy to get det(J| O ) > 0 and tr(J| O ) < 0. So the equilibrium O is locally stable. Particularly, when br − c < 0 and r − ωb = 0, System (9) has only one equilibrium O (0, 0). According to what we have proved before, the solution of system (9) is eventually uniformly bounded in Ω under this case. Since x = 0 is the isoclinic line, it is not possible in the whole plane to have a closed orbit which surrounds the equilibrium O. If system (9) has a limit cycle which exists in first quadrant, it must contain the equilibrium A 4 in the limit cycle. However, under this circumstance, system (9) has no other equilibria except O(0, 0), so we have the conclusion that system (9) has no limit cycle when br − c < 0 and r − ωb = 0. Like the proof of Theorem 1, we know that equilibrium O is globally asymptotically stable when br − c < 0 and r − ωb = 0; that is, all trajectories of the solution tend to equilibrium O if br − c < 0 and r − ωb = 0.
Secondly, we discuss the stability of the positive equilibrium A 2 . We substitute It follows from q 2 = k 2 − 4ωc and p = ωb − r that det(J| A 2 ) becomes det(J| A 2 ) = −αq(p + q) q − k .
Since p + q < 0 and q − k < 0, we get Apparently, when br − c < 0, r − ωb > 0 and q 2 > 0, the positive equilibrium A 2 is a saddle point. Finally, we analyze the stability of the positive equilibrium A 3 . We substitute Through simplification, we obtain It follows from the inequalities br − c < 0, r − ωb > 0, q 2 > 0, p − q < 0 and q + k > 0 that det(J| A 3 ) > 0. We also simplify tr(J| A 3 ) to get Suppose that β < 0, then tr(J| A 3 ) < 0, which means the positive equilibrium A 3 is locally stable. Under this circumstance, there are three equilibria for system (9). Therefore, we only consider the local situation of system (9). Since O(0, 0) is locally stable, and A 2 is always a saddle point, and A 3 is locally asymptotically stable, the trajectory of the solution tends to the equilibrium point A 3 when the initial value is selected in the attraction domain of equilibrium point A 3 , or the trajectory of the solution tends to the equilibrium point O when the initial value is selected in the attraction domain of equilibrium point O.
Suppose that β = 0, then tr(J| A 3 ) = 0, which means A 3 is a center-type singular point. Under this circumstance, the trajectory of the solution tends to the equilibrium point O or forms a closed trajectory near A 3 . In this case, system (9) admits the probability that a Hopf bifurcation occurs near A 3 . Since the form of A 3 is the same with A 1 , we use the result of Theorem 2 to obtain the first Lyapunov number α 2 of A 3 : holds, it follows that α 2 < 0. So we get that a supercritical Hopf bifurcation occurs near A 3 . Suppose that β > 0, then tr(J| A 3 ) > 0, which means the positive equilibrium A 3 is unstable. Under this circumstance, if O(0, 0) is locally stable, and A 2 is always a saddle point, and A 3 is unstable, then, there is a stable limit cycle generated by Hopf bifurcation in the neighborhood of A 3 when 0 < β 1. All trajectories of system (9) will be far away from A 3 and tend to O.
Next, we analyze the situation that A 2 and A 3 are merged into a single equilibrium point A 4 . We first give the Theorem 5 and then prove it.
Proof. We substitute A 4 (− p 2ω , − p 2ω ) into (10) to obtain which yields where k = r + ωb and p = ωb − r, it is easy to get det(J| A 4 ) = 0. It can easily be shown that this equilibrium point is a higher order singular point. Now we do this work. We shift this point A 4 (x 4 , y 4 ) of system (9) to the origin, then, system (9) can be written as where [x] 3 is a C ∞ function of at least third order about x. When r − ωb − 2α = 0, (14) changes into Then the eigenvalues of J| A 4 are λ 1 = 0 and λ 2 = r−ωb−2α 2 . At this time, A 4 is a codimension 1 singular point. Using x =x +ȳ, y =x + 2α r−ωbȳ and dt = 2 r−ωb−2α dτ, we also denotex,ȳ and τ by x, y and t, respectively. System (15) becomes Apply dy dt = 0 to get We take the above equality into dx dt and obtain where [x, y] 3 is a C ∞ function of at least third order about x and y. It follows that the right side of the above equality is a polynomial about x and y, and the order of it is at least two. By Theorem 7.1 in Chapter 2 of Zhang et al. [31], we know that A 4 is a saddle node. When r − ωb − 2α = 0, (14) changes into Then the eigenvalues of J| A 4 are λ 1 = 0 and λ 2 = 0. That is, A 4 is a codimension 1 singular point. System (15) becomes Taking the non-singular transformations x =x +ȳ, y =x −ȳ and dt = 1 2α dτ, we also denotex,ȳ and τ by x, y and t, respectively. System (16) becomes By a one-to-one transformation in a sufficiently small neighborhood of the origin x → x, y + P 2 (x, y) → y, We still denoteQ 2 (x, y) by Q 2 (x, y) with terms of order no less than two. We write where h(x), g(x) and f (x, y) are analytic functions, moreover, h(0) = g(0) = 0. Hence the origin of system (18) is a double singularity. By Theorem 7.3 in Chapter 2 and its Corollary of Zhang et al. [31], we find that the origin of system (17) is a degenerate critical point. That is, the equilibrium point A 4 of system (9) is a degenerate critical point.
According to what we have proved before, the solution of system (9) is eventually uniformly bounded in Ω under this case. Since x = 0 is the isoclinic line, it is not possible in the whole plane to have a closed orbit which surrounds the equilibrium O. If system (9) has a limit cycle which exists in first quadrant, there must contain the equilibrium A 4 in the limit cycle, and the index of it must be one. However, the index of equilibrium A 4 is zero whether it is a saddle node or a degenerate critical point, so we have the conclusion that system (9) has no limit cycle when br − c < 0, r − ωb > 0, q 2 = 0. The same to the proof of Theorem 1, we know that equilibrium O is globally asymptotically stable when br − c < 0, r − ωb > 0 and q 2 = 0. That is, all trajectories of the solution tend to equilibrium O when br − c < 0, r − ωb > 0 and q 2 = 0.

The Existence and Stability of the Equilibria for br − c = 0
Let us consider the existence and stability of the equilibria for the case br − c = 0.
Proof. Firstly, we analyze the stability of the equilibrium O(0, 0) of system (9). We take O(0, 0) into (10) and get When br − c = 0, we get det(J| O ) = 0. It can easily be shown that the equilibrium O is a higher order singular point. We change (9) into (19).
Performing the nonsingular transformations x =x, y =x −ȳ and dt = − 1 α dτ, we still denotex,ȳ and τ by x, y and t, and (19) When ω α − c αb 2 = 0, i.e., r − ωb = 0, the positive equilibrium A 5 does not exist under this condition, that is, the equilibrium O is a triple singularity (20) Let y + Ψ 1 (x, y) = 0, then by Taylor expansion at x = 0, we have After some calculations, we have By Theorem 7.1 in Chapter 2 of Zhang et al. [31], we know that O is an unstable node. According to what we have proved before, the solution of system (9) is eventually uniformly bounded in Ω under this case. Since x = 0 is the isoclinic line, it is not possible in the whole plane to have a closed orbit which surrounds the equilibrium O. In addition, if the closed orbit exists in some quadrant, then there must be an equilibrium in the closed orbit. So we have the conclusion that there is no closed orbit in the whole plane when br − c = 0 and r − ωb = 0. Meanwhile, the linear system of system (19) is as follows dx dt = 0, dy dt = αx − αy, and α 2 + (−α) 2 = 2α 2 = 0. It follows that system (19) has a singular line αx − αy = 0, i.e., y = x when br − c = 0 and r − ωb = 0. Then all the trajectories tend to the singular line y = x.
When ω α − c αb 2 = 0, i.e., r − ωb = 0, the positive equilibrium A 5 exists if r − ωb > 0, and the negative equilibrium B 5 exists if r − ωb < 0, that is, the equilibrium O is a double singularity. We rewrite (20) into the following equation Let y + Ψ 2 (x, y) = 0, then By some calculations, we derive By Theorem 7.1 in Chapter 2 of [31], we know that O is a saddle node.
Next, we are going to analyze the stability of the equilibrium A 5 . Obviously, the equilibrium O as a saddle node simultaneously exists with the equilibrium A 5 when br − c = 0 and r − ωb > 0. We first give the Theorem 7 and then prove it.

Theorem 7.
If br − c = 0, r − ωb > 0 and ξ < 0, the unique positive equilibrium point A 5 is locally stable. If br − c = 0, r − ωb > 0 and ξ > 0, the unique positive equilibrium A 5 is unstable and there exists at least one limit cycle in a small neighborhood of A 5 .
When ξ < 0, then tr(J| A 5 ) < 0, which means the positive equilibrium A 5 is locally stable. When ξ > 0, then tr(J| A 5 ) > 0, which means the positive equilibrium A 5 is unstable. In addition, the solution of system (9) is eventually uniformly bounded in Ω, so there exists at least one limit cycle in a small neighborhood of A 5 .
Applying the formula for the Lyapunov number of Andronov et al. [30], we have the first Lyapunov number α 3 of the equilibrium point A 5 of system (9) as follows ( f xxx + f xyy + g xxy + g yyy ) + 1 16m ( f xy ( f xx + f yy ) − g xy (g xx + g yy ) − f xx g xx + f yy g yy )} x=0,y=0 holds, then α 3 < 0. The Theorem 8 is proved.
Applying the Friedrich Theorem of Hopf bifurcation as before, we can know that there exists a unique stable limit cycle near A 5 when br − c = 0, r − ωb > 0 and condition (22) holds.

Numerical Simulation
In the following, applying dde23 in matlab, we verify our theoretical results by numerical simulations for system (9). According to the numerical constraints in the theorems, we first determine the values of b, r and c. Combined with numerical simulations, and try several times to find the value of ω that meets the conditions. For the value of α, it is easy to calculate with the given conditions. All figures correspond to the cases in the table above.

Conclusions
In this paper, we established a single population with distributed delay and nonlinear harvesting. Sufficient conditions of existence and stability of equilibria in three different cases have been derived. We presented the results in Table 1. We choose ω t −∞ α exp(−α(t − s))x(s)ds as a function for the impact of distributed delay, and choose cx(t) b+x(t) as a function for the harvest. Since this type of harvest exists in the ecosphere especially for vertebrates, and single species population model with distributed delay and Michaelis-Menten type harvesting has never been considered and studied by other scholars, our work is novel. Our results show that the qualitative conclusions obtained are rich and have certain guiding significance for some species with nonlinear harvest. In the future work, we will study single-species models with other harvesting types, such as

Conditions
Subconditions Equilibria Stability br − c > 0 O is a saddle, A 1 is globally asymptotically stable.
O is a saddle, A 1 is unstable, there exists at least one limit cycle in Ω.
O is a saddle, A 1 is unstable, there condition (H1) exists a unique limit cycle in Ω.
O is a saddle, A 1 is center-type, a supercritical condition (H1) Hopf bifurcation occurs near A 1 .
O is locally stable, A 2 is a saddle, A 3 is locally stable.
O is locally stable, A 2 is a saddle, A 3 is unstable.
O is locally stable, A 2 is a saddle, a supercritical and condition (H2) Hopf bifurcation occurs near A 3 .
O is a saddle node, A 5 is unstable, there exists at least one limit cycle near A 5 .
O is a saddle node, A 5 is locally stable.
O is a saddle node, A 5 is unstable, a supercritical and condition (H3) Hopf bifurcation occurs near A 5 .
r − ωb = 0 O O is an unstable node, system (9) has a singular line y = x. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: All data generated or analysed during this study are included in this article.