Hopf bifurcation analysis in a neutral predator – prey model with age structure in prey

We show a detailed study on the dynamics of a neutral delay differential population model with age structure in the prey species. By selecting the mature delay as a bifurcation parameter, we obtain the stability and Hopf bifurcations of the coexistence equilibrium. Moreover, by computing the normal form on the center manifold, we give the formulas determining the stability of periodic solutions and the direction of Hopf bifurcation. Finally, we give some numerical simulations to support and strengthen the theoretical results.


Introduction
In recent years, neutral-type differential equations have been applied extensively in mathematical biologies and population dynamics [1,3,5], and their dynamical behaviors are usually very complex.In this paper we will consider the neutral delay differential equation with age structure as follows: where V(t) denotes the number of mature individuals in the population at time t, b 1 and b 2 are the birth rate and the maximum possible per capita egg production rate, respectively.µ 0 and µ 1 are the per capita mortality rate of the immature population and the adults, respectively.τ is the generation time and 1/β is the population size at which the whole population reproduces at its maximum rate.One can derive this model (see [3]), from a structured population model for u(t, a) (the population density at time t and age a) as follows, ∂ t u(t, a) + ∂ a u(t, a) = −µ(a)u(t, a), with the death rate µ(a) = µ 1 , a > τ, µ 0 , a < τ.
Corresponding author.Email: gyx@hit.edu.cn Then the total population of the mature individuals is V(t) = ∞ τ u(t, a)da which fulfills Eq. (1.1).For more details on the derivation of (1.1), we refer the readers to [3].
By regarding V(t) as the prey population size, we incorporate a predation term g(V(t))U(t) into (1.1) to obtain the following general neutral prey-predator model: U(t) = −d 1 U(t) + αg(V(t))U(t), V(t) − b 2 e −µ 0 τ V(t − τ) = (b 1 e −βV(t−τ) + b 2 µ 1 )e −µ 0 τ V(t − τ) − µ 1 V(t) − cg(V(t))U(t). (1.2) Here the predator population U(t) consume mature population V(t) as food source, g : R + → R + is the functional response of predators, α represents the conversion rate from prey to predator, c is the per capita capturing rate of prey by a predator per unit time.For the case of linear functional response g( There are many research findings of predator-prey model similar to Eq. (1.3) with b 2 = 0, i.e., in the absence of neutral delay effect, see [6,12,17,18].For those models with neutral delay, Pielou [16] studied the effect of time delay in the population systems.Kuang [13] investigated a logistic model with the neutral term by introducing a time delay and obtained some sufficient conditions to ensure the local asymptotic stability of the positive equilibrium.In addition, there are few results about the delay neutral equation with age structure.Under assumption 0 < b 2 < 1, system (1.3) can be written as an abstract ordinary differential equation in suitable phase space by the theory on the decomposition of the phase space [8,9].Compared with a predator-prey model without neutral delay [11,15,21,22], we find neutral delay may easily induce bifurcations.We obtain the above results by analyzing the characteristic equation and the normal forms near Hopf bifurcations.This paper is organized as follows.In Section 2, we give the stability analysis and bifurcation results about system (1.3).We find Hopf bifurcation occurs when the delay passes through some critical values.In Section 3, we offer algorithms for determining the direction of Hopf bifurcation and the stability of bifurcating periodic solutions by using the normal form method and the center manifold theory.In Section 4, we carry out some numerical simulations to support our results.

Stability of the equilibria and local Hopf bifurcation
In this section, we mainly analyze the stability of positive equilibrium in the cases of τ = 0 and τ > 0. Regarding the stability and bifurcation of delay equation (1.3), the distribution analysis of roots of the characteristic equation plays an important role.Throughout this paper, we assume 0 < b 2 < 1 such that the neutral part of (1.1) defines a stable D-operator (see [8]).The equilibria of system (1.3) are the roots of the following equations: Then system (1.3) has three nonnegative equilibrium points: (i) E 0 = (0, 0), which corresponds to the total extinction of the predator and prey; , which corresponds to the extinction of predator.The boundary equilibrium point E 1 exists if and only if τ < 1 µ 0 ln (iii) The interior equilibrium point E * = (u * , v * ), which corresponds to the coexistence state of prey and predator, and . The dynamics of the predator-free equilibrium E 1 has been investigated in a previous study in [5].Biologically speaking, for model (1.3) we are interested to study the stability behavior of the interior equilibrium point E * .We first transform E * = (u * , v * ) of (1.3) to the origin via the translation u(t) = U(t) − u * and v(t) = V(t) − v * , then the linearization of system (1.3) around the origin takes the form where The characteristic equation of (2.1) is This is a characteristic equation with coefficients depending on time delay, so we employ the method in [2], rewrite (2.2) in the general form where P(λ, τ) Obviously, we have the following lemma.
When τ > 0, the distribution of the roots of the exponential polynomial equation with delay dependent parameters could not be analyzed in the usual way.Therefore, we introduce the geometric criterion established by Beretta and Kuang [2] (see also [14]).Firstly, we need to prove the following statements in order to discuss the distribution of the roots of (2.3) whose coefficients contain the delay τ.Theorem 2.2.For system (2.3), the following geometric stability switch criterions (i)-(v) are established. 2for each τ has at most a finite number of real zeros; (v) Each positive root ω(τ) of F(ω, τ) = 0 is continuous and differentiable in τ whenever it exists.
Proof.Obviously, by Eq. (2.3) we obtain thus the statements (i) and (ii) are established.In fact, therefore (iii) holds true.Notice that F(ω, τ) takes the following expression For each τ, F(ω, τ) = 0 has at most a finite number of real zeros.The statement (v) is also satisfied by the implicit function theorem.
(a) When F(ω, τ) = 0 has no positive roots, there is no stability switch.
(b) When F(ω, τ) = 0 has at least one positive root and each of them is simple, a finite number of stability switches may occur with the increase of τ.
Remark 2.4.Claim (i) implies that λ = 0 is not a characteristic root of (2.3) and (ii) implies that P(λ, τ) and Q(λ, τ) have no identical imaginary roots.Claim (iii) guarantees that the roots of the equation (2.3) with non-negative real parts are uniformly bounded.Claim (iv) guarantees that the equation (2.3) has at most a finite number of imaginary roots for the given τ, that is, there are only finite times for roots to cross the imaginary axis with the change of τ.Claim (v) is used to compute the derivative of the imaginary roots with respect to τ.
Combined with the conclusion of Lemma 2.3, we will study the existence of Hopf bifurcation for system (1.3).Assume that λ = iω (ω > 0) is the pure imaginary root of (2.2), then Separating the real and imaginary parts, we get where Thus, we further make the following assumptions: then there exist two positive real roots ω ± (τ) of (2.5), , then we have ω(τ)τ = θ(τ) + 2nπ.Hence, iω is a purely imaginary root of Eq. (2.2) if and only if τ is a zero of the function S n (τ), defined by We list the basic results to show the occurrence of Hopf bifurcation in [2] as follows.
When τ = τ k , λ = ±iω(τ k ) is a pair of conjugate purely imaginary roots of Eq. (2.3), and if δ(τ k ) > 0(< 0), the roots corresponding to the conjugate purely imaginary roots cross the imaginary axis from left (right) to right (left), where The next theorem is concerned with the stability for system (1.3) and the Hopf bifurcation.Remark 2.7.The stability of the positive equilibrium and the existence of Hopf bifurcation for system (1.3) have been investigated by analyzing the distribution of eigenvalues.It is noteworthy that (2.3) is an exponential polynomial equation with time delay in coefficients.In order to discuss the distribution of the roots of such equations, it is necessary to use the geometric criterion in [2] for the existence of purely imaginary roots, which is shown in Theorem 2.2 above.

Stability and direction of Hopf bifurcation
In this section, we will study the direction of Hopf bifurcation and the stability of bifurcating periodic solutions of model (1.3).The method to compute normal forms restricted on center manifolds is due to [4,7,19,20].For simplicity, let 3) is written as where and for φ = (ϕ 1 , ϕ 2 ) ∈ C, and where
The adjoint operator A * is defined by A * ψ = − dψ ds with domain Applying the formal adjoint theory in [8], we decompose C by Λ as C = P ⊕ Q, where P = span{Φ(θ)} and choose a basis Ψ for the adjoint space P * , such that Ψ, Φ = 1, where •, • is the bilinear form on C * × C defined by Let q ∈ C and q * ∈ C * be the eigenvectors of A(0) and A * corresponding to eigenvalues iω τ and −iω τ, respectively.Then we have where Using a computation process similar to [10], we compute the center manifold C 0 at µ = 0. Define z(t) = q * , u t , w(t, θ) = u t (θ) − 2 Re{z(t)q(θ)}.
On the center manifold C 0 we have w(t, θ) = w(z(t), z(t), θ), where z and z are local coordinates for center manifold C 0 in the direction of q * and q * .Notice that w is real if u t is real.
Theorem 3.1.µ 2 determines the direction of Hopf bifurcation: if µ 2 > 0 (< 0), then the bifurcating periodic solutions are forward (backward); β 2 determines the stability of the bifurcation periodic solutions: the bifurcation periodic solutions are orbitally stable (unstable) on the center manifold if

Numerical simulation
In this section, we shall carry out the numerical simulation on system (1.3).We choose For this set of parameter values, we observe that (H1) and (H2) hold, and τ max = 36.2195.The condition 0 ≤ τ < τ max is used to guarantee the existence of positive equilibrium E * .In particular, τ max > 0 means the birth rate of the prey population is far greater than the death rate of the adults based on the biological mechanism.The conditions (H1) and (H2) hold in the case of τ ∈ [0, τ max ), which provides a possibility for periodic oscillation of system (1.3).According to Theorem 2.6, we know m = 1, the graph of S 0 versus τ shown in Figure 4.1.

Conclusions
In this paper, we study a neutral predator-prey system with age structure in prey.Stability and Hopf bifurcation results about the inner equilibrium are obtained.From the theoretical analysis, one can find that neutral delay terms can alter the dynamics of the Lotka-Volterra system significantly.