Bifurcation analysis in delayed Nicholson equation with harvest term∗

In this paper, we investigate a delayed Nicholson equation with delay harvesting term which was proposed in open problems and conjectures formulated by Berezansky et al. (Applied Mathematical Modelling 34 (2010) 1405). The stability switching curves by taking two delays as parameters are obtained via the method introduced by An et al.(J. Differential Equations 266 (2019) 7073). The existence of Hopf singularity on a two-parameter plane is determined by the varying direction of two parameters. Furthermore, the normal form near the Hopf singularity is derived via applying the center manifolds theory and normal forms method of FDEs. Finally, some numerical simulations are carried out to illustrate the theoretical conclusions.


Introduction
Blowfly is one of the most important model organisms in biological research due to its short growth cycle [1,2]. In 1954, by conducting some experiments on the Australian sheep blowfly, Nicholson [3] observed oscillations of large amplitude in the blowfly population. Since then, there has been much research on the blowfly dynamics in both the biological and mathematical fields [4][5][6][7][8][9].
To investigate the blowfly population growing in an isolated laboratory, Gurney et. al. [10] proposed a simple time-delay model as follows: where u(t) denotes the population of mature adults at time t, p is the maximum per capita daily egg production rate; 1/a is the size at which the fly population reproduces at its maximum rate; δ represents the per capita adult death rate and τ is the maturation time required for an individual from birth to maturity. They found that system (1) gave good quantitative agreement with Nicholson's classic blowfly data. Besides, the "hump" relationship between future adult population and the current adult population can also be explained. System (1) has subsequently been considered in many fields by researchers, we refer the reader to read the review by Berezansky et al. [4] and the references therein. There has been many studies on model (1) with diffusion, see [11][12][13][14][15] and the references therein.
At the end of the review paper [4], the authors have formulated seven open problems and conjectures. The sixth problem is as following: Assume that a harvesting function is a function of the delayed estimate of the true population. Consider Eq. (1) with a linear harvesting term: with a delayed harvesting strategy. Here H denotes the capture rate, σ represents the capture delay. p, τ, δ are the same as that in Eq. (1). Clearly, Eq.(2) is a delay equation with two delays. A nature question is how do the delays affect the dynamical behaviour of Eq.(2)? The purpose of the paper is to answer this question. At first, we present the Hopf bifurcation curves on a two-parameter plane. We stress that for a possible Hopf singularity, the varying direction of these two delays are important as well as the transversal condition. Then, we calculate the normal form at a Hopf singularity along the varying direction (γ 1 µ, γ 2 µ). The calculation of the normal form with two parameters σ, τ is transformed into the calculation with one parameter µ. This dimension reduction method can reduce the amount of the calculation. Meanwhile, we prove that for some fixed σ, the stability of the positive equilibrium may switches finite times with the increase of τ .
The paper is organized as follows. In section 2, we perform the stability analysis of the positive equilibrium under certain conditions, and calculate the stability switching curves and present the Hopf bifurcation theorem on a two-parameter plane. In section 3, we calculate the normal form at a Hopf singularity along the varying direction (γ 1 µ, γ 2 µ) of (σ, τ ). In section 4, we carry out some numerical simulation to illustrate the theoretical results. Finally, we make a conclusion to sum up our paper.

Stability and existence of Hopf bifurcation
In this section, first we perform the stability analysis of the positive equilibrium, then we will investigate the existence of Hopf bifurcation on the two-parameter plane.
The equilibria of Eq.(2) can be obtained by solving the following equation: Clearly, x 1 = 0 is an extinction equilibrium, x * = ln p δ+H is a positive equilibrium when p > δ + H.
We make the following assumption to guarantee the existence of the positive equilibrium x * : (H 1 ) p > δ + H.
Throughout this paper, we focus on the dynamics analysis near the positive fixed point x * . So we assume that (H 1 ) is always true. The characteristic equation associated with the positive equilibrium x * is given by We can easily draw the following conclusion.
Furthermore, x * is asymptotically stable.
One can see that if ln( p δ + H ) = 1, then the equation (4) is reduced into a equation contain only one term with exponential function. In this situation, the fixed point is x * = 1. So in the following, we investigate the dynamics of (2) in the two cases of p = e(δ + H) and p ̸ = e(δ + H), respectively.
Then we conclude the following conclusions.
Theorem 1 shows that when the birth rate, death rate and the capture rate of the population meet certain relation, the mature delay τ has no influence on the stability of the population; When the capture rate is lower than the death rate, the positive equilibrium is locally stable as long as it exists; When the capture rate is higher than the death rate, the stability of the positive equilibrium is related to the capture delay σ.

