Electronic Journal of Differential Equations, Vol. 2021 (2021), No. 42, pp. 1–16. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu STABILITY AND BIFURCATION IN A DELAYED PREDATOR-PREY MODEL WITH HOLLING-TYPE IV RESPONSE FUNCTION AND AGE STRUCTURE YUTING CAI, CHUNCHENG WANG, DEJUN FAN Abstract. In this article, we study a predator-prey model with age structure, Holling-type IV response, and two time delays. By an algebraic method, we determine all the critical values for these two delays, such that the charac- teristic equation has purely imaginary roots. This provides a sharp stability region on the parameter plane of the positive equilibrium. Applying integrated semigroup theory and Hopf bifurcation theorem for abstract Cauchy problems with non-dense domain, we can show the occurrence of Hopf bifurcation as the time delays pass through these critical values. In particular, the phenomenon of stability switches can also be observed as the time delays vary. Numerical simulations are carried out to illustrate the theoretical results. 1. Introduction There is a long history of studies on the interaction between predator and prey in ecology [4, 10, 22, 25]. A typical mathematical model for describing this interaction is ẋ = p(x)− yq(x), ẏ = cyq(x)− dy, (1.1) where x(t) and y(t) are the population densities of the prey and predator, respec- tively; p(x) is the birth function of prey; q(x) is the intake rate of predator as a function of prey; c is the assimilation efficiency of predator; and d represents preda- tor’s death rate. Equation (1.1) has been extensively studied for various choices of p(x) and q(x). The logistic growth is frequently assumed for the increment of prey, that is, p(x) = rx(1− x/K), where r denotes intrinsic growth rate of prey, and K is the carrying capacity of the prey. Holling type functional responses are probably the most commonly used functions for q(x), and various dynamics, such as global attractivity of equilibrium and existence of limit cycle, can be observed for (1.1) with those functional responses, see [8, 9, 12]. To study the phenomenon of group defense, non-monotonic functional responses are also proposed. A typical choice of non-monotonic functional response is the so-called Holling type IV functional response, i.e. q(x) = mx x2+b1x+a , and we refer the readers to [5, 14] for more details. 2010 Mathematics Subject Classification. 37G10, 35F31. Key words and phrases. Age-structured model; Hopf bifurcation; Holling-type IV response. c©2021 Texas State University. Submitted January 15, 2020. Published May 14, 2021. 1 2 Y. CAI, C. WANG, D. FAN EJDE-2021/42 If a time delay τ is incorporated in the model (1.1), then it turns into dx(t) dt = p(x)− yq(x(t− τ)), dy(t) dt = cyq(x(t− τ))− dy. (1.2) Here, the delay τ takes into account the transformation time from prey quantities into predator populations. The dynamical behaviors of (1.2) are shown to be much more complicated than (1.1), see [6, 11, 21, 24] and references therein. In particular, the authors also consider (1.2) for logistic growth p(x) and Holling type IV q(x) with b1 = 0 [18, 28]. However, not much attention has been paid in the case of b1 6= 0 [?]. In this article, we consider a more generalized version of (1.2), namely ∂u(t, a) ∂t + ∂u(t, a) ∂a = −µu(t, a), dV (t) dt = rV (t) ( 1− V (t) K ) − V (t− τ1) ∫ +∞ 0 β(a)u(t, a)da V 2(t− τ1) + b1V (t− τ1) + 1 , u(t, 0) = V (t− τ1) ∫ +∞ 0 β(a)u(t, a)da V 2(t− τ1) + b1V (t− τ1) + 1 , u(0, ·) = u0 ∈ L1((0,+∞),R), V0 = φ ∈ C([−τ1, 0],R), (1.3) where u(t, a) is the density of predator of age a at time t; V (t) is the population of prey; β(a) is the maturation function which describes the effects of age on fecundity. If β(a) = 1[0,+∞)(a) and U(t) = ∫ +∞ 0 u(t, a)da, then (1.3) will degenerate to (1.2). Throughout this paper, the kernel β(a) is assumed to take the form β(a) := { β∗, a ≥ τ2, 0, a ∈ [0, τ2), such that ∫ +∞ 0 β(a)e−µada = M < +∞, where τ2 > 0 is the maturation period of predator. For age-structure model like (1.3), it is conjectured in [2] that the maturation period τ2 could induce periodic oscillation. However, this is not rig- orously proved until the development of Hopf bifurcation theory for an abstract non-densely defined Cauchy problem in [16]. This theory is established on the basis of centre manifold theorem [19] and normal form reduction [17], under the frame- work of integrated semigroups [26], and has been successfully applied to many age structured models [27, 15, 29]. Recall that the delay τ may also destablize the posi- tive equilibrium of (1.2), generating periodic solutions. Therefore, it is our interest to investigate the interactive impact of τ1 and τ2 on the dynamics of (1.3). By ana- lyzing the characteristic equation of (1.3) at the positive equilibrium, we determine all the values of (τ1, τ2) such that the characteristic equation has roots with zero real parts. This will give the sharp region on (τ1, τ2)-plane, where (1.3) has locally stable positive equilibrium. Furthermore, as (τ1, τ2) passes through the boundary of this region, we can show the existence of periodic solution with different period by Hopf bifurcation theorem. It should be mentioned that the characteristic equation of (1.3) takes the form P0(λ) + P1(λ)e−λτ1 + P2(λ)e−λτ2 + P3(λ)e−λ(τ1+τ2) = 0, (1.4) EJDE-2021/42 AN AGE-STRUCTURED PREDATOR-PREY MODEL 3 where Pk(λ), k = 1, 2, 3 are polynomials of λ. In [13], the authors present a sys- tematic analytic method for finding the crossing curves (on which (1.4) has purely imaginary roots) of (1.4) on (τ1, τ2)-plane, which will be employed here to study the characteristic equation of (1.3). When P3(λ) = 0, the authors in [7] propose a geometric method for detecting the crossing curves, and list all the possibility of their shapes. As a special case, equation (1.4) with τ1 = τ2 is also studied in [1], presenting the explicit formula for the critical values of Hopf bifurcation. This article is organized as follows: In Section 2, we formulate the initial age- structured model (1.3) as an abstract Cauchy problem. In Section 3, we show the existence of the positive steady state of system (1.3), and then derive its the characteristic equation. In Section 4, we analyze the distribution of the roots for characteristic equation, and obtain the the crossing curves on (τ1, τ2)-plane and crossing directions. These curves will not only determine the stability region for the positive steady state, but also tell where Hopf bifurcation could take place. In particular, for τ1 = τ2 = τ , we derive the explicit formula for Hopf bifurcation values of τ . In Section 5, we conduct some numerical simulations to verify the results. 2. Transformation to Cauchy problem Let â = a τ2 , t̂ = t τ2 and V̂ (t̂) = V (τ2t̂), û(t̂, â) = τ2u(τ2t̂, τ2â). Then, system (1.3) becomes, after dropping the hat, ∂u(t, a) ∂t + ∂u(t, a) ∂a = −µτ2u(t, a), dV (t) dt = τ2 [ rV (t) ( 1− V (t) K ) − V (t− τ1 τ2 ) ∫ +∞ 0 β(a)u(t, a) da V 2(t− τ1 τ2 ) + b1V (t− τ1 τ2 ) + 1 ] , u(t, 0) = τ2V (t− τ1 τ2 ) ∫ +∞ 0 β(a)u(t, a) da V 2(t− τ1 τ2 ) + b1V (t− τ1 τ2 ) + 1 , u(0, ·) = u0 ∈ L1((0,+∞),R), V0 = φ ∈ C([−τ1 τ2 , 0],R), (2.1) where the function β(a) is now defined by β(a) = { β∗, a ≥ 1, 0, 0 ≤ a < 1, with β∗ = δµMeδµτ2 . To derive the positive equilibrium and its characteristic equation, we need to formulate it into as an abstract non-densely defined Cauchy problem. This can be accomplished by the following steps: Step 1: Rewrite (2.1) as a system of first order partial differential equations. We denote by ρ(t, a) the density of prey of age a at time t. Let V (t) := ∫ +∞ 0 ρ(t, a) da. Then, from the second equation in system (2.1), we have ∂ρ(t, a) ∂t + ∂ρ(t, a) ∂a = −τ2dρ(t, a), ρ(t, 0) = G(u(t, a), ρ(t, a)), 4 Y. CAI, C. WANG, D. FAN EJDE-2021/42 ρ(0, a) = ρ0 ∈ L1((0,+∞),R). Here G(u(t, a), ρ(t, a)) =τ2 [ b ∫ +∞ 0 ρ(t, a) da ( 1− ∫ +∞ 0 ρ(t, a) da K ) + d( ∫ +∞ 0 ρ(t, a) da)2 K − ∫ +∞ 0 ρ(t− τ1 τ2 , a) da ∫ +∞ 0 β(a)u(t, a) da [ ∫ +∞ 0 ρ(t− τ1 τ2 , a) da]2 + b1 ∫ +∞ 0 ρ(t− τ1 τ2 , a) da+ 1 ] , and b, d represent the birth and death rate of prey respectively, such that r = b−d. Set w(t, a) = ( u(t, a) ρ(t, a) ) . We can rewrite (2.1) as ∂w(t, a) ∂t + ∂w(t, a) ∂a = −Dw(t, a), w(t, 0) =  τ2 ∫ +∞ 0 ρ(t− τ1τ2 ,a) da ∫ +∞ 0 β(a)u(t,a) da ( ∫ +∞ 0 ρ(t− τ1τ2 ,a) da) 2+b1 ∫ +∞ 0 ρ(t− τ1τ2 ,a) da+1 G (u(t, a), ρ(t, a))  , w(0, a) = w0 = ( u0 ρ0 ) ∈ L1((0,+∞),R2), (2.2) where D = ( τ2µ 0 0 τ2d ) . Step 2: Write (2.2) as an abstract delay equations. Consider the Banach space X := R2 × L1 ( (0,+∞),R2 ) , with the usual product norm ∥∥(ζ ϕ )∥∥ X = ‖ζ‖R2 + ‖ϕ‖L1 , for ( ζ ϕ ) ∈ X, and the space CA = {(ζ(·) φ(·) ) ∈ C ([ − τ1 τ2 , 0 ] , X ) : ζ(0) = 0 } . Define the linear operator L : D(L) ⊂ X → X by L ( 0 ϕ ) = ( −ϕ(0) −ϕ′ −Dϕ ) with D(L) = {0} ×W 1,1 ( (0,+∞),R2 ) , and the operator F : CA → X by F (( ζ(·) φ(·) )) = ( B (φ(·)) 0L1 ) , where B(φ(·)) =  τ2 ∫ +∞ 0 φ2(− τ1τ2 )(a) da ∫ +∞ 0 β(a)φ1(0)(a) da ( ∫ +∞ 0 φ2(− τ1τ2 )(a) da)2+b1 ∫ +∞ 0 φ2(− τ1τ2 )(a) da+1 G (φ1(0)(·), φ2(0)(·))  . Then, by setting y(t) = ( 0 w(t, a) ) , we can rewrite (2.2) as the Cauchy problem dy(t) dt = Ly(t) + F (yt), t ≥ 0, y(0) = ( 0 w0 ) ∈ CA. (2.3) EJDE-2021/42 AN AGE-STRUCTURED PREDATOR-PREY MODEL 5 Note that L is non-densely defined, since X0 := D(L) = {0} × L1 ( (0,+∞),R2 ) . Step 3: Transform (2.3) into an abstract ordinary differential equation. We define x ∈ C ( [0,+∞) × [− τ1τ2 , 0];X ) by x(t, θ) = y(t + θ), for t ≥ 0 and θ ∈ [− τ1τ2 , 0]. Therefore, x(t, θ) satisfies the equation ∂x(t, θ) ∂t − ∂x(t, θ) ∂θ = 0, θ ∈ [−τ1 τ2 , 0), ∂x(t, 0) ∂θ = Lx(t, 0) + F (x(t, ·)), θ = 0, x(0, ·) = y0 ∈ CA. (2.4) Let Z = X × C, C := C([− τ1τ2 , 0], X) with ( f φ ) = ‖f‖X + ‖φ‖C . If we set z(t) :=( 0 x(t) ) , then z satisfies d dt z(t) = Az(t) +H(z(t)), t > 0, z(0) = ( 0X y0 ) ∈ Z0, (2.5) where A : D(A) ⊂ Z → Z and H : Z0 → Z with Z0 = {0X} × CA are given by A ( 0X φ ) = ( −φ′(0) + Lφ(0) φ′ ) , H ( 0X φ ) = ( F (φ) 0CA ) . It is straightforward to show that (2.5) is an abstract non-densely defined Cauchy problem, because D(A) = {0X}×{φ ∈ C1([− τ1τ2 , 0], X), φ(0) ∈ D(L)} and D(A) = Z0 6= Z. The existence and uniqueness of a global and positive solution for system (2.5) follow directly from the results in [20]. 3. Equilibria and characteristic equation 3.1. Equilibria. Suppose that z = ( 0X ψ ) ∈ D(A) is an equilibrium of (2.1), where ψ = ( ζ(·) φ(·) ) ∈ C1 ( [−τ1 τ2 , 0];X ) , ψ(0) ∈ D(L), φ(·) = ( φ1(·) φ2(·) ) . Then, solving the system −ψ′(0) + Lψ(0) + F (ψ) = 0, ψ ′ = 0, (3.1) we arrive at the following conclusion on the existence of a positive equilibrium. 6 Y. CAI, C. WANG, D. FAN EJDE-2021/42 Lemma 3.1. System (2.5) always has the boundary equilibria z1 =  0X ζ1(·)( φ11(·) φ12(·) )  and z2 =  0X ζ2(·)( φ21(·) φ22(·) )  with ζ1(θ) = ζ2(θ) = 0R2 ,( φ11(θ)(a) φ12(θ)(a) ) = ( 0 0 ) , ( φ21(θ)(a) φ22(θ)(a) ) = ( 0 τ2dKe −τ2da ) . Furthermore, if (H1) K > 1 and b−M = −2 hold, then there exists a unique positive equilibrium of system (2.5), given by z =  0X ζ(·)( φ1(·) φ2(·) )  , ζ(θ) = 0R2 , ( φ1(θ)(a) φ2(θ)(a) ) = ( C1e −τ2µa C2e −τ2da ) , where C1 = τ2µMr(K − 1) K , C2 = τ2d. The linearized equation of (2.5) around the equilibrium z̄ is dz(t) dt = Az(t) +DH(z)z(t), where DH(z) ( 0X ψ ) = ( DF (ψ)(ψ) 0CA ) , ∀ ( 0X ψ ) ∈ D(A), ψ = ( ζ(·) φ(·) ) with DF (ψ)(ψ) = ( DB(φ)(φ) 0L1 ) and DB(φ)(φ) = ( τ2 2 dC2 C2 2+b1dτ2C2+τ2 2 d 2 0 0 0 )∫ +∞ 0 β(a)φ(0)(a) da + 0 C1Mτ2 2 d 2(τ2 2 d 2−C2 2 ) (C2 2+b1dτ2C2+τ2 2 d 2)2 0 − C1τ 2 2 d 2(τ2 2 d 2−C2 2 ) µ(C2 2+b1dτ2C2+τ2 2 d 2)2 ∫ +∞ 0 φ ( −τ1 τ2 ) (a) da + ( 0 0 − C2τ 2 2 d C2 2+b1dτ2C2+τ2 2 d 2 τ2bdK−2rC2 dK )∫ +∞ 0 φ(0)(a) da. 3.2. Characteristic equation. Now, we derive the characteristic equation, that determines the distribution of spectrum of A+DH(z), to study the local dynamical behavior of z̄. Denote ϑ := min{τ2d, τ2µ}, and Ω := {λ ∈ C : Re(λ) > −ϑ}. The next lemma from [3] will be used. Lemma 3.2. The operators L and A defined above satisfy the following statements: EJDE-2021/42 AN AGE-STRUCTURED PREDATOR-PREY MODEL 7 (i) If λ ∈ Ω, then λ ∈ ρ(L), and (λI − L)−1 ( δ ψ ) = ( 0 ϕ ) ⇔ ϕ(a) = e− ∫ a 0 (λI+D)dlδ + ∫ a 0 e− ∫ a s (λI+D)dlψ(s)ds with ( δ ψ ) ∈ X and ( 0 ϕ ) ∈ D(L). (ii) ρ(L) = ρ(A). Moreover, for each λ ∈ ρ(A), we also have the following explicit formula for the resolvent of A, (λI −A)−1 ( f φ ) = ( 0X φ ) ⇔ φ(θ) = eλθ(λI − L)−1[ψ(0) + f ] + ∫ 0 θ eλ(θ−s)ψ(s) ds. Denote G := DH(z). By Lemma 3.2, we know (λI − A) is invertible for any λ ∈ Ω. It follows from the identity [I − G(λI −A)−1] = (λI −A)−1[λI − (A+ G)], that λI − (A + G) is invertible if and only if I − G(λI − A)−1 is invertible. Now, we consider [I − G(λI −A)−1] ( δX ϕCA ) = ( γX ψCA ) , (3.2) where δX = ( δ1 δ2 ) , γX = ( γ1 γ2 ) , ψCA = ( ψ1(·) ψ2(·) ) ∈ C1 ( [−τ1 τ2 , 0]; X ) . Equation (3.2) is equivalent to the equation δX −DF (ψ) [ eλθ(λI − L)−1 (ϕCA(0) + δX) + ∫ 0 θ eλ(θ−s)ϕCA(s) ds ] = γX , ϕCA = ψCA . (3.3) From the first equation of (3.3), we obtain( I −DB(φ) ( eλθe− ∫ a 0 (λI+D) dl )) δ1 = δ1 +DB(φ) [ eλθ ∫ a 0 e− ∫ a s (λI+D) dlδ2(s) ds ] , δ2 = δ2 , which has a unique solution δ1 for any right hand side term if and only if ∆(λ) is invertible, where ∆(λ) = I −DB(φ) ( eλθe− ∫ a 0 (λI+D) dl ) . (3.4) Theorem 3.3. For the spectrum of A+ G, we have σ(A+ G) ∩ Ω = σp(A+ G) ∩ Ω = {λ ∈ Ω : det(∆(λ)) = 0}. Proof. Let λ ∈ Ω and det(∆(λ)) 6= 0. It then follows from (3.3) that (I − DF (ψ)(eλθ(λI − L)−1)) is invertible, and (I −DF (ψ)(eλθ(λI − L)−1))−1 ( δ1 δ2 ) = δX , where δX = ( (∆(λ))−1(δ1) +DB(φ) [ eλθ ∫ a 0 e− ∫ a s (λI+D) dlδ2(s) ds ] δ2 ) . 8 Y. CAI, C. WANG, D. FAN EJDE-2021/42 Thus, (I − G(λI −A)−1) is invertible, and (I − G(λI −A)−1)−1 ( γX ψCA ) = ( δX ψCA ) . Hence, we have {λ ∈ Ω : det(∆(λ)) 6= 0} ⊂ ρ(A+G)∩Ω, which implies σ(A+G)∩ Ω ⊂ {λ ∈ Ω : det(∆(λ)) = 0}. Suppose that λ ∈ Ω and det(∆(λ)) = 0. We are going to find ( 0 ψ ) ∈ D(A)\{0} such that (A+ G) ( 0 ψ ) = λ ( 0 ψ ) . From the above argument, it suffices to find ( α ϕ ) ∈ Z\{0} satisfying [I − G(λI − A)−1] ( α ϕ ) = 0, or equivalently to find ( α ϕ ) 6= 0 such that ∆(λ)α = 0, ϕ = 0. Since det(∆(λ)) = 0, there always exists α 6= 0 such that ∆(λ)α = 0. Thus, λ ∈ σP (A+ G). This proves {λ ∈ Ω : det(∆(λ)) = 0} ⊂ σP (A+ G). � After some symbolic manipulation, we obtain the characteristic equation of (2.1), det(∆(λ)) = λ2 + τ2Pλ+ τ22Q+ (τ2Sλ+ τ22R)e−λ + (τ2Nλ+ τ22Y )e− τ1 τ2 λ + τ22He − τ1τ2 λ · e−λ (λ+ τ2d)(λ+ τ2µ) =: f1(λ) f2(λ) = 0, (3.5) where P = µ− r + 2rF K , Q = µ(−r + 2rF K ), S = −FMµ G , R = rFMµ( K − 2F GK ), Y = E(1− F 2) G2 + EFM(1− F 2) G3 , N = E(1− F 2) µG , H = −EFM(1− F 2) G3 , E = C1 τ2 , F = C2 τ2d , G = C2 2 + bdC2τ2 + τ22 d 2 τ22 d 2 . Let λ = τ2ζ. Then f1(λ) = f1(τ2ζ) = τ22 g(ζ), where g(ζ) = ζ2 + Pζ +Q+ (Sζ +R)e−τ2ζ + (Nζ + Y )e−τ1ζ +He−(τ1+τ2)ζ . (3.6) EJDE-2021/42 AN AGE-STRUCTURED PREDATOR-PREY MODEL 9 4. Hopf bifurcation In this section, we analyze the distribution of purely imaginary roots of g(ζ) = 0. Specifically, the crossing curves of g(ζ) = 0 will be determined in the case of τ1 6= τ2, and as a special case of τ1 = τ2 = τ , we will derive the explicit formula for the critical Hopf bifurcation values. 4.1. The case τ1 6= τ2. We use the method in [13] to study g(ζ) = 0. Substituting ζ = iω, ω > 0 into g(ζ) = 0, we have (−ω2 + Piω +Q+ (Niω + Y )e−iωτ1) + (Siω +R+He−iωτ1)e−iωτ2 = 0. (4.1) It then follows from |e−iωτ2 | = 1 that | − ω2 + Piω +Q+ (Niω + Y )e−iωτ1 | = |Siω +R+He−iωτ1 |, which implies ω4 − (2Q− P 2 −N2 + S2)ω2 + (Q2 + Y 2 −R2 −H2) = 2A1(ω) cos(ωτ1)− 2B1(ω) sin(ωτ1), (4.2) where A1(ω) = (Y − PN)ω2 + (RH − Y Q), B1(ω) = (SH − PY +QN)ω −Nω3. Let φ1(ω) = arg{(Y − PN)ω2 +RH − Y Q+ (−Nω2 + SH − PY +QN)iω}. Then A1(ω) and A2(ω) can be written as A1(ω) = √ ((Y − PN)ω2 + (RH − Y Q))2 + ((SH − PY +QN)ω −Nω3)2 × cos(φ1(ω)), B1(ω) = √ ((Y − PN)ω2 + (RH − Y Q))2 + ((SH − PY +QN)ω −Nω3)2 × sin(φ1(ω)). (4.3) Substituting (4.3) into (4.2), we obtain ω4 − (2Q− P 2 −N2 + S2)ω2 + (Q2 + Y 2 −R2 −H2) = 2 √ ((Y − PN)ω2 + (RH − Y Q))2 + ((SH − PY +QN)ω −Nω3)2 × cos(φ1(ω) + ωτ1). (4.4) For (4.4) to have positive root ω, it is necessary that F (ω) = [ω4 − (2Q− P 2 −N2 + S2)ω2 + (Q2 + Y 2 −R2 −H2)]2 − 4[((Y − PN)ω2 + (RH − Y Q))2 + ((SH − PY +QN)ω −Nω3)2] < 0. (4.5) On the other hand, if ω > 0 satisfies (4.5), then we can always find τ1 such that (ω, τ1) is the root of (4.4). The set of all possible values of ω > 0 satisfying (4.5) is denoted by Ω. From [13, Lemma 3.2], it follows that Ω consists of a finite number of intervals of finite length Ωk; that is, Ω = ⋃N k=1 Ωk with |Ωk| <∞. Let ψ1(ω) = φ1(ω) + ωτ1 ∈ [0, π]. 10 Y. CAI, C. WANG, D. FAN EJDE-2021/42 Then cos(ψ1(ω)) = ω4 − (2Q− P 2 −N2 + S2)ω2 + (Q2 + Y 2 −R2 −H2) 2 √ ((Y − PN)ω2 + (RH − Y Q))2 + ((SH − PY +QN)ω −Nω3)2 , and therefore, τ±1,n1 (ω) = ±ψ1(ω)− φ1(ω) + 2n1π ω , n1 ∈ Z. (4.6) By an argument as above, if we define A2(ω = (Y H −RQ) + (R− PS)ω2 = √ A2(ω)2 +B2(ω)2 cos(φ2(ω)), B2(ω) = (NH − PR)ω + S(Q− ω2)ω = √ A2(ω)2 +B2(ω)2 sin(φ2(ω)), for φ2(ω) = arg{A2(ω) + iB2(ω)}, and ψ2(ω) = φ2(ω) + ωτ2 ∈ [0, π], then τ±2,n2 (ω) = ±ψ2(ω)− φ2(ω) + 2n2π ω , n2 ∈ Z. (4.7) According to [13, p.525], the crossing curves of g(ζ) = 0 are given in the following theorem. Theorem 4.1. The crossing curves of g(ζ) = 0 are given by Γ = ( ∪k=1,2,...,N Γ±kn1,n2 ) ∩ R2 +, (4.8) where Γ±kn1,n2 = { ( ±ψ1(ω)− φ1(ω) + 2n1π ω , ∓ψ2(ω)− φ2(ω) + 2n2π ω ) : ω ∈ Ωk } . (4.9) Note that these crossing curves may intersect τ1 or τ2 axis. One can easily identify these intersections along the following lines. For example, if τ1 = 0 and τ2 > 0, then g(ζ) = 0 turns into ζ2 + (P +N)ζ +Q+ Y + (Sζ +R+H)e−τ2ζ = 0. (4.10) Assume that (H2) P + S +N > 0, Q+H +R+ Y > 0, and Q+H −R− Y < 0. Then, all the roots of (4.10) with τ2 = 0 have strictly negative real parts. Let iω20, ω20 > 0 be the root of (4.10). Then, we can obtain −ω2 20 + (P +N)iω20 + Y +Q+ (Siω20 +R+H)e−iτ2ω20 = 0. Separating real and imaginary parts, we obtain ω2 20 + Y +Q = −(R+H) cosω20τ2 − Sω20 sinω20τ2, (P +N)ω20 = (R+H) sinω20τ2 − Sω20 cosω20τ2, from which it follows that f2(ω20) := ω4 20 +[(P +N)2−S2−2(Y +Q)]ω2 20 +(Y +Q)2− (R+H)2 = 0. (4.11) Equation (4.11) will have a unique positive root ω∗20, as long as (H3) Y +Q−R−H < 0 EJDE-2021/42 AN AGE-STRUCTURED PREDATOR-PREY MODEL 11 holds. Moreover, the associated critical values for τ2 are τ2n = 1 ω∗20 arccos[ (R+H)(ω∗220 − Y −Q) + S(P +N)ω∗220 (R+H)2 + S2ω∗220 ] + 2nπ ω∗20 , (4.12) for n = 0, 1, 2 . . . . Applying the implicit function theorem to (4.10), one has ( dζ dτ2 )−1 = −τ2 ζ + S ζ(R+H + Sζ) − 2ζ + (P +N) ζ(ζ2 + (P +N)ζ + Y +Q) . Since sign{d(Re ζ) dτ2 }−1τ2=τ2n = sign{Re S ζ(R+H + Sζ) +Re 2ζ + (P +N) ζ(ζ2 + (P +N)ζ + Y +Q) }τ2=τ2n = sign[ 2ω2 20 + S2 + (P +N)2 − 2(Y +Q) (R+H)2 + S2ω2 20 ]. To make sign{d(Reζ)dτ2 }τ2=τ2n > 0, we assume that (H4) S2 + (P +N)2 − 2(Y +Q) > 0, so that we have the following conclusions. Corollary 4.2. Assume that τ1 = 0 and (H1)–(H4) hold. Then the interior equi- librium z̄ is locally asymptotically stable for 0 < τ2 < τ20, and (2.1) undergoes a Hopf bifurcation at τ2 = τ20, where τ20 is defined by (4.12). Along the same lines as above, we have the following bifurcation results for the case τ2 = 0 and τ1 > 0. Corollary 4.3. Assume that τ2 = 0 and (H1), (H2), (H5) hold. Then, the interior equilibrium z̄ is locally asymptotically stable for 0 < τ1 < τ10, and (2.1) undergoes a Hopf bifurcation at τ1 = τ10, where τ10 = 1 ω∗10 arccos[ (Y +H)(ω∗210 − Y −Q) +N(P + S)ω∗210 (Y +H)2 +N2ω∗210 ], and ω∗10 is the unique positive root of f1(ω10) := ω4 10 + [(P +S)2−N2− 2(R+Q)]ω2 20 + (R+Q)2− (Y +H)2 = 0 (4.13) if and only if (H5) R+Q− Y −H < 0. Next, we discuss the direction in which the roots of g(ζ) = 0 cross the imaginary axis, as (τ1, τ2) deviates from the curve in Γ. As in [7], we call the direction of the curve corresponding to increasing ω ∈ Ωk the positive direction, and the region on the left-hand (right-hand) side as we head in the positive direction of the curve the region on the left (right). As a direct consequence of [7, Proposition 6.1], we have the following conclusions. Theorem 4.4. Let ω ∈ Ωk and (τ1, τ2) ∈ Γ±kn1,n2 such that iω is a simple solution of g(ζ) = 0. Then, as (τ1, τ2) moves from the region on the right to the region on the left of the crossing curve, a pair of complex roots of g(ζ) = 0 cross the imaginary axis to the right if R2I1 −R1I2 > 0, (4.14) 12 Y. CAI, C. WANG, D. FAN EJDE-2021/42 where R1 = Re{∂g(ζ, τ1, τ2) ∂τ1 } = Nω2 cosωτ1 − Y ω sinωτ1 − ωH sinω(τ1 + τ2), I1 = Im{∂g(ζ, τ1, τ2) ∂τ1 } = −Nω2 sinωτ1 − Y ω cosωτ1 − ωH cosω(τ1 + τ2), R2 = Re{∂g(ζ, τ1, τ2) ∂τ2 } = Sω2 cosωτ2 −Rω sinωτ2 − ωH sinω(τ1 + τ2), I2 = Im{∂g(ζ, τ1, τ2) ∂τ2 } = −Sω2 sinωτ2 −Rω cosωτ2 − ωH cosω(τ1 + τ2). The crossing is in the opposite direction if the (4.14) is reversed. 4.2. The case τ1 = τ2. When τ1 = τ2, g(ζ) = 0 becomes h(ζ) = ζ2 + Pζ +Q+ ((S +N)ζ + (R+ Y ))e−τζ +He−2τζ . (4.15) Suppose that iω(ω > 0) is a root of (4.15), that is, − ω2 + Piω +Q+ (S1iω +R1)e−τiω +He−2τiω = 0, (4.16) where R1 = R + Y and S1 = S +N . Applying equation (2.5) of page 6 in [1], one can conclude that ω must satisfy the polynomial equation ω8 + s1ω 6 + s2ω 4 + s3ω 2 + s4 = 0, (4.17) which has at least one positive solution for ω, by the third inequality of (H2). Here, s1 = 2P 2 − 4Q− S2 1 , s2 = 6Q2 − 2H2 − 4QP 2 −R2 1 + P 4 − P 2S2 1 + 2S2 1Q+ 2HS2 1 , s3 =2R2 1Q− P 2R2 1 − 4Q3 + 2Q2P 2 − S2 1Q 2 − 2QS2 1H + 4PS1R1H − 2R2 1H + 4QH2 − 2H2P 2 − S2 1H 2, s4 = [ −R2 1 + (Q+H)2 ] (Q−H)2. To present the explicit formula of the critical bifurcation values for τ as in [1], we also define θ = F (ω) D(ω) , (4.18) where D(ω) = 2 [ −ω4 + (2Q−R1 + P (−P + S1))ω2 + (−Q+R1 −H)(Q−H) ] , F (ω) = 2ω [ −R1P − S1ω 2 + (Q+H)S ] . Theorem 4.5. Assume that τ1 = τ2 = τ , (H1) and (H2) holds. Then, equation (4.17) has a finite number of positive roots, denoted by ωN , and (1.3) undergoes Hopf bifurcation when τ = τ jN , where τ jN = 2 arctan θ + 2jπ ωN , j ∈ Z, and θ is defined in (4.18). EJDE-2021/42 AN AGE-STRUCTURED PREDATOR-PREY MODEL 13 Suppose that λ(τ) = α(τ) + iω(τ) is the root of (4.17) such that α(τ jN ) = 0 and ω(τ jN ) = ωN . Denote G(ω, θ) = [R1(1 + θ2) + 2H(1− θ2)] [2ω(1− θ2) + 2Pθ] − [S1ω(1 + θ2)− 4Hθ] [P (1− θ2)− 4ωθ + S1(1 + θ2)]. (4.19) Therefore from [1, Lemma 2.10 ], we have the following result. Theorem 4.6. If G(ωN , θ) 6= 0, then iωN is a simple root of (4.15) for τ = τ jN , and dReλ(τ) dτ |τ=τjN = dα(τ) dτ |τ=τjN > 0, j ∈ Z, if G(ωN , θ) > 0, dReλ(τ) dτ |τ=τjN = dα(τ) dτ |τ=τjN < 0, j ∈ Z, if G(ωN , θ) < 0. (4.20) 5. Numerical simulations and Summary We illustrate the occurrence of Hopf bifurcation for system (2.1) by numerical simulations. Choose the following set of parameters: µ = 0.8, r = 1.5, K = 1.6, b1 = 0.5, M = 2.5, β(a) := { 2e0.8τ , a ≥ τ, 0, a ∈ (0, τ). (5.1) Then, one can compute P = 0.9667, Q = 0.1333, N = 0, Y = −0.5435, S = −0.8000, R = 1.2000, H = 0. By plotting the graph of F (ω) (see Figure (1) (left)), we have Ω = [0.7134, 1.3606], and thus, the crossing curves of g(ζ) = 0 can be obtained, according to (4.6) and (4.7). Moreover, the quantity R2I1 − R1I2, that determines the crossing directions, can be also calculated by Theorem 4.4, see Figure (1) (right). These curve will determine the stable region of the positive steady state of (1.3), i.e., the bottom-left region in (τ1, τ2)-plane, where the point B locates at. In addition, when the parameters (τ1, τ2) passing through the crossing curves, a periodic solution will be bifurcated from the positive steady state through Hopf bifurcation. It is also observed from Figure (1) (right) that the phenomenon of stability switches takes place, as the parameters (τ1, τ2) moves from the point A (periodic solution) to B (stable steady state) and then to C (periodic solution again). The phase portraits of (1.3), with (τ1, τ2) determined by A, B and C, are shown in Figure 2. When τ1 = τ2, solving (4.17), we obtain the unique solution ω1 ≈ 1.0857, and therefore, owing to Theorem 4.5, we further obtain the first Hopf bifurcation value τ01 ≈ 0.1481. Using (1.4), one has G(ω1, θ1) ≈ 2.6257 > 0. Hence, the solution will approach the positive equilibrium as t → ∞ for τ < τ01 , and increasing τ destabilizes the equilibrium and a periodic solution is bifurcated, see Figure (3). This paper mainly investigates the dynamics of a delayed predator-prey model with the Holling-type IV response and age structure, from the Hopf bifurcation point of view. By converting the model (1.3) into an abstract non-densely defined Cauchy problem, we derive the characteristic equation at the positive steady state, which is a transcendental equation involving two time delays. The crossing curves in (τ1, τ2)-plane, on which the characteristic equation has purely imaginary roots, are obtained. From these curves, we show that the model could exhibit rich dynamics including Hopf bifurcation and stability switches. It also should be mentioned 14 Y. CAI, C. WANG, D. FAN EJDE-2021/42 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ω -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 F (ω ) 0 2 4 6 8 10 12 14 16 18 20 τ1 0 2 4 6 8 10 12 14 16 18 20 τ 2 .A .B .C Figure 1. Left: the graph of F (ω). Right: the crossing curves of g(ζ) = 0. Arrows represent the crossing directions, that is, the region on the end of an arrow has two more characteristic roots with positive real parts. Here, the parameter values are given by (5.1), and the functions V (0) = 1 and u(0, a) = 1.3333e−0.8a are assigned to the initial values. τ1 = 0.3, τ2 = 0.7 τ1 = 1.5, τ2 = 1 τ1 = 2.5, τ2 = 0.8 Figure 2. A solution of (1.3) with different choice of time delays. Figure 3. Positive equilibrium of (1.3) is stable for τ = 0.1 < τ01 (left), and (1.3) has periodic solutions for τ = 0.5 > τ01 (right). The other parameters are the same as Figure 1. that double Hopf bifurcation may occur if the crossing curve intersects with itself or another crossing curve, since the characteristic equation will have two pairs of roots on the imaginary axis. EJDE-2021/42 AN AGE-STRUCTURED PREDATOR-PREY MODEL 15 Acknowledgments. C. Wang was partially supported by the NSFC (11671110) and by Heilongjiang NSF (LH2019A010). D. Fan was partially supported by the Shandong NSF (ZR2020MA010). References [1] S. Chen, J. Shi, J. Wei; Time delay-induced instabilities and Hopf bifurcations in general reaction-diffusion systems, J. Non. Sci., 23 (2013), 1–38. [2] J. M. Cushing, M. Saleem; A predator prey model with age structure, J. Math. Biol., 14 (1982), 231–250. [3] A. Ducrot, P. Magal, S. Ruan; Projectors on the generalized eigenspaces for partial differential equations with time delay, Inf. Dimen. Dyn. Syst., 64 (2013), 353–390. [4] H. Freedman; Deterministic mathematical models in population ecology, Biometrics, 22 (1980), 219–236. [5] H. I. Freedman, G. S. K. Wolkowicz; Predator-prey systems with group defence: the paradox of enrichment revisited, Bull. Math. Biol., 48 (1986), 493–508. [6] K. Gopalsamy; Delayed responses and stability in two-species systems, J. Austr. Math. Soc., 25 (1984), 473–500. [7] K. Gu, S. Niculescu, J. Chen; On stability crossing curves for general systems with two delays, J. Math. Anal. Appl., 311 (2005), 231–252. [8] C. S. Holling; The functional response of predators to prey density and its role in mimicry and population regulation, Mem. Ent. Soc. Can., 46 (1965), 1–60. [9] S. B. Hsu, T. W. Huang; Global stability for a class of predator-prey systems, SIAM J. Appl. Math., 55 (1995), 763–783. [10] M. Kot; Elemnets of mathematical ecology, Cambridge, (2001). [11] Y. Kuang; Delay differential equations with applications in population dynamics, Academic Press, (1993). [12] P. H. Leslie, J. C. Gower; The properties of a stochastic model for the predator-prey type of interaction between two species, Biometrika, 47 (1960), 34–219. [13] X. Lin, H. Wang; Stability analysis of delay differential equations with two discrete delays, Can. Appl. Math. Q., 20 (2012), 519–533. [14] H. Liu; Hopf bifurcation in a delayed predator-prey model with a Holling-type IV functional response, International Conference on Artificial Intelligence & Computational Intelligence, IEEE, 248 (2009), 482–490. [15] Z. Liu, N. Li; Stability and bifurcation in a predator-prey model with age structure and delays, J. Non. Sci., 25 (2015), 937–957. [16] Z. Liu, P. Magal, S. Ruan; Hopf bifurcation for non-densely defined Cauchy problems, Z. Angew. Math. Phys., 62 (2001), 1–15. [17] Z. Liu, P. Magal, S. Ruan; Normal forms for semiliear equation with non-dense domain with applications to age structure models, J. Diff. Equ., 257 (2014), 921–1011. [18] Z. Liu, R. Yuan; Bifurcations in predator-prey systems with nonmonotonic functional re- sponse, Nonl. Anal. RWA, 6 (2005), 187–205. [19] P. Magal, S. Ruan; Center manifolds for semilinear equations with non-dense domain and apllications on Hopf bifurcation in age structured models, Mem. Amer. Math. Soc., 202 (2009), 1–76. [20] P. Magal, S. Ruan; On semilinear Cauchy problems with non-dense domain, Adv. Differ. Equ., 14 (2009), 1041–1084. [21] R. M. May; Time delay versus stability in population models with two or three tropic levels, Ecology, 54 (1973), 315–325. [22] J. Murray; Mathematical biology, Springer, (2003). [23] J. Ren, X. Li; Bifurcations in a seasonally forced predator-prey model with generalized holling type IV functional response, Inter. J. Bif. Chaos, 26 (2016), 1650203. [24] S. Ruan; Absolute stability, conditional stability and bifurcation in Kolmogorov-type predator- prey systems with discrete delays, Quart. Appl. Math., 3 (2001), 159–173. [25] Z. Y. Sun, J. F. Wang; Dynamics and pattern formation in diffusive predator-prey models with predator-axis, Electron. J. Differential Equations, 2020 (2020), no. 36, 1–14. [26] H. R. Thieme; Integrated semigroups and integrated solutions to abstract Cauchy problems, J. Math. Ana. Appl., 152 (1990), 416–447. 16 Y. CAI, C. WANG, D. FAN EJDE-2021/42 [27] Z. Wang, Z. Liu; Hopf bifurcation of an age-structured compartmental pest-pathogen model, J. Math. Anal. Appl., 385 (2012), 1134–1150. [28] D. Xiao, S. Ruan; Multiple bifurcations in a delayed predator-prey system with nonmonotonic functional response, J. Diff. Equ., 176 (2001), 494–510. [29] S. Xu, C. Wang, D. Fan; Stability and bifurcation in an age-structured model with stocking rate and time delays, Disc. Cont. Dyn. Syst. B, 24 (2019), 2535-2549. Yuting Cai School of Mathematics, Harbin Institute of Technology, Harbin, Heilongjiang, 150001, China Email address: 514935056@qq.com Chuncheng Wang (corresponding author) School of Mathematics, Harbin Institute of Technology, Harbin, Heilongjiang, 150001, China Email address: wangchuncheng@hit.edu.cn Dejun Fan School of Mathematics, Harbin Institute of Technology (Weihai) Weihai, Shandong, 264209, China. Email address: dejun fan@163.com 1. Introduction 2. Transformation to Cauchy problem 3. Equilibria and characteristic equation 3.1. Equilibria 3.2. Characteristic equation 4. Hopf bifurcation 4.1. The case 1=2 4.2. The case 1=2 5. Numerical simulations and Summary Acknowledgments References