
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
[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−βu2−auv1+bu,dvdt=θauv1+bu−d2v, | (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, limv→∞g(η,α,v)=η, g(0,α,v)=αα+v, g(η,0,v)=η, g(η,α,0)=1, g(1,α,v)=1, ∂g∂η>0, ∂g∂α>0, ∂g∂v<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(b−cva2+(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−γu−vp+hu+mv),dvdt=v(1−v1+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
{∂u∂t=d1Δu+u(θ1+Kv+r−γu−vp+hu+mv),x∈Ω,t>0,∂v∂t=d2Δv+v(1−v1+qu),x∈Ω,t>0,∂u∂n=∂v∂n=0,x∈∂Ω,t>0,u(x,0)≥0,v(x,0)≥0,x∈Ω, | (1.4) |
where d1≥0 and d2≥0 are, respectively, the diffusion rate of prey and predator. Laplacian operator Δ=∂2∂x2 is in the one-dimensional diffusion, Δ=∂2∂x2+∂2∂y2 is in the two-dimensional diffusion, Δ=∂2∂x2+∂2∂y2+∂2∂z2 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[1−v(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+:0≤u(t)≤r+θγ,0≤v(t)≤1+q(r+θ)γ} as t→+∞.
Define thresholds
R0=θ+r,θ∗=(1+K)[1−r(p+m)]p+m. |
From model (1.3), we denote
{H1(u,v)=u(θ1+Kv+r−γu−vp+hu+mv),H2(u,v)=v(1−v1+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)]+Kq2−Krq(h+mq),A2=γ(1+K)(p+m)+q(1+2K)−r[(1+K)(h+mq)+Kq(p+m)]−θ(h+mq),A3=(1+K)[1−r(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=A21−3A0A2,Δ2=A21A22−27A20A23−4A31A3−4A0A32+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 1−r(p+m)>0, which can be obtained that A3>0 when θ<θ∗, A3<0 when θ>θ∗, or A3=0 when θ=θ∗; if 1p+m≤r, we know that 1−r(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 1−r(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) A3≥0 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+m≤r.
From Lemma 2.2(1), we know that A3≥0 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:
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:
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).
Condition | Equilibria | ||
R0<0 | E0, E2 | ||
R0>0 | θ≤θ∗1p+m>r | Δ2>0 | E0, E1, E2, E∗2, E∗3 |
Δ2=0 | E0, E1, E2, E∗ | ||
Δ2<0 | E0, E1, E2 | ||
θ≤θ∗1p+m>r | Δ2>0 | E0, E1, E2, E∗1, E∗2, E∗3 | |
Δ2=0, Δ1>0 | E0, E1, E2, E∗, E∗1(orE∗3) | ||
Δ2=0, Δ1=0 | E0, E1, E2, E∗∗ | ||
Δ2<0 | E0, E1, E2, E∗3 | ||
r>01p+m≤r | Δ2>0 | E0, E1, E2, E∗1, E∗2, E∗3 | |
Δ2=0, Δ1>0 | E0, E1, E2, E∗, E∗1(orE∗3) | ||
Δ2=0, Δ1=0 | E0, E1, E2, E∗∗ | ||
Δ2<0 | E0, E1, E2, E∗3 |
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 E∗i(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 E∗i(i=1,2,3);
(2) if Δ2=0, and
(i) Δ1>0, then model (1.3) has two positive equilibria E∗ and E∗1(orE∗3);
(ii) Δ1=0, then model (1.3) has a unique positive equilibrium E∗∗=(−A13A0,1−qA13A0);
(3) if Δ2<0, then model (1.3) has a unique positive equilibrium E∗3.
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+m≤r,
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+r−2γu−vp+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=1−2v1+qu. |
We can obtain the following characteristic equation:
λ2−tr(JE)λ+det(JE)=0, |
where
tr(JE)=a11+a22,det(JE)=a11a22−a12a21. |
(Ⅰ) 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+r−1p+m−1,det(JE2)=1p+m−θ1+K−r. |
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{χ∗,χ∗+1−qb12}, then E∗ is stable;
(2) if χ<min{χ∗,χ∗+1−qb12}, then E∗ is unstable;
(3) if χ=χ∗ and qb12>1, then Hopf bifurcation occurs.
Proof. For E∗(u∗,v∗), we obtain
JE∗=(χ∗−χ+1−b12q−1), | (2.3) |
where
χ∗=θ1+Kv∗+r−2γu∗−1,χ=(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∗
λ2−tr(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{χ∗,χ∗+1−qb12}, while E∗ is unstable when χ<min{χ∗,χ∗+1−qb12}.
When χ=χ∗ and qb12>1, we have tr(JE∗)=0 and det(JE∗)>0. Then, we have roots λ1,2=±i√det(JE∗) at χ=χ∗. Choosing χ as a bifurcation parameter (in fact, θ is the control parameter), we have
[d(Reλ(χ))dχ]|λ=i√det(JE∗),χ=χ∗=12d(tr(JE∗))dχ|χ=χ∗=−12≠0. |
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)=u−u∗−u∗lnuu∗+v−v∗−v∗lnvv∗. |
Then, we obtain
dVdt=−γ(u−u∗)2+hv∗(u−u∗)2+(p+hu∗)(u−u∗)(v−v∗)(p+hu∗+mv∗)(p+hu+mv)−Kθ(u−u∗)(v−v∗)(1+Kv∗)(1+Kv)+qv∗(u−u∗)(v−v∗)−(1+qu∗)(v−v∗)2(1+qu∗)(1+qu)≤−(γ−hpm)|u−u∗|2−11+q|v−v∗|2+(Kθ+1p+q+q2(θ+r)γ)|u−u∗||v−v∗|=−B1(|u−u∗|2−B32B1|u−u∗||v−v∗|)−B2|v−v∗|2=−B1(|u−u∗|2−B32B1|u−u∗||v−v∗|−B234B21|v−v∗|2)−(B2−B234B1)|v−v∗|2=−B1(|u−u∗|−B32B1|v−v∗|)2−(B2−B234B1)|v−v∗|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=u−u∗ and V=v−v∗, 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+v∗p+h(u+u∗)+m(v+v∗)),dvdt=(v+v∗)(1−v+v∗1+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∗+2hmu∗v∗(p+hu∗+mv∗)3),g02=θK2u∗(1+Kv∗)3+mu∗(p+hu∗)(p+hu∗+mv∗)3,g21=h(p2+hpu∗−m2v∗2+2hmu∗v∗)(p+hu∗+mv∗)3,g12=θK2(1+Kv∗)3+m(p2−h2u∗2+mpv∗+2hmu∗v∗)(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ξ1−iξ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+(h11−g20)ξ1+(h02−g11)ξ21−g02ξ31]x2+[h11+(2h02−g11)ξ1−2g02ξ21]xy+(h02−g02ξ1)ξ2y2−g03ξ1ξ22y3+[h21+(2h12−g21)ξ1−2g12ξ21−3g03ξ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+(2h12−g21)ξ1−6g12ξ21−3g03ξ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+(h11−g20)ξ1+(h02−g11)ξ21−g02ξ31],Ψxy(0,0,χ∗)=h11+(2h02−g11)ξ1−2g02ξ21,Ψyy(0,0,χ∗)=2(h02−g02ξ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 mi≥1 and μ1 is the smallest eigenvalue. Denote the real-valued Sobolev space S={(u,v)∈H2(0,lπ)×H2(0,lπ):∂u∂n|∂Ω=∂v∂n|∂Ω=0}, and let SC=S⊕Si={s1+s2i:s1,s2∈S}. Denote N0=N∪{0}.
The linearization equation of model (1.4) at E∗ is
(∂u∂t∂v∂t)=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
λ2−Tr(k)λ+Det(k)=0, | (3.2) |
where
Tr(k)=−(d1+d2)k2+χ∗−χ,Det(k)=d1d2k4+[d1−(χ∗−χ+1)d2]k2−χ∗+χ−1+qb12. |
When condition χ>max{χ∗,χ∗+1−qb12} 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 d1d2≥max{0,χ∗−χ+1}, we can obtain Det(k)>0;
(2) if d1d2<χ∗−χ+1, we need to discuss the sign of Δ=(d1−(χ∗−χ+1)d2)2+4d1d2(χ∗−χ+1−qb12),
(ⅰ) 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)+2qb12−2√−qb12(χ∗−χ+1−qb12)<d1d2<χ∗−χ+1, and Det(k)<0 when −(χ∗−χ+1)+2qb12−2√−qb12(χ∗−χ+1−qb12)>d1d2>0.
Hence, the Theorem 3.1 is shown by the properties of E∗ in model (1.4).
Theorem 3.1. When χ>max{χ∗,χ∗+1−qb12} hold, and
(1) if
d1d2≥max{0,χ∗−χ+1}, |
or
−(χ∗−χ+1)+2qb12−2√−qb12(χ∗−χ+1−qb12)<d1d2<χ∗−χ+1, |
then E∗ is locally asymptotically stable;
(2) if
0<d1d2<−(χ∗−χ+1)+2qb12−2√−qb12(χ∗−χ+1−qb12), |
then E∗ is unstable, and the Turing instability occurs.
Remark 3.1. In Figure 3, the blue line represents the equation
d2=1(χ∗−χ+1)+2qb12−2√−qb12(χ∗−χ+1−qb12)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.
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−γu−vp+hu+mv)(u−u∗)]dx+∫Ω[(1−v1+qu)(v−v∗)]dx−d1u∗∫Ω|∇u|2u2dx−d2v∗∫Ω|∇v|2v2dx≤−∫Ω[γ−hpm−12(Kθ+1p+q+q2(θ+r)γ)](u−u∗)2dx−d1u∗∫Ω|∇u|2u2dx−∫Ω[11+q−12(Kθ+1p+q+q2(θ+r)γ)](v−v∗)2dx−d2v∗∫Ω|∇v|2v2dx, |
which is used in the inequality a2+b2≥2ab.
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 dVdt≤0, which implies that the proof of Theorem 3.2 is completed.
Let X=u−u∗ and Y=v−v∗, and we can transform model (1.4) (for brevity, u and v represent X and Y again, respectively) to:
{∂u∂t=d1Δu+(u+u∗)(θ1+K(v+v∗)+r−γ(u+u∗)−v+v∗p+h(u+u∗)+m(v+v∗)),∂v∂t=d2Δv+(v+v∗)(1−v+v∗1+q(u+u∗)). | (3.3) |
Then, model (3.3) can be rewritten as
{∂u∂t=d1Δu+g(u,v,χ),∂v∂t=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 −uxx→u and its corrsponding eigenfunction is ϕn(x)=cosnlx. Then, (3.2) becomes
Λ2−Tn(χ)Λ+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 n∈N0 such that
Tn(χ0)=0,Dn(χ0)>0,Tj(χ0)≠0,Dj(χ0)≠0forj≠n, | (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{χ∗,χ∗+1−qb12}.
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 j≥1, and Dk(χH0)>0 for any k∈N.
(2) When n≥1, since b11(χH0)=0, b′11(χ)<0 for χ0<χ∗, we have 0<b11(χ)<b11(0):=χ∗+1 for χ0<χ∗.
Denote
ln=n√d1+d2χ∗+1,n∈N. |
Then for ln<l<ln+1, and 1≤j≤n, let χHj be the root of χ∗−χ=(d1+d2)n2l2, then it satisfies χH0>χHj>0. Moreover, by b′11(χ)<0 for χ0<χ∗, we obtain
0<χHn<⋯<χH3<χH2<χH1<χH0<χ∗ |
and
Tj(χHj)=0,Ti(χHj)≠0fori≠j,1≤j≤n. |
Now, it is demonstrated that Dn(χHj)>0 for j≠n. 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
d1d2≥max{0,χ∗−χ+1}, |
or
−(χ∗−χ+1)+2qb12−2√−qb12(χ∗−χ+1−qb12)<d1d2<χ∗−χ+1. |
Theorem 3.3. Assume that −(χ∗−χ+1)+2qb12−2√−qb12(χ∗−χ+1−qb12)<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).
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 |
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.
We draw the bifurcation diagram of model (1.3) with the various parameters \theta , r , K , h , p , q , m (see Figures 7–9) 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.
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.
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.
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.
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 4–7).
(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 10–13).
(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. |
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 |
Condition | Equilibria | ||
R_0 < 0 | E_0 , E_2 | ||
R_0 > 0 | \begin{array}{l} \theta \leq \theta^* \\ \frac{1}{p+m}>r \end{array} | \Delta_2 > 0 | E_0 , E_1 , E_2 , E^*_2 , E^*_3 |
\Delta_2=0 | E_0 , E_1 , E_2 , E_* | ||
\Delta_2 < 0 | E_0 , E_1 , E_2 | ||
\begin{array}{l} \theta \leq \theta^* \\ \frac{1}{p+m}>r \end{array} | \Delta_2 > 0 | E_0 , E_1 , E_2 , E^*_1 , E^*_2 , E^*_3 | |
\Delta_2=0 , \Delta_1 > 0 | E_0 , E_1 , E_2 , E_* , E^*_1(or \; E^*_3) | ||
\Delta_2=0 , \Delta_1=0 | E_0 , E_1 , E_2 , E_*^* | ||
\Delta_2 < 0 | E_0 , E_1 , E_2 , E^*_3 | ||
\begin{array}{l} r>0 \\ \frac{1}{p+m} \leq r \end{array} | \Delta_2 > 0 | E_0 , E_1 , E_2 , E^*_1 , E^*_2 , E^*_3 | |
\Delta_2=0 , \Delta_1 > 0 | E_0 , E_1 , E_2 , E_* , E^*_1(or \; E^*_3) | ||
\Delta_2=0 , \Delta_1=0 | E_0 , E_1 , E_2 , E^*_* | ||
\Delta_2 < 0 | E_0 , E_1 , E_2 , E^*_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 |
Condition | Equilibria | ||
R_0 < 0 | E_0 , E_2 | ||
R_0 > 0 | \begin{array}{l} \theta \leq \theta^* \\ \frac{1}{p+m}>r \end{array} | \Delta_2 > 0 | E_0 , E_1 , E_2 , E^*_2 , E^*_3 |
\Delta_2=0 | E_0 , E_1 , E_2 , E_* | ||
\Delta_2 < 0 | E_0 , E_1 , E_2 | ||
\begin{array}{l} \theta \leq \theta^* \\ \frac{1}{p+m}>r \end{array} | \Delta_2 > 0 | E_0 , E_1 , E_2 , E^*_1 , E^*_2 , E^*_3 | |
\Delta_2=0 , \Delta_1 > 0 | E_0 , E_1 , E_2 , E_* , E^*_1(or \; E^*_3) | ||
\Delta_2=0 , \Delta_1=0 | E_0 , E_1 , E_2 , E_*^* | ||
\Delta_2 < 0 | E_0 , E_1 , E_2 , E^*_3 | ||
\begin{array}{l} r>0 \\ \frac{1}{p+m} \leq r \end{array} | \Delta_2 > 0 | E_0 , E_1 , E_2 , E^*_1 , E^*_2 , E^*_3 | |
\Delta_2=0 , \Delta_1 > 0 | E_0 , E_1 , E_2 , E_* , E^*_1(or \; E^*_3) | ||
\Delta_2=0 , \Delta_1=0 | E_0 , E_1 , E_2 , E^*_* | ||
\Delta_2 < 0 | E_0 , E_1 , E_2 , E^*_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 |