Stability analysis when p ̸ = e(δ + H)
In this subsection, we shall perform the dynamics analysis near the positive equilibrium x * when p ̸ = e(δ + H) under the assumption (H 1 ). In this situation, there are two terms with exponent function in (4). In order to analyze the distribution of the roots of (4), the method introduced by An et al. [16] shall be used. We first state the following lemma due to Ruan and Wei [18].

Lemma 2.
As (σ, τ ) varies continuously in R 2 + , the number of zeros (with multiplicity counted) of D(λ; σ, τ ) on C + can change only if a zero appears or cross the imaginary axis.
By direct calculation, we get that the root of Eq.(4) is given by when (τ, σ) = (0, 0), and λ = 0 is not a root of Eq.(4) always. Now we are going to seek purely imaginary roots of (4) to explore the stability switching curves of x * .
We present the following lemma to clarify that in a finite region, the number of stability switching curves is finite.

Lemma 3. The set Ω consists of a finite number of intervals of finite length.
One can verify that τ = τ + k2 when σ = σ + k1 , and τ = τ − k2 when σ = σ − k1 , respectively. Therefore, and the set of all stability switching curves can be written as Proposition 1. T is the set of all stability switching curves on σ − τ plane for Eq. (4). For any (σ, τ ) ∈ T , Eq.(4) has at least a pair of purely imaginary roots ±iω, with ω ∈ Ω.

Hopf bifurcation theorem on two-parameter plane
Now we conclude the following Hopf bifurcation theorem.

Theorem 2.
For any p ∈ T and for any smooth curve S intersecting with T transversely at p, we define the tangent of S at p by ⃗ l. If ∂Reλ ∂ ⃗ l p ̸ = 0, and all other eigenvalues of (4) have non-zero real parts at p, then system (2) undergoes a Hopf bifurcation at p when parameters (σ, τ ) cross T at p along S.
For a better understanding of Theorem 2, we use Fig.2 to illustrate it as follows. If (σ, τ ) cross T at p along the curve S 1 with the tangent vector ⃗ l 1 , and the transversal condition is dReλ d ⃗ l1 | p = 0, then there's no Hopf bifurcation at p along the curve S 1 ; If (σ, τ ) cross T at p along the curve S 2 with the tangent vector ⃗ l 2 , and the transversal condition is dReλ d ⃗ l2 | p ̸ = 0, then system undergoes a Hopf bifurcation at p along the curve S 2 .
For points on stability switching curves in a two parameter plane, the following case may occur. When a point p(σ, τ ) on T vary along a direction ⃗ l 1 , dReλ d ⃗ l1 | p = 0; When this point vary along another direction ⃗ l 2 , dReλ d ⃗ l2 | p ̸ = 0. Therefore, whether p is a Hopf singularity depends on the direction in which (σ, τ ) vary. Now we can summarize the following theorem about the stability of x * .
Substituting σ = σ * + γ 1 µ, τ = τ * + γ 2 µ into (22), then we can rewrite system (22) in the space C as where The linearized system of (24) at the origin is in the following form Then system (24) can be written as where Denote the enlarged space BC of C as formulating (27) in the extended Banach space BC as an abstract ordinary differential equationẋ where A 0 denotes the infinitesimal generators associated with the linear equation , and X 0 is given by X 0 (θ) = 0 for θ ∈ [−1, 0) and X 0 (0) = 1.
Using the Riesz representation theorem, we see that there exists a bounded variation function η(θ) (θ ∈ [−1, 0]) such that In fact, we can choose Let A * denotes the formal adjoint of A 0 on C * := C([0, 1], R * ) as under the bilinear form Then q(θ) = e iω0τ * θ and q * (s) = De −iω0τ * s are the eigenvectors of A 0 and A * corresponding to iω 0 τ * , respectively. And < q * , q >= 1, where We denote Φ = (q(θ), q(θ)). Now we can decompose BC into a center subspace and its orthocomplement, i.e., BC = P ⊕ Kerπ, where π : C → P is the projection defined by For simplicity, we denote According to (31), x t can be composed as where w ∈ BC ∩ Kerπ := Q 1 for any t. Then system (28) is equivalent to the systemż where A 1 is the restriction of A 0 on

