Processing math: 57%
Research article Special Issues

Correlated random walks in heterogeneous landscapes: Derivation, homogenization, and invasion fronts

  • Many models for the movement of particles and individuals are based on the diffusion equation, which, in turn, can be derived from an uncorrelated random walk or a position-jump process. In those models, individuals have a location but no well-defined velocity. An alternative, and sometimes more accurate, model is based on a correlated random walk or a velocity-jump process, where individuals have a well defined location and velocity. The latter approach leads to hyperbolic equations for the density of individuals, rather than parabolic equations that result from the diffusion process. Almost all previous work on correlated random walks considers a homogeneous landscape, whereas diffusion models for uncorrelated walks have been extended to spatially varying environments. In this work, we first derive the equations for a correlated random walk in a one-dimensional spatially varying environment with either smooth variation or piecewise constant variation. Then we show how to derive the so-called parabolic limit from the resulting hyperbolic equations. We develop homogenization theory for the hyperbolic equations, and show that taking the parabolic limit and homogenization are commuting actions. We illustrate our results with two examples from ecology: the persistence and spread of a population in a patchy heterogeneous landscape.

    Citation: Frithjof Lutscher, Thomas Hillen. Correlated random walks in heterogeneous landscapes: Derivation, homogenization, and invasion fronts[J]. AIMS Mathematics, 2021, 6(8): 8920-8948. doi: 10.3934/math.2021518

    Related Papers:

    [1] Ruizhi Yang, Dan Jin . Dynamics in a predator-prey model with memory effect in predator and fear effect in prey. Electronic Research Archive, 2022, 30(4): 1322-1339. doi: 10.3934/era.2022069
    [2] Mengting Sui, Yanfei Du . Bifurcations, stability switches and chaos in a diffusive predator-prey model with fear response delay. Electronic Research Archive, 2023, 31(9): 5124-5150. doi: 10.3934/era.2023262
    [3] Fengrong Zhang, Ruining Chen . Spatiotemporal patterns of a delayed diffusive prey-predator model with prey-taxis. Electronic Research Archive, 2024, 32(7): 4723-4740. doi: 10.3934/era.2024215
    [4] Yujia Xiang, Yuqi Jiao, Xin Wang, Ruizhi Yang . Dynamics of a delayed diffusive predator-prey model with Allee effect and nonlocal competition in prey and hunting cooperation in predator. Electronic Research Archive, 2023, 31(4): 2120-2138. doi: 10.3934/era.2023109
    [5] Zhili Zhang, Aying Wan, Hongyan Lin . Spatiotemporal patterns and multiple bifurcations of a reaction- diffusion model for hair follicle spacing. Electronic Research Archive, 2023, 31(4): 1922-1947. doi: 10.3934/era.2023099
    [6] Xiaowen Zhang, Wufei Huang, Jiaxin Ma, Ruizhi Yang . Hopf bifurcation analysis in a delayed diffusive predator-prey system with nonlocal competition and schooling behavior. Electronic Research Archive, 2022, 30(7): 2510-2523. doi: 10.3934/era.2022128
    [7] Wenbin Zhong, Yuting Ding . Spatiotemporal dynamics of a predator-prey model with a gestation delay and nonlocal competition. Electronic Research Archive, 2025, 33(4): 2601-2617. doi: 10.3934/era.2025116
    [8] Weiyu Li, Hongyan Wang . Dynamics of a three-molecule autocatalytic Schnakenberg model with cross-diffusion: Turing patterns of spatially homogeneous Hopf bifurcating periodic solutions. Electronic Research Archive, 2023, 31(7): 4139-4154. doi: 10.3934/era.2023211
    [9] Mengxin He, Zhong Li . Dynamic behaviors of a Leslie-Gower predator-prey model with Smith growth and constant-yield harvesting. Electronic Research Archive, 2024, 32(11): 6424-6442. doi: 10.3934/era.2024299
    [10] Yichao Shao, Hengguo Yu, Chenglei Jin, Jingzhe Fang, Min Zhao . Dynamics analysis of a predator-prey model with Allee effect and harvesting effort. Electronic Research Archive, 2024, 32(10): 5682-5716. doi: 10.3934/era.2024263
  • Many models for the movement of particles and individuals are based on the diffusion equation, which, in turn, can be derived from an uncorrelated random walk or a position-jump process. In those models, individuals have a location but no well-defined velocity. An alternative, and sometimes more accurate, model is based on a correlated random walk or a velocity-jump process, where individuals have a well defined location and velocity. The latter approach leads to hyperbolic equations for the density of individuals, rather than parabolic equations that result from the diffusion process. Almost all previous work on correlated random walks considers a homogeneous landscape, whereas diffusion models for uncorrelated walks have been extended to spatially varying environments. In this work, we first derive the equations for a correlated random walk in a one-dimensional spatially varying environment with either smooth variation or piecewise constant variation. Then we show how to derive the so-called parabolic limit from the resulting hyperbolic equations. We develop homogenization theory for the hyperbolic equations, and show that taking the parabolic limit and homogenization are commuting actions. We illustrate our results with two examples from ecology: the persistence and spread of a population in a patchy heterogeneous landscape.



    Recently, one of the hot issues in the ecosystem being studied is the impact of fear effect on the dynamics of prey and predators in predator-prey models. The fear effect is an inherent psychological reaction of the organisms to increase alertness and respond to danger[1]. It can trigger various anti-predation responses, such as changing the reproduction capacity and strategies [2], changing foraging behaviors and selecting new habitats [3], and reducing the birth and survival rate of offspring [4]. In 2020, Sarkara and Khajanchi [5] proposed the following predator-prey model that introduces the cost of fear into the prey

    {dudt=r0u(η+α(1η)α+v)d1uβu2auv1+bu,dvdt=θauv1+bud2v, (1.1)

    where u and v, respectively, represent the densities of prey and predator. The meaning of the parameters can refer to [5]. In addition, in [5], g(η,α,v)=η+α(1η)α+v stands for a fear function which describes that the prey is affected by the fear of predator. So, we know some characteristics of g(η,α,v) showing limαg(η,α,v)=1, limvg(η,α,v)=η, g(0,α,v)=αα+v, g(η,0,v)=η, g(η,α,0)=1, g(1,α,v)=1, gη>0, gα>0, gv<0, which imply that the inhibitory effect of fear on the birth rate of prey will increase with the increase of predator population and will decrease with the enhancement of the ability to identify the capture of predator.

    In model (1.1), au1+bu stands for the Holling Ⅱ functional response [6], which expresses the prey consumed by each predator per unit time, which only depends on the density of prey without being disturbed by the predator. It is of great significance to use different functional response functions to describe the relationship between prey and predator, which is caused by the difficulty of capturing the prey by the predator and the different absorption conversion rates of the predator. As we all know, the Beddington-DeAngelis (B-D) functional response function [7,8] is described as u1+au+bv, which depends on both the density of prey and predator populations. Assume that b=0, the B-D functional response becomes the Holling Ⅱ functional response, which implies that it is more comprehensive and accurate in describing the interference and handling of populations. In 1960, Leslie and Gower [9] proposed a novel functional response vcu to describe the conversion rate of prey to predator, called the Leslie-Gower (L-G) functional response. Compared with the Holling Ⅱ functional response, the L-G functional response is affected by prey and interfered with by predators [10]. Therefore, the Holling Ⅱ functional response function in model (1.1) can be replaced by the B-D functional response and L-G functional response, respectively, to study further the impact of different functional response functions on populations. Therefore, we consider the following predator-prey model, in which prey is subject to B-D functional response and predator is subject to modified L-G functional response

    {dudt=r0u(η+α(1η)α+v)d0uβu2(1δ)uva1+(1δ)u+e(1δ)v,dvdt=v(bcva2+(1δ)u), (1.2)

    where the meaning of the parameters are the same as model (1.1), b is the growth rate of predator, and δ stands for the strength of prey refuge [11] with δ[0,1). ua1+u+ev represents the B-D functional response. cva2+u represents a modified L-G functional response, where c is the maximum value of the diminishment rate of predator due to prey, and a2 measures the extent to which environment provides predator.

    To explore how the fear effect acts on the populations, many scholars have constructed and studied a large number of models with different functional response functions. For example, Wang et al. [12] as well as Sarkara and Khajanchi [5], both proposed a predator-prey model with fear and Holling Ⅱ. They all found that prey and predator populations were affected by fear effect. Pal et al. [13] constructed a B-D predator-prey model with fear, which found that fear has a destructive effect on stability. Wang et al. [14] investigated an improved L-G predator-prey model with fear. They found that as the degree of fear increases, it will lead to a decrease in population density and the extinction of prey.

    For simplicity, let x=ca1u,y=ca2bv,τ=bt, and model (1.1) can transform into the following simplified model (for simplicity, u, v, t represent x, y, τ again, respectively):

    {dudt=u(θ1+Kv+rγuvp+hu+mv),dvdt=v(1v1+qu), (1.3)

    where r=r0ηd0b, γ=a1βbc, θ=r0(1η)b, K=a2bcα, p=a1ca2(1δ), q=a1(1δ)a2c, h=a1a2, m=be. Here, r can be positive or negative. Now, we give an analysis of the nature of r as follows:

    (1) When r0η>d0, that is, 1>η>d0r0, we have r>0, which means that when the cost of minimum fear is high, the birth rate of prey affected by fear is higher than the mortality rate of prey.

    (2) When r0η<d0, that is, d0r0>η>0, one has r<0, which means that when the cost of minimum fear is small, the birth rate of prey affected by fear effect is lower than the mortality rate of prey.

    (3) When η=1, i.e., without the effect of the fear effect, we have r0>d0, which means that the birth rate is higher than the mortality rate, which is consistent with the actual situation of biological species in the ecosystem.

    (4) θ+r represents the natural growth rate of prey, which is greater than 0, meaning that prey will survive for a long time.

    Then, some natural questions arise from model (1.3):

    ● What are the conditions for the existence and stability of the equilibria?

    ● What are the conditions for the occurrence of Hopf bifurcation, and, if so, how to determine its stability and direction?

    On the other hand, species not only evolve on the timescale but also move randomly on the spatial scale. So, it is inevitable to consider the issues of species in time and space. In 1952, Turing [15] first described the movement of species on time and space scales via using the reaction-diffusion equations. He found that the steady state equilibrium in the spatial model is stable in the absence of diffusion but transforms unstable in the presence of diffusion, which means that diffusion can induce instability of populations [16,17]. In addition to the instability driven by diffusion, the reaction-diffusion system also can be triggered by other mechanisms, such as steady states solutions [18], Hopf bifurcation [19], pattern formation and etc [20]. Han et al. [21] considered a modified L-G predator-prey model with cross-diffusion and indirect predation effect, which shows that cross-diffusion can drive Turing instability and pattern formation. Tiwari et al. [22] proposed and analyzed a B-D predator-prey interaction model with fear. They investigated some properties of bifurcation, such as Hopf bifurcation and pitchfork bifurcation.

    Motivated by these works, we construct a reaction-diffusion predator-prey model as follows

    {ut=d1Δu+u(θ1+Kv+rγuvp+hu+mv),xΩ,t>0,vt=d2Δv+v(1v1+qu),xΩ,t>0,un=vn=0,xΩ,t>0,u(x,0)0,v(x,0)0,xΩ, (1.4)

    where d10 and d20 are, respectively, the diffusion rate of prey and predator. Laplacian operator Δ=2x2 is in the one-dimensional diffusion, Δ=2x2+2y2 is in the two-dimensional diffusion, Δ=2x2+2y2+2z2 is in the three-dimensional diffusion, ΩRn is a bounded domain with smooth boundary Ω, and n is the outward unit normal vector of the boundary Ω as in [23]. There is no population flux across the boundaries owing to homogeneous Neumann boundary conditions. Then, we ask:

    ● What are the critical conditions that determine Turing instability?

    ● What are the conditions for the occurrence of the spatially homogeneous and spatial inhomogeneous periodic solutions?

    ● If spatial homogeneous and spatial inhomogeneous periodic solutions occur, what are the conditions for determining the stability and direction?

    The rest of the paper is organized as follows. In Section 2, we focus on discussing the existence and stability of the equilibria and give the Hopf bifurcation analysis of the nonspatial model (1.3). In Section 3, we investigate the Hopf bifurcation analysis of the spatial model (1.4). In Section 4, we present a series of numerical simulations to reveal the theoretical analysis.

    From the nonspatial model (1.3), one has

    u(t)=u(0)exp(t0[θ1+Kv(s)+rγu(s)v(s)p+hu(s)+mv(s)]ds),v(t)=v(0)exp(t0[1v(s)1+qu(s)]ds).

    We know u(t),v(t)0 based on the above two expressions. Then, R2+={u>0,v>0} is positively invariant for the nonspatial model (1.3).

    Lemma 2.1. All solutions (u(t),v(t)) of model (1.3) are contained in the region U={(u,v)R2+:0u(t)r+θγ,0v(t)1+q(r+θ)γ} as t+.

    Define thresholds

    R0=θ+r,θ=(1+K)[1r(p+m)]p+m.

    From model (1.3), we denote

    {H1(u,v)=u(θ1+Kv+rγuvp+hu+mv),H2(u,v)=v(1v1+qu). (2.1)

    Clearly, model (1.3) has three boundary equilibrium points: E0=(0,0), E1=(θ+rγ,0), E2=(0,1), where E0 and E2 always exist, and E1 existes when R0>0.

    By solving (2.1), we obtain that v=1+qu and

    A0u3+A1u2+A2u+A3=0, (2.2)

    where

    A0=γKq(h+mq),A1=γ[(1+K)(h+mq)+Kq(p+m)]+Kq2Krq(h+mq),A2=γ(1+K)(p+m)+q(1+2K)r[(1+K)(h+mq)+Kq(p+m)]θ(h+mq),A3=(1+K)[1r(p+m)]θ(p+m).

    It is obvious that Eq (2.2) is a third-order algebraic equation, which has one, two, or three positive roots. Hence, discussing the number of positive equilibria of model (1.3) is equivalent to discussing the number of positive roots of Eq (2.2).

    First, let

    f(u)=A0u3+A1u2+A2u+A3,Δ1=A213A0A2,Δ2=A21A2227A20A234A31A34A0A32+18A0A1A2A3.

    From Eq (2.2), we know A0>0, and the signs of A1, A2 and A3 are uncertain. First, we discuss the sign of A3.

    (1) When R0>0 and r>0 hold, that is, θ+r>0 and 1>η>d0r0, if 1p+m>r, we know that 1r(p+m)>0, which can be obtained that A3>0 when θ<θ, A3<0 when θ>θ, or A3=0 when θ=θ; if 1p+mr, we know that 1r(p+m)0, which we can get that A3<0.

    (2) When R0>0 holds, if r<0, that is, d0r0>η>0, which implies that 1r(p+m)>0, then we get that A3>0 when θ<θ, A3<0 when θ>θ, or A3=0 when θ=θ.

    (3) When R0>0 and r=0 hold, we know that A3>0 when 1+Kp+m>θ>0, that is, θ<θ, A3<0 when θ>θ, or A3=0 when θ=θ.

    Therefore, we get the following lemma:

    Lemma 2.2. Regarding the sign of A3, we have

    (1) A30 if (S1):R0>0,θθ,1p+m>r;

    (2) A3<0 if (S2):R0>0,θ>θ,1p+m>r;

    (3) A3<0 if (S3):R0>0,r>0,1p+mr.

    From Lemma 2.2(1), we know that A30 when (S1) holds, which means that Eq (2.2) has at most two positive roots (see Figure 1). Next, by discussing the sign of Δ2, we obtain the number of positive roots of Eq (2.2), seeing the following lemma:

    Figure 1.  The positive roots of f(u)=0 when A30.

    Lemma 2.3. Suppose that (S1) holds. Then, we have

    (1) if Δ2<0, then Eq (2.2) has no positive root;

    (2) if Δ2=0, then Eq (2.2) has a unique positive root;

    (3) if Δ2>0, then Eq (2.2) has two different positive roots.

    From Lemma 2.2(2), (3), we know that A3<0 when (S2) or (S3) holds, which means that Eq (2.2) has at least one positive root and at most three positive roots (see Figure 2). Then, by discussing the sign of Δ1 and Δ2, one has the following lemma:

    Figure 2.  The positive roots of f(u)=0 when A3<0.

    Lemma 2.4. Suppose that (S2) or (S3) holds. We obtain

    (1) if Δ2<0, then Eq (2.2) has a unique positive root;

    (2) if Δ2=0, and

    (i) Δ1=0, then Eq (2.2) has a unique positive root;

    (ii) Δ1>0, then Eq (2.2) has two positive roots;

    (3) if Δ2>0, then Eq (2.2) has three different positive roots.

    Therefore, summarizing the above lemmas, we obtain the following theorem (for convenience, we summarize Theorem 2.1 in Table 1).

    Table 1.  The existence of the equilibria of model (1.3).
    Condition Equilibria
    R0<0 E0, E2
    R0>0 θθ1p+m>r Δ2>0 E0, E1, E2, E2, E3
    Δ2=0 E0, E1, E2, E
    Δ2<0 E0, E1, E2
    θθ1p+m>r Δ2>0 E0, E1, E2, E1, E2, E3
    Δ2=0, Δ1>0 E0, E1, E2, E, E1(orE3)
    Δ2=0, Δ1=0 E0, E1, E2, E
    Δ2<0 E0, E1, E2, E3
    r>01p+mr Δ2>0 E0, E1, E2, E1, E2, E3
    Δ2=0, Δ1>0 E0, E1, E2, E, E1(orE3)
    Δ2=0, Δ1=0 E0, E1, E2, E
    Δ2<0 E0, E1, E2, E3

     | Show Table
    DownLoad: CSV

    Theorem 2.1. (Ⅰ) Model (1.3) has a unique trivial equilibrium E0=(0,0) and one semi-trivial equilibrium E2=(0,1);

    (Ⅱ) if R0>0, model (1.3) also has a semi-trivial equilibrium E1=(u1,0)=(θ+rγ,0);

    (Ⅲ) when (S1) holds, model (1.3) has at most two positive equilibria. Moreover,

    (1) if Δ2>0, then model (1.3) has two different positive equilibria Ei(i=2,3);

    (2) if Δ2=0, then model (1.3) has a unique positive equilibrium E;

    (3) if Δ2<0, then model (1.3) has no positive equilibrium;

    (Ⅳ) when (S2) or (S3) holds, model (1.3) has at least one positive equilibrium and at most three positive equilibria. Moreover,

    (1) if Δ2>0, then model (1.3) has three different positive equilibria Ei(i=1,2,3);

    (2) if Δ2=0, and

    (i) Δ1>0, then model (1.3) has two positive equilibria E and E1(orE3);

    (ii) Δ1=0, then model (1.3) has a unique positive equilibrium E=(A13A0,1qA13A0);

    (3) if Δ2<0, then model (1.3) has a unique positive equilibrium E3.

    For model (1.3), there is at least one positive equilibrium and at most three. Consequently, there are numerous rich and complex dynamical characteristics, which increases the difficulty in studying the dynamical properties of this model. Therefore, in order to alleviate the complexity of this study, we will only consider the case with a unique positive equilibrium, denoted as E. So in the subsequent research, we consistently assume that model (1.3) has a unique positive equilibrium, namely, E(u,v).

    Theorem 2.2. (Ⅰ) E0 and E1 are always hyperbolic saddle;

    (Ⅱ) (1) if R0>0, θθ, 1p+m>r, then E2 is a stable node;

    (2) if one of the conditions hold

    (i) R0>0,θ>θ,1p+m>r;

    (ii) R0>0,1p+mr,

    then E2 is a hyperbolic saddle;

    (3) if R0>0, θ=θ, 1p+m>r, then E2 is a degenerate equilibrium.

    Proof. The corresponding Jacobian matrix for equilibria E=(u,v) of model (1.3) is easily obtained as follows:

    JE=(a11a12a21a22),

    where

    a11=θ1+Kv+r2γuvp+hu+mv+huv(p+hu+mv)2,a12=(θKu(1+Kv)2+u(p+hu)(p+hu+mv)2),a21=qv2(1+qu)2,a22=12v1+qu.

    We can obtain the following characteristic equation:

    λ2tr(JE)λ+det(JE)=0,

    where

    tr(JE)=a11+a22,det(JE)=a11a22a12a21.

    (Ⅰ) For equilibrium E0 and E1, it is obvious that one of the eigenvalues is λ=1>0, then E0 and E1 are hyperbolic saddle.

    (Ⅱ) For equilibrium E2, we obtain

    tr(JE2)=θ1+K+r1p+m1,det(JE2)=1p+mθ1+Kr.

    It is obvious that the sign of det(JE2) is equivalent to the sign of A3. Then, from Lemma 2.2, we know the sign of det(JE2). Hence, when R0>0, θθ, 1p+m>r, we have tr(JE2)<0 and det(JE2)>0, which imply that E2 is a stable node. When (S2) or (S3) holds, we have det(JE2)<0, which implies that E2 is a hyperbolic saddle. When R0>0, θ=θ, 1p+m>r, we have det(JE2)=0, which implies that E2 is a degenerate equilibrium.

    Theorem 2.3. Assume that E exists, and

    (1) if χ>max{χ,χ+1qb12}, then E is stable;

    (2) if χ<min{χ,χ+1qb12}, then E is unstable;

    (3) if χ=χ and qb12>1, then Hopf bifurcation occurs.

    Proof. For E(u,v), we obtain

    JE=(χχ+1b12q1), (2.3)

    where

    χ=θ1+Kv+r2γu1,χ=(p+mv)v(p+hu+mv)2,b12=θKu(1+Kv)2+u(p+hu)(p+hu+mv)2.

    Then we obtain the characteristic equation of E

    λ2tr(JE)λ+det(JE)=0,

    where

    tr(JE)=χχ,det(JE)=χ+χ1+qb12.

    By using the Routh-Hurwitz criterion, the conditions for the stability of E are established, that is to say, E is stable when χ>max{χ,χ+1qb12}, while E is unstable when χ<min{χ,χ+1qb12}.

    When χ=χ and qb12>1, we have tr(JE)=0 and det(JE)>0. Then, we have roots λ1,2=±idet(JE) at χ=χ. Choosing χ as a bifurcation parameter (in fact, θ is the control parameter), we have

    [d(Reλ(χ))dχ]|λ=idet(JE),χ=χ=12d(tr(JE))dχ|χ=χ=120.

    Hence, model (1.3) undergoes a Hopf bifurcation at χ=χ.

    Remark 2.1. For positive equilibrium E, assume that qb12>1, and

    (1) if χ>χ, then E is stable;

    (2) if χ<χ, then E is unstable;

    (3) if χ=χ, then E emerges Hopf bifurcation.

    That means that χ can determine the stability of E.

    Remark 2.2. (1) When m=0 and p=0, the B-D functional response of model (1.3) changes to linear approximations. In this case, the model has a unique positive equilibrium and exhibits limit cycles under fear interference.

    (2) When m=0, the B-D functional response of model (1.3) becomes the Holling Ⅱ functional response. In this case, the model still has a unique positive equilibrium and exhibits limit cycles under fear interference, but shows complexity compared to case (1). Readers can refer to the research results in [14].

    Theorem 2.4. Suppose that R0>0, γ>hpm, and 41+q(γhpm)>(1p+q+rq2γ) hold. If

    γKγ+q2[41+q(γhpm)(1p+q+rq2γ)]>θ,

    then E is globally asymptotically stable.

    Proof. Denote

    V(u,v)=uuulnuu+vvvlnvv.

    Then, we obtain

    dVdt=γ(uu)2+hv(uu)2+(p+hu)(uu)(vv)(p+hu+mv)(p+hu+mv)Kθ(uu)(vv)(1+Kv)(1+Kv)+qv(uu)(vv)(1+qu)(vv)2(1+qu)(1+qu)(γhpm)|uu|211+q|vv|2+(Kθ+1p+q+q2(θ+r)γ)|uu||vv|=B1(|uu|2B32B1|uu||vv|)B2|vv|2=B1(|uu|2B32B1|uu||vv|B234B21|vv|2)(B2B234B1)|vv|2=B1(|uu|B32B1|vv|)2(B2B234B1)|vv|2,

    where B1=γhpm, B2=11+q, B3=Kθ+1p+q+q2(θ+r)γ.

    Since conditions R0>0, γ>hpm, and 41+q(γhpm)>(1p+q+rq2γ) hold, we obtain B2>B234B1 when

    (γhpm)>1+q4(Kθ+1p+q+q2(θ+r)γ)2,

    which imply that dVdt<0. Then, the proof is finished.

    Theorem 2.5. Assume that R0>1, qb12>1, and χ=χ hold, periodic solutions bifurcated from Hopf bifurcation are stable (unstable) and the direction is subcritical (supercritical) when Υ(χ)<0(>0).

    Proof. Let U=uu and V=vv, and we can transform model (1.3) (for brevity, u and v stand for U and V, respectively) to

    {dudt=(u+u)(θ1+K(v+v)+rγ(u+u)v+vp+h(u+u)+m(v+v)),dvdt=(v+v)(1v+v1+q(u+u)). (2.4)

    Rewrite model (2.4) as

    (utvt)=JE(uv)+(g(u,v,χ)h(u,v,χ)), (2.5)

    where JE is denoted in (2.3) and

    {g(u,v,χ)=g20u2+g11uv+g02v2+g30u3+g21u2v+g12uv2+g03v3+O(|(u,v)|4),h(u,v,χ)=h20u2+h11uv+h02v2+h30u3+h21u2v+h12uv2+h03v3+O(|(u,v)|4),

    with

    g20=γ+hχp+hu+mv,g11=(θK(1+Kv)2+p2+hpu+mpv+2hmuv(p+hu+mv)3),g02=θK2u(1+Kv)3+mu(p+hu)(p+hu+mv)3,g21=h(p2+hpum2v2+2hmuv)(p+hu+mv)3,g12=θK2(1+Kv)3+m(p2h2u2+mpv+2hmuv)(p+hu+mv)4,g30=h2χ(p+hu+mv)2,g03=θK3u(1+Kv)4+m2u(p+hu)(p+hu+mv)4,h20=q21+qu,h11=q1+qu,h02=11+qu,h30=q3(1+qu)2,h21=q2(1+qu)2,h12=2q(1+qu)2,h03=0.

    Assuming that JE has two characteristic roots, it can be written as λ1,2(χ)=φ(χ)±iψ(χ), where

    φ(χ)=χχ2,ψ(χ)=χ+χ1+qb12φ2.

    Obviously, eigenvalues λ1(χ) and λ2(χ) are complex conjugate if χ+χ1+qb12>φ2. Then, the eigenvectors of JE corresponding to the eigenvalues of λ(χ)=iψ(χ) at χ=χ are given by

    ξ=[1ξ1iξ2],

    where ξ1=χχ+22b12,ξ2=ψ(χ)b12.

    Based on the transformation (u,v)T=(10ξ1ξ2)(x,y)T, model (2.5) becomes

    [dxdtdydt]=[φ(χ)ψ(χ)ψ(χ)φ(χ)][xy]+[Φ(x,y)Ψ(x,y)], (2.6)

    where

    Φ(x,y)=[g20+g11ξ1+g02ξ21]x2+[g11+2g02ξ1]ξ2xy+g02ξ22y2+[g30+g21ξ1+g12ξ21+g03ξ31]x3+[g12+3g03ξ1]ξ22xy2,Ψ(x,y)=1ξ2[h20+(h11g20)ξ1+(h02g11)ξ21g02ξ31]x2+[h11+(2h02g11)ξ12g02ξ21]xy+(h02g02ξ1)ξ2y2g03ξ1ξ22y3+[h21+(2h12g21)ξ12g12ξ213g03ξ31]x2y.

    Next, we calculate the 1st Lyapunov coefficient as follows:

    Υ(χ)=116[Φxxx+Φxyy+Ψxxy+Ψyyy](0,0,χ)+116ψ(χ)[Φxy(Φxx+Φyy)Ψxy(Ψxx+Ψyy)ΦxxΨxx+ΦyyΨyy](0,0,χ),

    where

    Φxxx(0,0,χ)=6(g30+g21ξ1+g12ξ21+g03ξ31),Φxyy(0,0,χ)=2(g12+3g03ξ1)ξ22,Ψxxy(0,0,χ)=2[h21+(2h12g21)ξ16g12ξ213g03ξ31],Ψyyy(0,0,χ)=6g03ξ1ξ22,Φxx(0,0,χ)=2(g20+g11ξ1+g02ξ21),Φxy(0,0,χ)=(g11+2g02ξ1)ξ2,Φyy(0,0,χ)=2g02ξ22,Ψxx(0,0,χ)=2ξ2[h20+(h11g20)ξ1+(h02g11)ξ21g02ξ31],Ψxy(0,0,χ)=h11+(2h02g11)ξ12g02ξ21,Ψyy(0,0,χ)=2(h02g02ξ1)ξ2,

    owing to

    [(Reφ(χ))χ]χ=χ=12<0.

    Therefore, one gets the 1st Lyapunov coefficient Λ=Υ(χ)φ(χ), which determines the stability and direction of Hopf bifurcating periodic solution.

    Let 0=μ0<μ1<μ2<<μi< be the eigenvalues for the operator Δ subject to the homogeneous Neumann boundary condition on Ω, where μi has multiplicity mi1 and μ1 is the smallest eigenvalue. Denote the real-valued Sobolev space S={(u,v)H2(0,lπ)×H2(0,lπ):un|Ω=vn|Ω=0}, and let SC=SSi={s1+s2i:s1,s2S}. Denote N0=N{0}.

    The linearization equation of model (1.4) at E is

    (utvt)=D(ΔuΔv)+JE(uv), (3.1)

    where D=diag(d1,d2) and JE is denoted in (2.3).

    Then, the characteristic equation of (3.1) is given by

    λ2Tr(k)λ+Det(k)=0, (3.2)

    where

    Tr(k)=(d1+d2)k2+χχ,Det(k)=d1d2k4+[d1(χχ+1)d2]k2χ+χ1+qb12.

    When condition χ>max{χ,χ+1qb12} holds, then we have Tr(0)<0 and Det(0)>0, which imply that Tr(k)<0 always holds.

    Next, we analyze the sign of Det(k). (1) if d1d2max{0,χχ+1}, we can obtain Det(k)>0;

    (2) if d1d2<χχ+1, we need to discuss the sign of Δ=(d1(χχ+1)d2)2+4d1d2(χχ+1qb12),

    (ⅰ) if Δ<0, which implies that Det(k)>0;

    (ⅱ) if Δ>0, which implies that Det(k) must have negative roots, then by calculating, we have Det(k)>0 when (χχ+1)+2qb122qb12(χχ+1qb12)<d1d2<χχ+1, and Det(k)<0 when (χχ+1)+2qb122qb12(χχ+1qb12)>d1d2>0.

    Hence, the Theorem 3.1 is shown by the properties of E in model (1.4).

    Theorem 3.1. When χ>max{χ,χ+1qb12} hold, and

    (1) if

    d1d2max{0,χχ+1},

    or

    (χχ+1)+2qb122qb12(χχ+1qb12)<d1d2<χχ+1,

    then E is locally asymptotically stable;

    (2) if

    0<d1d2<(χχ+1)+2qb122qb12(χχ+1qb12),

    then E is unstable, and the Turing instability occurs.

    Remark 3.1. In Figure 3, the blue line represents the equation

    d2=1(χχ+1)+2qb122qb12(χχ+1qb12)d1,

    and we can see that the red region is the region where Turing instability occurs, and the green region is the region where the steady state solution E is stable. On the other hand, it also shows that when the diffusion of predator is fixed, if the diffusion of prey is smaller, it is more difficult to maintain the stability of populations.

    Figure 3.  Diagram for Turing instability on d1d2 in model (1.4) with θ=4.

    Theorem 3.2. When R0>0 and γ>hpm hold, if

    12(Kθ+1p+q+q2(θ+r)γ)<min{11+q,γhpm},

    then E is globally asymptotically stable.

    Proof. Denote a Lyapunov function

    V=Ω(uuξuξdξ+vvζvζdζ)dx.

    Then

    dVdt=Ω[(θ1+Kv+rγuvp+hu+mv)(uu)]dx+Ω[(1v1+qu)(vv)]dxd1uΩ|u|2u2dxd2vΩ|v|2v2dxΩ[γhpm12(Kθ+1p+q+q2(θ+r)γ)](uu)2dxd1uΩ|u|2u2dxΩ[11+q12(Kθ+1p+q+q2(θ+r)γ)](vv)2dxd2vΩ|v|2v2dx,

    which is used in the inequality a2+b22ab.

    Then, under the Neumann boundary conditions and 11+q>12(Kθ+1p+q+q2(θ+r)γ), γ>hpm+12(Kθ+1p+q+q2(θ+r)γ), we obtain that dVdt0, which implies that the proof of Theorem 3.2 is completed.

    Let X=uu and Y=vv, and we can transform model (1.4) (for brevity, u and v represent X and Y again, respectively) to:

    {ut=d1Δu+(u+u)(θ1+K(v+v)+rγ(u+u)v+vp+h(u+u)+m(v+v)),vt=d2Δv+(v+v)(1v+v1+q(u+u)). (3.3)

    Then, model (3.3) can be rewritten as

    {ut=d1Δu+g(u,v,χ),vt=d2Δv+h(u,v,χ).

    The linearized operator of model (3.3) at (0,0,χ) is given by

    P(χ)=(d1Δ+b11(χ)b12(χ)qd2Δ1),

    where the domain of P(χ) is SC and

    b11(χ)=χχ+1,b12(χ)=θKu(1+Kv)2+u(p+hu)χv(p+mv).

    Let k=n2l2, where n2l2 is the eigenvalue of uxxu and its corrsponding eigenfunction is ϕn(x)=cosnlx. Then, (3.2) becomes

    Λ2Tn(χ)Λ+Dn(χ)=0,n=0,1,2,,

    where

    Tn(χ)=(d1+d2)n2l2+χχ,Dn(χ)=d1d2n4l4+[d1(χχ+1)d2]n2l2χ+χ1+qb12(χ), (3.4)

    and we obtain the eigenvalues

    Λ(χ)=Tn(χ)±T2n(χ)4Dn(χ)2,n=0,1,2,.

    Let χ0 be the possible Hopf bifurcation point satisfying conditions (H1): there exists nN0 such that

    Tn(χ0)=0,Dn(χ0)>0,Tj(χ0)0,Dj(χ0)0forjn, (3.5)

    and for the unique pair of complex eigenvalues ρ(χ)±ω(χ)i near the imaginary axis

    ρ(χ0)0, (3.6)

    where

    ρ(χ)=(d1+d2)n22l2+χχ2,ω(χ)=Dn(χ)ρ2(χ).

    From (3.4), when χ>max{χ,χ+1+qb12}, it is obvious that Tn(χ0)<0 and Dn(χ0)>0, which imply that (u,v) is stable. Thus, any potential bifurcation points shall be in χ<min{χ,χ+1qb12}.

    By calculating, we get

    ρ(χ)=12, (3.7)

    so the transversality condition (3.6) is always held.

    (1) When n=0, let χH0=χ, then we find that χH0 is always a Hopf bifurcation point for l>0 due to T0(χH0)=0, Tj(χH0)<0 for any j1, and Dk(χH0)>0 for any kN.

    (2) When n1, since b11(χH0)=0, b11(χ)<0 for χ0<χ, we have 0<b11(χ)<b11(0):=χ+1 for χ0<χ.

    Denote

    ln=nd1+d2χ+1,nN.

    Then for ln<l<ln+1, and 1jn, let χHj be the root of χχ=(d1+d2)n2l2, then it satisfies χH0>χHj>0. Moreover, by b11(χ)<0 for χ0<χ, we obtain

    0<χHn<<χH3<χH2<χH1<χH0<χ

    and

    Tj(χHj)=0,Ti(χHj)0forij,1jn.

    Now, it is demonstrated that Dn(χHj)>0 for jn. For χ(0,χH0], we get

    Dn(χ)=d1d2n4l4+[d1(χχ+1)d2]n2l2χ+χ1+qb12(χ).

    Obviously, Dn(χ)=0 is a quadratic function of variable n2l2, so we have Dn(χHj)>0 when

    d1d2max{0,χχ+1},

    or

    (χχ+1)+2qb122qb12(χχ+1qb12)<d1d2<χχ+1.

    Theorem 3.3. Assume that (χχ+1)+2qb122qb12(χχ+1qb12)<d1d2<χχ+1 or \frac{d_1}{d_{2}}\geq\max\{0, \chi^*-\chi+1\} hold, and for any l\in(l_n, l_{n+1}] , there are n points \chi^H_j(l) , j\in [1, n] satisfying

    0 < \chi^H_n < \cdots < \chi^H_3 < \chi^H_2 < \chi^H_1 < \chi^H_0 < \chi^*

    such that Hopf bifurcation arises at each \chi = \chi_j^H .

    We apply the normal form theory and center manifold theorem [24,25] to prove Theorems 3.4 and 3.5.

    Theorem 3.4. Assume the condition of Theorem 3.3 is satisfied, the bifurcating periodic solutions of spatial homogeneous are stable (unstable), and the direction is subcritical (supercritical) when \mathrm{Re}(c_1(\chi_0^H)) < 0 (>0) .

    Proof. Let

    \begin{equation*} \begin{aligned} \mathfrak{q} = \left(\mathfrak{a}_0, \mathfrak{b}_0\right)^T = \left(\frac{1+\omega_0\mathrm{i}}{q}, 1 \right)^T, \; \mathfrak{q}^* = \left( \mathfrak{a}_0^*, \mathfrak{b}_0^*\right)^T = \left(\frac{q}{2l\pi\omega_0}\mathrm{i}, \frac{\omega_0-\mathrm{i}}{2l\pi\omega_0} \right)^T, \end{aligned} \end{equation*}

    where \omega_0 = \sqrt{-1+qb_{12}(\chi^H_0)} , such that {P}(\chi_0)\mathfrak{q} = \omega_0\mathfrak{q}\mathrm{i} , P^*(\chi^H_0)\mathfrak{q}^* = -\omega \mathfrak{q}^*\mathrm{i} , \langle \mathfrak{q}^*, \mathfrak{q}\rangle = 1 , and \langle \mathfrak{q}^*, \overline{\mathfrak{q}}\rangle = 0 , where \langle \mathfrak{p}, \mathfrak{q}\rangle = \int_0^{l\pi} \overline{\mathfrak{p}}^T\mathfrak{q } \rm{d}x .

    According to [25], we obtain that

    \begin{equation*} \begin{aligned} \mathcal{E}_{\mathfrak{q}\mathfrak{q}} = \left(\mathfrak{c}_0, \mathfrak{d}_0\right)^T, \; \mathcal{E}_{\mathfrak{q}\bar{\mathfrak{q}}} = \left(\mathfrak{e}_0, \mathfrak{f}_0\right)^T, \; \mathcal{J}_{\mathfrak{q}\mathfrak{q}\bar{\mathfrak{q}}} = \left(\mathfrak{g}_0, \mathfrak{h}_0\right)^T, \end{aligned} \end{equation*}

    where

    \begin{equation*} \begin{aligned} \mathfrak{c}_{0} = &\frac{g_{uu}(1-\omega_0^2)}{q^2} +\frac{2g_{uv}}{q}+g_{vv}+\frac{2(g_{uu}+g_{uv})\omega_0\mathrm{i}}{q}, \\ \mathfrak{d}_{0} = &\frac{h_{uu}(1-\omega_0^2)}{q^2} +\frac{2h_{uv}}{q}+h_{vv}+\frac{2(h_{uu}+h_{uv})\omega_0\mathrm{i}}{q}, \\ \mathfrak{e}_{0} = &\frac{g_{uu}(1-\omega_0^2)}{q^2}+\frac{2g_{uv}}{q}+g_{vv}, \quad \mathfrak{f}_{0} = \frac{h_{uu}(1-\omega_0^2)}{q^2}+\frac{2h_{uv}}{q}+h_{vv}, \\ \mathfrak{g}_{0} = &\frac{g_{uuu}(1-\omega_0^2)}{q^3} +\frac{3g_{uuv}(1-\omega_0^2)}{q^2} +\frac{3g_{uvv}}{q}+g_{vvv}+\bigg[\frac{g_{uuu}(1-\omega_0^2)}{q^3} +\frac{2g_{uuv}}{q^2} +\frac{g_{uvv}}{q}\bigg]\omega_0\mathrm{i}, \\ \mathfrak{h}_{0} = &\frac{h_{uuu}(1-\omega_0^2)}{q^3} +\frac{3h_{uuv}(1-\omega_0^2)}{q^2} +\frac{3h_{uvv}}{q}+\bigg[\frac{h_{uuu}(1-\omega_0^2)}{q^3} +\frac{2h_{uuv}}{q^2} +\frac{h_{uvv}}{q}\bigg]\omega_0\mathrm{i}, \end{aligned} \end{equation*}

    with

    \begin{equation*} \begin{aligned} &g_{uu} = -2\gamma+\frac{2h\chi^H_0}{p+hu^*+mv^*}, \;g_{uuu} = -\frac{6h^2\chi^H_0}{(p+hu^*+mv^*)^2}, \quad h_{uu} = -\frac{2q^2}{(1+qu^*)}, \\ &g_{uv} = -\left(\frac{\theta K}{(1+Kv^*)^2}+\frac{(p^2+hpu^*+pm{v^*}^2+2hmu^*v^*)}{(p+hu^*+mv^*)^3}\right), \quad h_{uv} = \frac{q}{(1+qu^*)}, \\ &g_{vv} = \frac{2\theta K^2u^*}{(1+Kv^*)^3}+\frac{2mu^*(p+ hu^*)\chi^H_0}{v^*(p+mv^*)(p+hu^*+mv^*)}, \quad h_{vv} = -\frac{2}{1+qu^*}, \\ &g_{uuv} = \frac{2h(p^2+hpu^*-m^2{v^*}^2+2hmu^*v^*)}{(p+hu^*+mv^*)^4}, \quad h_{uuu} = \frac{6q^3}{(1+qu^*)^2}, \\ &g_{uvv} = \frac{2\theta K^2}{(1+Kv^*)^3}+\frac{2m(p^2-h^2{u^*}^2+mpv^*+2hmu^*v^*)}{(p+hu^*+mv^*)^4}, \quad h_{uuv} = -\frac{2q^2}{(1+qu^*)^2}, \\ &g_{vvv} = -\left(\frac{6\theta K^3u^*}{(1+Kv^*)^4}+\frac{6m^2u^*(p+ hu^*)\chi^H_0}{v^*(p+mv^*)(p+hu^*+mv^*)^2}\right), \quad h_{uvv} = \frac{2q}{(1+qu^*)^2}. \end{aligned} \end{equation*}

    Then, through straightforward calculations, we can get that

    \begin{equation*} \begin{aligned} \langle \mathfrak{q}^*, \mathcal{E}_{\mathfrak{q}\mathfrak{q}}\rangle = & g_{uu}+g_{uv}+\frac{h_{vv}}{2}+\frac{h_{uu}(1-\omega_0^2- 2q)}{2q^2}-\frac{1}{2\omega_0}\bigg[2g_{uv}+q g_{vv}-h_{vv}\\ &-\frac{h_{uu}(1-\omega_0^2)}{ q^2}+\frac{g_{uu}(1-\omega_0^2)-2h_{uu}\omega_0^2-2h_{uv}(1+\omega_0^2)}{q} \bigg]\mathrm{i}, \\ \end{aligned} \end{equation*}
    \begin{aligned} \left\langle\mathfrak{q}^*, \mathcal{E}_{\mathfrak{q} \bar{\mathfrak{q}}}\right\rangle = & \frac{h_{v v}}{2}+\frac{h_{u u}\left(1-\omega_0^2\right)}{2 q^2}+\frac{h_{u v}}{q} \\ & -\frac{1}{2 \omega_0}\left[2 g_{u v}+q g_{v v}-h_{v v}-\frac{h_{u u}\left(1-\omega_0^2\right)}{q^2}+\frac{g_{u u}\left(1-\omega_0^2\right)-2 h_{u v}}{q}\right] \mathrm{i}, \\ \left\langle\mathfrak{q}^*, \mathcal{J}_{\mathfrak{q} \mathfrak{q} \bar{\mathfrak{q}}}\right\rangle = & -\left[\frac{g_{u v v}}{2}+\frac{g_{u u u}\left(1-\omega_0^2\right)-h_{u u v}\left(1-3 \omega_0^2\right)}{2 q^2}+\frac{g_{u u v}-h_{u v v}}{q}\right] \\ & +\frac{1}{2 \omega_0}\left[3 g_{u v v}+q g_{v v v}+\frac{h_{u u u}}{q^3}+\frac{g_{u u u}\left(1-\omega_0^2\right)+h_{u u v}\left(3-\omega_0^2\right)}{q^2}\right. \\ & \left.+\frac{3 g_{u u v}\left(1-\omega_0^2\right)+h_{u v v}\left(3+\omega_0^2\right)}{q}\right] \mathrm{i} . \end{aligned}

    According to [25], we have

    \begin{equation*} \begin{aligned} \mathrm{Re}(c_1(\chi^H_0)) = &\mathrm{Re}\left(\frac{\mathrm{i}}{2\omega_0}\left(f_{20}f_{11}-2|f_{11}|^2-\frac{|f_{02}|^2}{3}\right)+\frac{f_{21}}{2}\right)\\ = &\mathrm{Re}\left(\frac{\mathrm{i}}{2\omega_0}\langle \mathfrak{q}^*, \mathcal{E}_{\mathfrak{q}\mathfrak{q}}\rangle\cdot\langle \mathfrak{q}^*, \mathcal{E}_{\mathfrak{q}\bar{\mathfrak{q}}}\rangle+\frac{1}{2}\langle \mathfrak{q}^*, \mathcal{J}_{\mathfrak{q}\mathfrak{q}\bar{\mathfrak{q}}}\rangle\right)\\ = &\frac{1}{4\omega_0^2} \left[g_{uu}+g_{uv}+h_{vv}+\frac{h_{uu}(1-\omega_0^2)}{q^2} +\frac{2h_{uv}-h_{uu}}{q}\right]\\ &\quad\times\bigg[2g_{uv}+q g_{vv}-h_{vv}-\frac{h_{uu}(1-\omega_0^2)}{q^2} +\frac{g_{uu}(1-\omega_0^2)-2h_{uv}}{q}\bigg]\\ &-\frac{h_{uu}+h_{uv}}{4q} \bigg[h_{vv}+\frac{h_{uu}(1-\omega_0^2)}{q^2}+\frac{2h_{uv}}{q}\bigg]\\ &-\frac{1}{4}\bigg[g_{uvv} +\frac{{g_{uuu}(1-\omega_0^2)}-h_{uuv}(1-3\omega_0^2)}{q^2} +\frac{2(g_{uuv}-h_{uvv})}{q} \bigg]. \end{aligned} \end{equation*}

    From (3.7), we have \rho'(\chi) < 0 , then the sign of \mathrm{Re}(c_1(\chi^H_0)) can determine the direction and stability of Hopf bifurcation periodic solutions.

    Theorem 3.5. Assume the conditions of Theorem 3.3 are satisfied, the bifurcating periodic solutions of spatial inhomogeneous are stable (unstable), and the direction is subcritical (supercritical) when \mathrm{Re}(c_1(\chi_j^H)) < 0 (>0) .

    Proof. Let

    \begin{equation*} \begin{aligned} &\mathfrak{q} = \left(\mathfrak{a}_j, \mathfrak{b}_j\right)^T\cos\frac{j}{l}x = \left( \frac{1}{q}\left(1+d_2\frac{j^2}{l^2}+{\omega_j}\mathrm{i}\right), 1 \right)^T\cos\frac{j}{l}x, \\ &\mathfrak{q}^* = \left(\mathfrak{a}_j^*, \mathfrak{b}_j^*\right)^T\cos\frac{j}{l}x = \left(\frac{q}{\omega_j l\pi}\mathrm{i}, \frac{1}{l\pi} -\left(\frac{{d_2j^2}}{\omega_jl^3\pi}+\frac{1}{\omega_jl\pi}\right)\mathrm{i} \right)^T\cos\frac{j}{l}x, \end{aligned} \end{equation*}

    where \omega_j = \sqrt{-1+qb_{12}(\chi^H_j)-2d_2\frac{j^2}{l^2}-d_2^2\frac{j^4}{l^4}} , which satisfies {P}(\chi_j)\mathfrak{q} = \omega_j\mathfrak{q}\mathrm{i} , P^*(\chi^H_j)\mathfrak{q}^* = -\omega \mathfrak{q}^*\mathrm{i} , \langle \mathfrak{q}^*, \mathfrak{q}\rangle = 1 , and \langle \mathfrak{q}^*, \overline{\mathfrak{q}}\rangle = 0 .

    From [25], we obtain that

    \begin{equation*} \begin{aligned} \mathcal{E}_{\mathfrak{q}\mathfrak{q}} = \left(\mathfrak{c}_j, \mathfrak{d}_j\right)^T\cos^2\frac{j}{l}x, \; \mathcal{E}_{\mathfrak{q}\bar{\mathfrak{q}}} = \left(\mathfrak{e}_j, \mathfrak{f}_j\right)^T\cos^2\frac{j}{l}x, \; \mathcal{J}_{\mathfrak{q}\mathfrak{q}\bar{\mathfrak{q}}} = \left(\mathfrak{g}_j, \mathfrak{h}_j\right)^T\cos^3\frac{j}{l}x, \end{aligned} \end{equation*}

    where

    \begin{equation*} \begin{aligned} \mathfrak{c}_{j} = &\frac{g_{uu}}{q^2} \left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right] +\frac{2g_{uv} }{q}\left(1+d_2\frac{j^2}{l^2}\right)+g_{vv} +\frac{2\omega_j}{q}\left[\frac{g_{uu}}{q}\left(1+d_2\frac{j^2}{l^2}\right) +g_{uv} \right]\mathrm{i}, \\ \mathfrak{d}_{j} = &\frac{h_{uu}}{q^2} \left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right] +\frac{2h_{uv} }{q}\left(1+d_2\frac{j^2}{l^2}\right)+h_{vv} +\frac{2\omega_j}{q}\left[\frac{h_{uu}}{q}\left(1+d_2\frac{j^2}{l^2}\right) +h_{uv}\right]\mathrm{i}, \\ \mathfrak{e}_{j} = &\frac{g_{uu}}{q^2} \left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right] +\frac{2g_{uv} }{q}\left(1+d_2\frac{j^2}{l^2}\right)+g_{vv}, \\ \mathfrak{f}_{j} = &\frac{h_{uu}}{q^2} \left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right] +\frac{2h_{uv} }{q}\left(1+d_2\frac{j^2}{l^2}\right)+h_{vv}, \\ \mathfrak{g}_{j} = &\frac{g_{uuu}}{q^3}\left(1+d_2\frac{j^2}{l^2}\right) \left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right]+\frac{3g_{uuv}}{q^2} \left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right]+g_{vvv}\\ &+\frac{3g_{uvv}}{q}\left(1+d_2\frac{j^2}{l^2}\right)+\frac{\omega_j}{q}\left(\frac{g_{uuu}}{q^2} \left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right] +\frac{2g_{uuv}}{q}\left(1+d_2\frac{j^2}{l^2}\right) +3g_{uvv}\right)\mathrm{i}, \\ \mathfrak{h}_{j} = &\frac{h_{uuu}}{q^3}\left(1+d_2\frac{j^2}{l^2}\right) \left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right] +\frac{3h_{uuv}}{q^2}\left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right]\\ &+\frac{3h_{uvv}}{q}\left(1+d_2\frac{j^2}{l^2}\right)+\frac{\omega_j}{q}\left(\frac{h_{uuu}}{q^2} \left[\left(1+d_2\frac{j^2}{l^2}\right)^2-\omega_j^2\right] +\frac{2h_{uuv}}{q}\left(1+d_2\frac{j^2}{l^2}\right) +3h_{uvv}\right)\mathrm{i}, \end{aligned} \end{equation*}

    with

    \begin{equation*} \begin{aligned} &g_{uu} = -2\gamma+\frac{2h\chi^H_j}{p+hu^*+mv^*}, \;g_{uuu} = -\frac{6h^2\chi^H_j}{(p+hu^*+mv^*)^2}, \quad h_{uu} = -\frac{2q^2}{1+qu^*}, \\ &g_{uv} = -\left(\frac{\theta K}{(1+Kv^*)^2}+\frac{(p^2+hpu^*+pm{v^*}^2+2hmu^*v^*)}{(p+hu^*+mv^*)^3}\right), \quad h_{uv} = \frac{q}{1+qu^*}, \\ &g_{vv} = \frac{2\theta K^2u^*}{(1+Kv^*)^3}+\frac{2mu^*(p+ hu^*)\chi^H_j}{v^*(p+mv^*)(p+hu^*+mv^*)}, \quad h_{vv} = -\frac{2}{1+qu^*}, \\ &g_{uuv} = \frac{2h(p^2+hpu^*-m^2{v^*}^2+2hmu^*v^*)}{(p+hu^*+mv^*)^4}, \quad h_{uuu} = \frac{6q^3}{(1+qu^*)^2}, \\ &g_{uvv} = \frac{2\theta K^2}{(1+Kv^*)^3}+\frac{2m(p^2-h^2{u^*}^2+mpv^*+2hmu^*v^*)}{(p+hu^*+mv^*)^4}, \quad h_{uuv} = -\frac{2q^2}{(1+qu^*)^2}, \\ \end{aligned} \end{equation*}
    \begin{equation*} \begin{aligned} &g_{vvv} = -\left(\frac{6\theta K^3u^*}{(1+Kv^*)^4}+\frac{6m^2u^*(p+ hu^*)\chi^H_j}{v^*(p+mv^*)(p+hu^*+mv^*)^2}\right), \quad h_{uvv} = \frac{2q}{(1+qu^*)^2}. \end{aligned} \end{equation*}

    According to [25], we have

    \begin{equation*} \begin{aligned} &{[(2\omega_j \mathrm{i}I-P_{2j} (\chi^H_j))]}^{-1} = \frac{\varsigma_1-\varsigma_2\mathrm{i}}{\varsigma_1^2+\varsigma_2^2}\left( \begin{array}{cc} 2\omega_j\mathrm{i}+4d_2\frac{j^2}{l^2}+1&-b_{12}(\chi^H_j)\\ q& 2\omega_j\mathrm{i}+(3d_1-d_2)\frac{j^2}{l^2}-1 \end{array} \right), \\ &{[(2\omega_j\mathrm{i}I-P_{0} (\chi^H_j))]}^{-1} = \frac{\varsigma_3-\varsigma_4\mathrm{i}}{\varsigma_3^2+\varsigma_4^2}\left( \begin{array}{cc} 2\omega_j\mathrm{i}+1&-b_{12}(\chi^H_j)\\ q& 2\omega_j\mathrm{i}-(d_1+d_2)\frac{j^2}{l^2}-1 \end{array} \right), \end{aligned} \end{equation*}

    where

    \begin{equation*} \begin{aligned} &\varsigma_1 = 3d_2(4d_1-d_2)\frac{j^4}{l^4}+3(d_1-d_2)\frac{j^2}{l^2}-3\omega_j^2, &&\varsigma_2 = 6\omega_j(d_1+d_2)\frac{j^2}{l^2}, \\ &\varsigma_3 = d_2^2\frac{j^4}{l^4}-(d_1-d_2)\frac{j^2}{l^2}-3\omega_j^2, &&\varsigma_4 = -2\omega_j(d_1+d_2)\frac{j^2}{l^2}. \end{aligned} \end{equation*}

    Then, we have

    \begin{equation*} \begin{aligned} w_{20} = &\frac{1}{2}\left([(2\omega_j\mathrm{i}I-P_{2j}(\chi^H_j))]^{-1}\cos\frac{2j}{l}x+[(2\omega_j \mathrm{i}I-P_{0}(\chi^H_j))]^{-1}\right)\cdot\left( \begin{array}{cc} \mathfrak{c}_j\\ \mathfrak{d}_j \end{array} \right)\\ = &\frac{\varsigma_1-\varsigma_2\mathrm{i}}{2(\varsigma_1^2+\varsigma_2^2)}\left( \begin{array}{cc} \left(2\omega_j\mathrm{i}+4d_2\frac{j^2}{l^2}+1\right)\mathfrak{c}_j-b_{12}(\chi^H_j)\mathfrak{d}_j\\ q\mathfrak{c}_j+\left(2\omega_j\mathrm{i}+(3d_1-d_2)\frac{j^2}{l^2}-1\right)\mathfrak{d}_j \end{array} \right)\cos\frac{2j}{l}x\\ &+\frac{\varsigma_3-\varsigma_4\mathrm{i}}{2(\varsigma_3^2+\varsigma_4^2)}\left( \begin{array}{cc} (2\omega_j\mathrm{i}+1)\mathfrak{c}_j-b_{12}(\chi^H_j)\mathfrak{d}_j\\ q\mathfrak{c}_j+\left(2\omega_j\mathrm{i}-(d_1+d_2)\frac{j^2}{l^2}-1\right)\mathfrak{d}_j \end{array} \right). \end{aligned} \end{equation*}

    Since

    \begin{equation*} \begin{aligned} &{[(-P_{2j}(\chi^H_j))]}^{-1} = \frac{1}{\varsigma_5}\left( \begin{array}{cc} 4d_2\frac{j^2}{l^2}+1&-b_{12}(\chi^H_j)\\ q& (3d_1-d_2)\frac{j^2}{l^2}-1 \end{array} \right), \\ &{[(-P_{0}(\chi^H_j))]}^{-1} = \frac{1}{\varsigma_6}\left( \begin{array}{cc} 1&-b_{12}(\chi^H_j)\\ q& -(d_1+d_2)\frac{j^2}{l^2}-1 \end{array} \right), \end{aligned} \end{equation*}

    where

    \varsigma_5 = 3d_2(4d_1-d_2)\frac{j^4}{l^4}+3(d_1-d_2)\frac{j^2}{l^2}+\omega_j^2, \; \varsigma_6 = d_2^2\frac{j^4}{l^4}-(d_1-d_2)\frac{j^2}{l^2}+\omega_j^2.

    Then, we get

    \begin{equation*} \begin{aligned} w_{11} = &\frac{1}{2}\left([-P_{2j}(\chi^H_j)]^{-1}\cos\frac{2j}{l}x+[-P_{0}(\chi^H_j)]^{-1}\right)\cdot\left( \begin{array}{cc} \mathfrak{e}_j\\ \mathfrak{f}_j \end{array} \right)\\ = &\frac{1}{2\varsigma_5}\left( \begin{array}{cc} \left(4d_2\frac{j^2}{l^2}+1\right) \mathfrak{e}_j-b_{12}(\chi^H_j)\mathfrak{f}_j\\ q\mathfrak{e}_j+\left((3d_1-d_2)\frac{j^2}{l^2}-1\right)\mathfrak{f}_j \end{array} \right)\cos\frac{2j}{l}x\\ &+\frac{1}{2\varsigma_6}\left( \begin{array}{cc} \mathfrak{e}_j-b_{12}(\chi^H_j)\mathfrak{f}_j\\ q\mathfrak{e}_j-\left((d_u+d_v)\frac{j^2}{l^2}+1\right)\mathfrak{f}_j \end{array} \right). \end{aligned} \end{equation*}

    Then,

    \begin{equation*} \begin{aligned} \mathcal{E}_{w_{20}q} = &\left( \begin{array}{cc} g_{uu}\bar{\mathfrak{a}}_j\tilde{\xi}+g_{uv}(\bar{\mathfrak{a}}_j\tilde{\zeta} +\tilde{\xi})+g_{vv}\tilde{\zeta}\\ h_{uu}\bar{\mathfrak{a}}_j\tilde{\xi}+h_{uv}(\bar{\mathfrak{a}}_j\tilde{\zeta} +\tilde{\xi})+h_{vv}\tilde{\zeta} \end{array} \right)\cos\frac{2j}{l}x\cos\frac{j}{l}x\\ &+\left( \begin{array}{cc} g_{uu}\bar{\mathfrak{a}}_j\tilde{\gamma}+g_{uv}(\bar{\mathfrak{a}}_j\tilde{\chi} +\tilde{\gamma})+g_{vv}\tilde{\chi}\\ h_{uu}\bar{\mathfrak{a}}_j\tilde{\gamma}+h_{uv}(\bar{\mathfrak{a}}_j\tilde{\chi} +\tilde{\gamma})+h_{vv}\tilde{\chi} \end{array} \right)\cos\frac{j}{l}x, \\ \mathcal{E}_{w_{11}q} = &\left( \begin{array}{cc} g_{uu}\mathfrak{a}_j\bar{\xi}+g_{uv}(\mathfrak{a}_j\bar{\zeta}+\bar{\xi})+g_{vv}\bar{\zeta}\\ h_{uu}\mathfrak{a}_j\bar{\xi}+h_{uv}(\mathfrak{a}_j\bar{\zeta}+\bar{\xi})+h_{vv}\bar{\zeta} \end{array} \right)\cos\frac{2j}{l}x\cos\frac{j}{l}x\\ &+\left( \begin{array}{cc} g_{uu}\mathfrak{a}_j\bar{\gamma}+g_{uv}(\mathfrak{a}_j\bar{\chi} +\bar{\gamma})+g_{vv}\bar{\chi}\\ h_{uu}\mathfrak{a}_j\bar{\gamma}+h_{uv}(\mathfrak{a}_j\bar{\chi} +\bar{\gamma})+h_{vv}\bar{\chi} \end{array} \right)\cos\frac{j}{l}x, \end{aligned} \end{equation*}

    where

    \begin{equation*} \begin{aligned} \tilde{\xi} = &\frac{\varsigma_1-\varsigma_2\mathrm{i}}{2(\varsigma_1^2+\varsigma_2^2)} \left[\left(2\omega_j\mathrm{i}+4d_2\frac{j^2}{l^2}+1\right)\mathfrak{c}_j-b_{12}(\chi^H_j)\mathfrak{d}_j\right], \\ \tilde{\zeta} = &\frac{\varsigma_1-\varsigma_2\mathrm{i}}{2(\varsigma_1^2+\varsigma_2^2)} \left[q\mathfrak{c}_j+\left(2\omega_j\mathrm{i}+(3d_1-d_2)\frac{j^2}{l^2}-1\right)\mathfrak{d}_j \right], \\ \tilde{\gamma} = &\frac{\varsigma_3-\varsigma_4\mathrm{i}}{2(\varsigma_3^2+\varsigma_4^2)}\left[ (2\omega_j\mathrm{i}+1)\mathfrak{c}_j-b_{12}(\chi^H_j)\mathfrak{d}_j \right], \\ \tilde{\chi} = &\frac{\varsigma_3-\varsigma_4\mathrm{i}}{2(\varsigma_3^2+\varsigma_4^2)} \left[q\mathfrak{c}_j+\left(2\omega_j\mathrm{i}-(d_1+d_2)\frac{j^2}{l^2}-1\right)\mathfrak{d}_j \right], \\ \bar{\xi} = &\frac{1}{2\varsigma_5}\left[ \left(4d_2\frac{j^2}{l^2}+1\right) \mathfrak{e}_j-b_{12}(\chi^H_j)\mathfrak{f}_j \right], \\ \bar{\zeta} = &\frac{1}{2\varsigma_5}\left[ q\mathfrak{e}_j+\left((3d_1-d_2)\frac{j^2}{l^2}-1\right)\mathfrak{f}_j \right], \\ \bar{\gamma} = &\frac{1}{2\varsigma_6}\left[ \mathfrak{e}_j-b_{12}(\chi^H_j)\mathfrak{f}_j \right], \\ \bar{\chi} = &\frac{1}{2\varsigma_6}\left[ q\mathfrak{e}_j-\left((d_u+d_v)\frac{j^2}{l^2}+1\right)\mathfrak{f}_j \right]. \end{aligned} \end{equation*}

    Note that

    \begin{equation*} \begin{aligned} &\int_0^{l\pi}\cos^2\frac{j}{l}x \rm{d}x = \frac{l\pi}{2}, && \int_0^{l\pi}\left(\cos\frac{2j}{l}x\cos^2\frac{j}{l}x\right) \rm{d}x = \frac{l\pi}{4}, \\ &\int_0^{l\pi}\cos^3\frac{j}{l}x \rm{d}x = 0, && \int_0^{l\pi}\cos^4\frac{j}{l}x \rm{d}x = \frac{3l\pi}{8}. \end{aligned} \end{equation*}

    Thus

    \begin{equation*} \begin{aligned} \langle {\mathfrak{q}}^*, \mathcal{E}_{\mathfrak{q}{\mathfrak{q}}}\rangle = &\langle {\mathfrak{q}}^*, \mathcal{E}_{\mathfrak{q}\bar{\mathfrak{q}}}\rangle = 0, \;\langle {\mathfrak{q}}^*, \mathcal{J}_{\mathfrak{q}\mathfrak{q}\bar{\mathfrak{q}}}\rangle = \frac{3l\pi}{8}(\mathfrak{g}_{j}\bar{\mathfrak{a}}_j^*+\mathfrak{h}_j\bar{\mathfrak{b}}_j^*), \\ \langle {\mathfrak{q}}^*, \mathcal{E}_{w_{20}\bar{\mathfrak{q}}}\rangle = & \frac{l\pi}{4}[\bar{\mathfrak{a}}_j^* (g_{uu}\bar{\mathfrak{a}}_j\tilde{\xi}+g_{uv}(\bar{\mathfrak{a}}_j\tilde{\zeta} +\tilde{\xi})+g_{vv}\tilde{\zeta})+ \bar{\mathfrak{b}}_j^*(h_{uu}\bar{\mathfrak{a}}_j\tilde{\xi}+h_{uv}(\bar{\mathfrak{a}}_j\tilde{\zeta} +\tilde{\xi})+h_{vv}\tilde{\zeta})]\\ &+\frac{l\pi}{2}[\bar{\mathfrak{a}}_j^*(g_{uu}\bar{\mathfrak{a}}_j\tilde{\gamma}+g_{uv}(\bar{\mathfrak{a}}_j\tilde{\chi} +\tilde{\gamma})+g_{vv}\tilde{\chi})+ \bar{\mathfrak{b}}_j^*(h_{uu}\bar{\mathfrak{a}}_j\tilde{\gamma}+h_{uv}(\bar{\mathfrak{a}}_j\tilde{\chi} +\tilde{\gamma})+h_{vv}\tilde{\chi})], \\ \langle {\mathfrak{q}}^*, \mathcal{E}_{w_{11}{\mathfrak{q}}}\rangle = & \frac{l\pi}{4}[\bar{\mathfrak{a}}_j^* (g_{uu}\mathfrak{a}_j\bar{\xi}+g_{uv}(\mathfrak{a}_j\bar{\zeta}+\bar{\xi})+g_{vv}\bar{\zeta})+ \bar{\mathfrak{b}}_j^*(h_{uu}\mathfrak{a}_j\bar{\xi}+h_{uv}(\mathfrak{a}_j\bar{\zeta}+\bar{\xi})+h_{vv}\bar{\zeta})]\\ &+\frac{l\pi}{2}[\bar{\mathfrak{a}}_j^*(g_{uu}\mathfrak{a}_j\bar{\gamma}+g_{uv}(\mathfrak{a}_j\bar{\chi} +\bar{\gamma})+g_{vv}\bar{\chi})+ \bar{\mathfrak{b}}_j^*(h_{uu}\mathfrak{a}_j\bar{\gamma}+h_{uv}(\mathfrak{a}_j\bar{\chi} +\bar{\gamma})+h_{vv}\bar{\chi})]. \end{aligned} \end{equation*}

    According to [25], we have

    \begin{equation*} \begin{aligned} \mathrm{Re}(c_1(\chi^H_j)) = &\mathrm{Re}\left(\frac{\mathrm{i}}{2\omega_0}\langle \mathfrak{q}^*, \mathcal{E}_{\mathfrak{q}\mathfrak{q}}\rangle\cdot\langle \mathfrak{q}^*, \mathcal{E}_{\mathfrak{q}\bar{\mathfrak{q}}}\rangle+\langle \mathfrak{q}^*, \mathcal{E}_{w_{11}{\mathfrak{q}}}\rangle+\frac{1}{2}\langle \mathfrak{q}^*, \mathcal{E}_{w_{20}\bar{\mathfrak{q}}}\rangle+\frac{1}{2}\langle \mathfrak{q}^*, \mathcal{J}_{\mathfrak{q}\mathfrak{q}\bar{\mathfrak{q}}}\rangle\right)\\ = &\mathrm{Re}\langle \mathfrak{q}^*, \mathcal{E}_{w_{11}{\mathfrak{q}}}\rangle+\frac{1}{2}\mathrm{Re}\langle \mathfrak{q}^*, \mathcal{E}_{w_{20}\bar{\mathfrak{q}}}\rangle+\frac{1}{2}\mathrm{Re}\langle \mathfrak{q}^*, \mathcal{J}_{\mathfrak{q}\mathfrak{q}\bar{\mathfrak{q}}}\rangle. \end{aligned} \end{equation*}

    From (3.7), we have \rho'(\chi) < 0 , then the direction and stability of the periodic solution of Hopf bifurcation can be determined by the sign of \mathrm{Re}(c_1(\chi^H_j)) .

    We perform a series of numerical simulation on models (1.3) and (1.4) to illustrate our results. Let initial values (u(0), v(0)) = (2, 3) and parameter values in Table 2. Then, we obtain that R_0 = 4.6 > 0 , \theta^* = 2.1833 > \theta = 1.2 , \chi = 0.6467 < \chi^* = 1.4034 , which imply that model (1.3) has a stable node E_2(23, 0) , an unstable positive equilibrium E^*(u^*, v^*) = (4.1418, 5.9701) , and a stable limit cycle (see Figure 4). When \theta = 4 , we have \theta^* = 2.1833 < \theta = 4 and \chi = 0.5132 > \chi^* = 0.4843 , which imply that a stable node E_2(23, 0) becomes an unstable and the unstable equilibrium becomes E^*(u^*, v^*) = (7.1403, 9.5684) which is stable; moreover, the limit cycle is disappears (see Figure 5).

    Table 2.  The parameter values of model (1.3).
    Parameter Value Parameter Value Parameter Value Parameter Value
    r 3.4 \gamma 0.2 K 0.34 p 0.2
    \theta 1.2 h 0.38 m 0.04 q 1.2

     | Show Table
    DownLoad: CSV
    Figure 4.  Phase portrait and time sequence diagram of model (1.3). In (a), red 'o' represents stable equilibrium E_2 = (23, 0) and green 'o' represents unstable equilibrium E^*(u^*, v^*) = (4.1418, 5.9701) for \chi = 0.6467 < \chi^* = 1.4034 ; the red line represents a stable limit cycle.
    Figure 5.  Phase portrait and time sequence diagram of model (1.3). In (a), green 'o' represents unstable equilibrium E_2 = (23, 0) and red 'o' represents stable equilibrium E^*(u^*, v^*) = (7.1403, 9.5684) for \chi = 0.5132 > \chi^* = 0.4843 . \theta = 4 .

    According to Figures 4 and 5, we find that as \theta ( 1.2\rightarrow 4 ) increases, the positive equilibrium E^*(u^*, v^*) of model (1.3) changes from unstable to stable, and the limit cycle changes from existence to nonexistence, in addition, it also led to an increase in population size. This fully demonstrates that \theta can affect the stability and the population size of model (1.3).

    Note that \theta = \frac{r_0(1-\eta)}{b} , that is, \theta and \eta are negatively correlated. Therefore, if the minimum fear cost \eta is used as a key parameter to explore the dynamic behavior and the trend of populations, it can be replaced by \theta . We find from Figures 6 and 7 that as \theta increases, that is, as the minimum fear cost \eta decreases, the model (1.3) exhibits a limit cycle, which first changes from a large and stable range to a small and unstable range and gradually disappears, while the positive equilibrium gradually changes from an unstable to a globally stable. This means that a lower cost of fear can maintain population stability, causing a decrease in population size but not leading to extinction, while a higher cost of fear can disrupt stability and cause periodic changes in population size.

    Figure 6.  Bifurcation diagrams of model (1.3) in u-v (left) and \theta-u-v (right) parametric space. Hopf critical point H: (u, v, \theta) = (6.9341, 9.3209, 3.7210) , the 1st Lyapunov coefficient = -5.1964\times 10^{-3} ; LP: (u, v, \theta) = (1.3198, 2.5838, 0.1392) ; Neutral saddle point H: (u, v, \theta) = (0.2607, 1.3128, 0.5587) ; BP: (u, v, \theta) = (0, 1, 1.0273) .
    Figure 7.  Bifurcation diagram shows that \theta = \frac{r_0(1-\eta)}{b} can stabilize the model (1.3).

    We draw the bifurcation diagram of model (1.3) with the various parameters \theta , r , K , h , p , q , m (see Figures 79) and observe that the various parameters can stabilize the oscillating model via Hopf bifurcation. It is revealed that: (ⅰ) for \theta , r , m , p , the small parameter values will destroy the stability and produce Hopf bifurcation; (ⅱ) for K , the large parameter value will destroy the stability and produce Hopf bifurcation; (ⅲ) for h , q , the bubble phenomenon [26] appears, which means that the appropriate parameter values will destroy the stability and produce Hopf bifurcation, but the small and large parameter values cannot destroy the stability of populations. Therefore, different parameters can have an impact on the dynamic behavior of model (1.3), mainly interfering with the stability of the equilibria and the existence, stability, and direction of Hopf bifurcations.

    Figure 8.  Bifurcation diagram for prey of model (1.3) with various parameters.
    Figure 9.  Bifurcation diagram for predator of model (1.3) with various parameters.

    For model (1.4), let initial values be (u_0(x), v_0(x)) = (4.5297+0.2\sin(x), 7.0543+0.2\sin(x)) , d_2 = 0.2 , \Omega = (0, 40) , and change the different parameter values \theta and d_1 to discuss the influence of fear and diffusion rate of prey on the dynamics of model (1.4).

    Let \theta = 4 , d_1 = 0.1 , and according to Theorem 3.1, we get that \chi = 0.5132 > \chi^* = 0.4843 and -(\chi^*-\chi+1)+2qb_{12}-2\sqrt{-qb_{12}(\chi^*-\chi+1-qb_{12})} = 0.0969 < \frac{d_1}{d_{2}} = 0.5 < \chi^*-\chi+1 = 0.9711, which imply that the steady state solution of model (1.4) is locally asymptotically stable. If d_1 is selected as 0.01 and 0.001, respectively, according to Theorem 3.1, we have

    \; 0 < \frac{d_1}{d_{2}} = 0.05 < -(\chi^*-\chi+1)+2qb_{12}-2\sqrt{-qb_{12}(\chi^*-\chi+1-qb_{12})} = 0.0969, \\ 0 < \frac{d_1}{d_{2}} = 0.005 < -(\chi^*-\chi+1)+2qb_{12}-2\sqrt{-qb_{12}(\chi^*-\chi+1-qb_{12})} = 0.0969,

    which imply that the Turing instability occurs. This further reflects that diffusion can lead to the occurrence of Turing patterns, providing ideas for studying the morphology of populations and effectively enriching the explosion of species diversity.

    Next, select \theta as 6 and 10, respectively, then we can also get that

    \chi = 0.4793 > \chi^* = 0.1598, \;-(\chi^*-\chi+1)+2qb_{12}-2\sqrt{-qb_{12}(\chi^*-\chi+1-qb_{12})} = 0.0403, \\ \;\chi = 0.444 > \chi^* = -0.3495, \;-(\chi^*-\chi+1)+2qb_{12}-2\sqrt{-qb_{12}(\chi^*-\chi+1-qb_{12})} = 0.003.

    By Theorem 3.1, the steady state solution of model (1.4) is locally asymptotically stable when \theta = 6 or \theta = 10 .

    By comparing Figures 10 and 11 in detail, we find that (ⅰ) the larger \theta and d_1 , the steady state solution of model (1.4) is more stable, indicating that the smaller the fear effect, or the larger the diffusion rate of prey (predator), the more conducive to the stability of prey (predator) population; (ⅱ) the smaller \theta and d_1 , the steady state solution of model (1.4) is more unstable, indicating that the greater the fear effect, or the smaller the diffusion rate of prey, the less conducive to the stability of prey(predator) population. Thus, the conclusion of Theorem 3.1 and Remark 3.1 is verified.

    Figure 10.  Spatiotemporal evolution of prey of model (1.4) with different \theta and d_1 . x = 40 , t = 150 .
    Figure 11.  Spatiotemporal evolution of predator of model (1.4) with different \theta and d_1 . x = 40 , t = 150 .

    In Figure 12, we give the time evolution of model (1.4) with different \theta and d_1 at x = 10 , which corresponds to Figures 10 and 11. These indicated that when the solution of model (1.4) tends to the equilibria or periodic solution, its changing form is calm; but when its changing form undergoes a big sudden change, it means that Turing instability has occurred.

    Figure 12.  Time evolution diagram of model (1.4) with different \theta and d_1 at x = 10 . Red line represents prey population and green line represents predator population.

    In Figure 13(a)(c), the dynamics of the solution of model (1.4) is a smooth oscillatory. In Figure 13(d), the approximations have evolved into the spatially homogeneous steady states u^* and v^* . Then, we have a conclusion that the larger \theta can keep the solution of model (1.4) to the stationary distribution, while a small \theta can make the solution tend to the smooth oscillatory, that is, the smaller the fear effect, the more stable the prey (predator) population will be.

    Figure 13.  Numerical simulation on the dynamics of the solution of model (1.4) with different \theta . d_1 = 0.001 . x = 80 , t = 40 . (a): \theta = 4 , (b): \theta = 5 , (c): \theta = 6 , (d): \theta = 10 .

    In this paper, a predator-prey model with B-D functional response for prey and a modified L-G functional response for predator are established to explore how fear and diffusion affect the spatiotemporal dynamics of the model. For the nonspatial model (1.3):

    (1) The model exhibits rich dynamic properties. In terms of the existence of equilibrium points, the model has three boundary equilibria and at most three different positive equilibria (see Theorem 2.1 or Table 2).

    (2) Discussing the local stability of the positive equilibria (see Theorems 2.2 and 2.3) and the global stability of E^* (see Theorem 2.4), which shows that under the assumption of qb_{12} > 1 , when \chi > \chi^* , E^* is stable; when \chi < \chi^* , E^* is unstable; when \chi = \chi^* , the Hopf bifurcation will occur.

    (3) Exploring the impact of fear on Hopf bifurcation. Fear determines the direction and stability of periodic solutions bifurcated from Hopf bifurcation that are investigated (see Theorem 2.5). The lower cost of fear facilitates the occurrence of Hopf bifurcation, leading to the emergence of limit cycles (see Figures 47).

    (4) Our results indicate that fear can enrich and destroy the dynamic properties of the model.

    For the spatial model (1.4):

    (1) The diffusion of populations led to the occurrence of Turing patterns. When the ratio of the diffusion rate of prey to the diffusion rate of predator exceeds a certain critical value, Turing instability occurs (see Theorem 3.1 and Figure 3), which favors the formation of biodiversity.

    (2) Discussing the stability and direction of spatially homogeneous and inhomogeneous periodic solutions (see Theorems 3.4 and 3.5), which are mainly influenced by the decisive effects of fear and diffusion (see Figures 1013).

    (3) Our results reveal that the larger fear and lower diffusion rate of prey both can destroy the stability of populations and further promote the emergence of periodic solutions, but to a certain extent, and it is advantageous to improve the survival of the species and the formation of biodiversity.

    (4) Our results enrich and develop the dynamics of predator-prey models with diffusion and different functional responses, such as the Holling Ⅱ term [5,12,24], L-G term [27], and B-D term [13,22].

    In the future, we will continue to study the effects of fear and diffusion on the three population model, food web model, partial differential population models with cross diffusion, prey-taxis, and spatial memory. In addition, it can also study the influence of external environmental noise, impulsive, and delay on population dynamics.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This work is supported by the National Natural Science Foundation of China (Nos. 12171004, 12301616), the Doctor Research Project Foundation of Liaoning Province (No. 2023-BS-210), the Liaoning Key Laboratory of Development and Utilization for Natural Products Active Molecules (Nos. LZ202302, LZ202303), the Basic scientific research project of Liaoning Provincial Department of Education (No. JYTMS20231706), and the Doctor Research Project Foundation of Anshan Normal University (No. 23b08).

    The authors declare there is no conflicts of interest.



    [1] Y. Alqawasmeh, F. Lutscher, Movement behaviour of fish, harvesting-induced habitat degradation and the optimal size of marine reserves, Theor. Ecol., 12 (2019), 453–466. doi: 10.1007/s12080-019-0411-x
    [2] Y. Alqawasmeh, F. Lutscher, Persistence and spread of stage-structured populations in heterogeneous landscapes, J. Math. Biol., 78 (2019), 1485–1527. doi: 10.1007/s00285-018-1317-8
    [3] A. Bensoussan, J.-L. Lions, G. Papanicolaou, Asymptotic Analysis for Periodic Structures, Providence: AMS Chelsea Publishing, 2010.
    [4] R. S. Cantrell, C. Cosner, Spatial Ecology via Reaction-Diffusion Equations, Hoboken: Wiley, 2003.
    [5] B. Choi, Y.-J. Kim, Diffusion of biological organisms: Fickian and Fokker–Planck type diffusions, SIAM J. Appl. Math., 79 (2019), 1501–1527. doi: 10.1137/18M1163944
    [6] J. Chung, Y. J. Kim, O. Kwong, C. W. Yoon, Biological advection and cross diffusion with parameter regimes, AIMS Mathematics, 4 (2020), 1721–1744.
    [7] Y. Dolak, T. Hillen, , Cattaneo models for chemotaxis, numerical solution and pattern formation, J. Math. Biol., 46 (2003), 153–170. doi: 10.1007/s00285-002-0173-7
    [8] J. P. Duncan, R. N. Rozum, J. A. Powell, K. M. Kettenring, Multi-scale methods predict invasion speeds in variable landscapes, Theor. Ecol., 10 (2017), 287–303. doi: 10.1007/s12080-017-0329-0
    [9] A. Einstein, Zur Theorie der Brownschen Bewegung, Ann. Phy., 19 (1906), 371–381.
    [10] L. Fahrig, Effect of habitat fragmentation on the extinction threshold: a synthesis, Ecol. Appl., 12 (2002), 346–353.
    [11] R. Filliger, M.-O. Hongler, Supersymmetry in random two-velocity processes, Physica A, 332 (2004), 141–150. doi: 10.1016/j.physa.2003.09.048
    [12] R. A. Fisher, The advance of advantageous genes, Ann. Eugenics, 7 (1937), 355–369. doi: 10.1111/j.1469-1809.1937.tb02153.x
    [13] J. F. Fryxell, Predictive modelling of patch use by terrestrial herbivores. In H.H.T. Prins and F. van Langevelde, editors, Dynamics of Foraging Resource Ecology: Spatial and Temporal, chapter 6A, pages 105–123. Springer, 2008.
    [14] M. J. Garlick, J. A. Powell, M. B. Hooten, L. R. McFarlane, Homogenization of large-scale movement models in ecology, Bull. Math. Biol., 73 (2011), 2088–2108. doi: 10.1007/s11538-010-9612-6
    [15] S. Goldstein, On diffusion by discontinuous movements and the telegraph equation, Quart. J. Mech. Appl. Math., 4 (1951), 129–156. doi: 10.1093/qjmam/4.2.129
    [16] K. P. Hadeler, Nonlinear propagation in reaction transport systems. In S. Ruan and G. Wolkowicz, editors, Differential Equations with Applications to Biology. The Fields Institute Lecture Series, AMS, 1998.
    [17] K. P. Hadeler, Reaction transport systems in biological modelling. In V. Capasso and O. Diekmann, editors, Mathematics Inspired by Biology, Lect. Notes Math. 1714, pages 95–150, Heidelberg, 1999. Springer Verlag.
    [18] K. P. Hadeler, Topics in Mathematical Biology, Heidelberg, Springer, 2018.
    [19] K. P. Hadeler, R. Illner, P. van den Driesche, A disease transport model. In G. Lumer and L. Weiss, editors, Evolution equations and their applications in physical and life sciences, pages 369–385, New York, 2000. Marcel Dekker.
    [20] T. Hillen, A Turing model with correlated random walk, J. Math. Biol., 35 (1996), 49–72. doi: 10.1007/s002850050042
    [21] T. Hillen, Invariance principles for hyperbolic random walk systems, J. Math. Ana. Appl., 210 (1997), 360–374. doi: 10.1006/jmaa.1997.5411
    [22] T. Hillen, Hyperbolic models for chemosensitive movement, Math. Mod. Meth. Appl. Sci., 12 (2002), 1007–1034. doi: 10.1142/S0218202502002008
    [23] T. Hillen, On the {L}^2-moment closure of transport equations: The general case, Disc. Cont. Dyn. Syst. B, 5 (2005), 299–318. doi: 10.3934/dcdsb.2005.5.299
    [24] T. Hillen, Existence theory for correlated random walks on bounded domains, Canad. Appl. Math. Quart., 18 (2010), 1–40.
    [25] E. E. Holmes, Are diffusion models too simple? A comparison with telegraph models of invasion, Am. Nat., 142 (1993), 779–795. doi: 10.1086/285572
    [26] D. D. Joseph, L. Preziosi, Heat waves, Rev. Mod. Phys., 61 (1998), 41–73.
    [27] M. Kac, A stochastic model related to the telegrapher's equation, Rocky MT J. Math., 4 (1956), 497–509.
    [28] F. Lutscher, Modeling alignment and movement of animals and cells, J. Math. Biol., 45 (2002), 234–260. doi: 10.1007/s002850200146
    [29] F. Lutscher, A. Stevens, Emerging patterns in a hyperbolic model for locally interacting cell systems, J. Nonlin. Sci., 12 (2002), 619–640.
    [30] G. A. Maciel, F. Lutscher, How individual response to habitat edges affects population persistence and spatial spread, Am. Nat., 182 (2013), 42–52. doi: 10.1086/670661
    [31] G. A. Maciel, F. Lutscher, Allee effects and population spread in patchy landscapes, J. Biol. Dyn., 9 (2015), 109–123. doi: 10.1080/17513758.2015.1027309
    [32] G. A. Maciel, F. Lutscher, Movement behavior determines competitive outcome and spread rates in strongly heterogeneous landscapes, Theor. Ecol., 11 (2018), 351–365. doi: 10.1007/s12080-018-0371-6
    [33] J. D. Murray, Mathematical Biology, New York: Springer, 1989.
    [34] H. G. Othmer, A continuum model for coupled cells, J. Math. Biol., 17 (1983), 351–369.
    [35] H. G. Othmer, S. R. Dunbar, W. Alt, Models of dispersal in biological systems, J. Math. Biol., 26 (1988), 263–298. doi: 10.1007/BF00277392
    [36] O. Ovaskainen, S. J. Cornell, , Biased movement at a boundary and conditional occupancy times for diffusion processes, J. Appl. Prob., 40 (2003), 557–580. doi: 10.1239/jap/1059060888
    [37] J. R. Potts, T. Hillen, M. A. Lewis, , Edge effects and the spatio-temporal scale of animal movement decisions, Theor. Ecol., 9 (2015), 233–247.
    [38] J. Powell, N. E. Zimmermann, Multiscale analysis of active seed dispersal contributed to resolving Reid's paradox, Ecology, 85 (2004), 490–506. doi: 10.1890/02-0535
    [39] H. Schwetlick, Travelling fronts for multidimensional nonlinear transport equations, Ann. I H Poincarè (C) Ana. Nonlin., 17 (2000), 523–550.
    [40] N. Shigesada, K. Kawasaki, E. Teramoto, Traveling periodic waves in heterogeneous environments, Theor. Popul. Biol., 30 (1986), 143–160. doi: 10.1016/0040-5809(86)90029-8
    [41] P. Turchin, Quantitative Analysis of Movement: measuring and modeling population redistribution of plants and animals, Sunderland: Sinauer Assoc., 1998.
    [42] B. Yurk, C. Cobbold, Homogenization techniques for population dynamics in strongly heterogeneous landscapes, J. Biol. Dyn., 12 (2018), 171–193. doi: 10.1080/17513758.2017.1410238
    [43] E. Zauderer, Correlated random walks, hyperbolic systems and Fokker-Planck equations, Math. Comupt. Model., 17 (1993), 43–47.
  • This article has been cited by:

    1. Boumediene Abdellaoui, Giovanni Siclari, Ana Primo, Fujita exponent for non-local parabolic equation involving the Hardy–Leray potential, 2024, 24, 1424-3199, 10.1007/s00028-024-00984-5
    2. Fatima Zohra Bengrine, Ana Primo, Giovanni Siclari, Existence and non-existence results for parabolic systems with an Hardy-Leray potential, 2025, 550, 0022247X, 129533, 10.1016/j.jmaa.2025.129533
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3847) PDF downloads(186) Cited by(2)

Figures and Tables

Figures(5)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog