Bifurcation analysis in delayed Nicholson blowflies equation with delayed harvest

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. (Appl. Math. Model. 34: 1405–1417, 2010). The stability switching curves by taking two delays as parameters are obtained via the method introduced by An et al.(J. Differ. Equ. 266: 7073–7100, 2019). 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 natural question is how do the delays affect the dynamical behavior of Eq.(2)? The purpose of the paper is to answer this question from the view of bifurcation. At first, we present the Hopf bifurcation curves on a two-parameter plane. We stress that for a possible Hopf singularity, both the varying direction of these two delays and the transversal condition are important. 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 switch finite times with the increase of τ . In this paper, we mainly discuss the case when p is a constant. It should be noted that if p depends on τ , the stability switching curves can be calculated according to the method in [16], thus the corresponding stability analysis of the system can also be obtained.
The paper is organized as follows. In Sect. 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 Sect. 3, we calculate the normal form at a Hopf singularity along the varying direction (γ 1 μ, γ 2 μ) of (σ, τ ). In Sect. 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, we shall perform the stability analysis of the positive equilibrium and investigate the exis-tence of Hopf bifurcation on the two-parameter plane. At first, we would like to mention that the harvest delay σ affects the positivity (or nonnegativity) of solutions of Eq.(2) with positive initial value satisfies x(ϕ, t) > 0 for t ≥ 0. This means that when the harvest delay σ = 0, the solution of Eq.(2) with positive initial value x 0 = ϕ(θ) is positive. Without loss of generality, we assume that 0 < σ ≤ τ . Now we consider the following initial value problem: where ϕ > 0 for θ ∈ [−τ, 0].
We know that max x≥0 xe −x = e −1 . Then it follows that It follows that if H (e δt 0 − 1) , then x(t 0 , ϕ) < 0. This shows that the Claim is true. In the following, we will focus on the stability and bifurcation analysis of this model.
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 * : Throughout this paper, we focus on the dynamical 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. One can see that if ln( p δ + H ) = 1, then the equation (5) is reduced into a equation containing 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.

Stability analysis when p = e(δ + H )
We already know that when p = e(δ + H ), the equilibrium x * = 1, and Eq. (5) becomes 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 maturation delay τ has no influence on the local 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 local stability of the positive equilibrium is related to the capture delay σ .

Stability analysis when p = e(δ + H )
In this subsection, we shall perform the dynamical 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 (5). In order to analyze the distribution of the roots of (5), 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.

The case of p being a constant
By direct calculation, we get that the root of Eq. (5) is given by . Now we are going to seek purely imaginary roots of (5) to explore the stability switching curves of x * . Let λ = iω (ω > 0) be a root of Eq. (5), and substitute it into Eq. (5), we have For simplicity, denote Thus, we have which equals the existence of nonnegative σ satisfying (11) if and only if Define and denote the set of all positive ω's satisfying (12) as 1 . Let Then we have Similarly, we get that where and ψ ∈ (−π, 2π ] ∩ Arg{δ + iω} such that Here, the condition on ω is as follows: Denote the set of all positive ω's satisfying Eq. (17) as 2 . We can easily prove that (12) is equivalent to (17) by squaring both sides of the two inequalities. Thus, 1 = 2 : . 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 τ = τ + k 2 when σ = σ + k 1 , and τ = τ − k 2 when σ = σ − k 1 , 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. (5). For any (σ, τ ) ∈ T , Eq. (5) has at least a pair of purely imaginary roots ±iω, with ω ∈ .
If (29) has zeros for some ω ∈ I k , say τ i± (ω), i = 1, 2, ..., lying in I k ω , then we can get the corresponding σ value: will determine the values of τ and σ , for which Eq.(5) has a pair of purely imaginary roots ±iω. Therefore, the set of all such points defines the stability switching curves

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 (5) 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
Substituting σ = σ * + γ 1 μ, τ = τ * + γ 2 μ into (34), then we can rewrite system (34) in the space C as where The linearized system of (36) at the origin is in the following form Then system (36) can be written as where Denote the enlarged space BC of C as formulating (39) in the extended Banach space BC as an abstract ordinary differential equatioṅ where A 0 denotes the infinitesimal generators associated with the linear Eq. (38), defined by , and X 0 is given by X 0 (θ ) = 0 for θ ∈ [−1, 0) and X 0 (0) = 1.
For simplicity, we denote According to (43), x t can be composed as where w ∈ BC ∩ Kerπ := Q 1 for any t. Then system (40) is equivalent to the systeṁ where A 1 is the restriction of A 0 on Consider the formal Taylor expansion Notice that combining (43) and the Taylor expansion of exponential function, we have where Besides, the term G 3 (0, z, w) can be obtained as follows where Now Eq. (45) can be written aṡ where w ∈ Q 1 , B = diag{iωτ * , −iωτ * } and f j = ( f 1 j , f 2 j ) T , j ≥ 2, are defined by Similarly to [19,20], use the notations in it and define the operator By a recursive transformations of variables we can obtain the normal forms of (49). We conclude that this recursive process transforms of (49) into the following equationṡ where g j = (g 1 j , g 2 j ), j ≥ 2 have the following form stand for the terms of order j in (z, w), which are obtained after the computation of normal forms up to order j-1. Since the first three terms in the normal form determine the properties of the system, the higher-order terms won't need to be calculated. 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 (57) as followṡ where The dynamics of (58) reflects the dynamics of (57), 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.

Numerical simulations
In this section, we carry out some simulations to illustrate the above theoretical results.
Obviously, when μ > 0, system (2) undergoes Hopf bifurcation, the direction of the bifurcation is forward, and the bifurcating periodic solutions are orbitally stable. We choose μ = 0.00062, (σ, τ ) becomes P 4 (0.48, 3.456), we illustrate the dynamics of system (2) at P 4 in Fig. 7d, where the initial value is For Q 2 , we can draw ReK 11 < 0, ReK 21 < 0 from Table 1. The simplified system of (57) is in the following forṁ 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 coincided 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 maturation 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 maturation 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 overcapture, which may result in the population extinction.
By taking two delays as parameters, the stability switching curves of x * are calculated. Then the twoparameter 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 switch finite times with the varying of the τ . This indicates that the sensitivity of τ has great influence on the stability of the system. with K 21 = L 21 + 3 2 (M 21 + N 21 ) .