Consider the formal Taylor expansion
combining (31) and the Taylor expansion of exponential function, we have We can also get f 2 (0, Φz + w) = f 2 (0, z, w) in (34) as follows where Besides, the term G 3 (0, z, w) can be obtained as follows where

Now Eq.(33) can be written aṡ
Similar to [19,20], use the notations in it, define the operator By a recursive transformations of variables , we can obtain the normal forms of (37). We conclude that this recursive process transforms of (37) into the following equationsż where g j = (g 1 j , g 2 j ), j ≥ 2 have the following form where f j = (f 1 j , f 2 j ) stand for the terms of order j in (z, w), which are obtained after the computation of normal forms up to order j − 1. The normal form truncated to the third order has the following form: Here, g 1 3 (z, 0, 0) = Proj Ker(M 1 3 ) f 1 3 (z, 0, 0) and The calculations of g 1 2 (z, 0, σ) and g 1 3 (z, 0, 0) are in Appendix A. After the calculations, the normal form truncated to the third order on the center manifold for the Hopf singularity is obtained as follows: Make the polar coordinate transformation we obtain the simplified system of (45) as followṡ where a 0 = Re{K 11 }, b 0 = Re{K 21 }.
The dynamics of (46) reflects the dynamics of (45), correspondingly, the dynamics near the Hopf bifurcation of (2) is obtained. Once the varying direction ⃗ l(γ 1 µ, γ 2 µ) of (σ, τ ) is fixed, whether the possible Hopf singularity (σ * , τ * ) is a Hopf singularity will also be determined. Furthermore, the normal form near a Hopf singularity can be calculated, thus the dynamics of the system near a Hopf singularity can be obtained.
Obviously, when µ < 0, system (2) undergoes Hopf bifurcation, the direction of the bifurcation is backward, and the bifurcating periodic solutions are orbitally stable; when µ > 0 and µ is sufficiently small, there's no bifurcating periodic solutions and the positive equilibrium is locally stable. This conclusion is coincide with the simulation of P 3 .
To sum up the above conclusion, we get that when σ is fixed at 0.48, x * is locally asymptotically stable when τ ∈ [0, 1.52883); then x * becomes unstable and system exhibits periodic oscillation when τ ∈ (1.52883, 2.93784); and then x * turns to be locally asymptotically stable again when τ ∈ (2.93784, 3.45538) and finally, x * is unstable for enough large τ . These indicate that for some fixed σ * , the stability of x * switches finite times with the increasing of τ .

Conclusion
In this paper, a scalar model with mature delay and harvest delay is considered. We conclude that when the birth rate, death rate and capture rate of the population meet a certain relation, the mature delay τ has no effect on the stability of the population. Furthermore, when the capture rate is greater than the death rate, the locally stable positive equilibrium would become unstable with the increase of capture delay σ. This indicates that large capture delay can affect people's accurate prediction of population size and lead to over-capture, which may result in the population extinction.
By taking two delays as parameters, the stability switching curves of x * are calculated. Then the two-parameter Hopf bifurcation theorem at a singularity along a direction ⃗ l is presented. Analysis indicates that whether a point is a Hopf singularity depends on the varying direction of two delays.
We emphasize that the varying direction of (σ, τ ) is set to be (γ 1 µ, γ 2 µ) instead of (µ 1 , µ 2 ), thus the calculation of the normal form near the Hopf singularity with two parameters is transformed into one parameter. This dimension reduction method greatly reduces the amount of the calculation. Besides, this method can also be applied to other multi-parameter bifurcation studies.
By aid of the numerical simulations, we verify the theoretical analysis by choosing the varying direction as ⃗ l(γ 1 µ, γ 2 µ) = (0, µ). We obtain that for a fixed σ, the stability of the positive equilibrium may switches finite times with the varying of the τ . This indicates that the sensitivity of τ has great influence on the stability of the system.