# Spatiotemporal Evolution Characteristics of Time-Delay Ecological Competition Systems with Food-Limited and Diffusion

Authors: Feilong Wang; Min Xiao; Zhengxin Wang; Jing Zhao; Gong Chen; Jinde Cao
Journal: Complexity
 (2022)

---

## Abstract

In this paper, we put forward a time-delay ecological competition system with food restriction and diffusion terms under Neumann boundary conditions. For the case without delay, the conditions for local asymptotic stability and Turing instability are constructed. For the case with delay, the existence of Hopf bifurcation is demonstrated by analyzing the root distribution of the corresponding characteristic equations. Furthermore, by using the normal form theory and the center manifold reduction of partial functional differential equations, explicit formulas are obtained to determine the direction of bifurcations and the stability of bifurcating periodic solutions. Finally, some simulation examples are provided to substantiate our analysis.


---

## Body

## 1. Introduction

Based on the diversity of biological populations and ecosystems, various ecological competition systems have been established and widely studied (see [1–5]). In ecological competition systems, the interaction between populations is usually reflected by functional response functions. Results have shown that the predation relationship between populations greatly affects the dynamic behavior of predator-prey systems (see [6–8]). In fact, with the development of the economy, humans will harvest biological populations and develop related biological resources to obtain economic benefits, such as fisheries, forestry, and wildlife management (see [9–11]). In recent years, many scholars have introduced the harvesting terms into biological systems to study the system modeling and related dynamic characteristics [12–15]. May et al. [16] established the following model to analyse the interaction between populations under different harvesting strategies:(1)dUdt=r1U1−UK−αUV−H1,dVdt=r2V1−VβU−H2,where U and V represent the prey and predator densities, respectively. r1 and r2 indicate the intrinsic growth rates. K represents the internal growth limit of the population without predators. α stands for the maximum predation coefficient of predator. β is the quality standard for measuring prey as food. H1 and H2 represent the impact of human harvesting for predator and prey populations. On this basis, scholars have studied the relevant dynamic characteristics of different harvesting terms in system (1). In [17], the authors studied the case of H1=r1h1u and H2=r2h2v in system (1), that is, the constant effort harvesting to both prey and predator population. Particularly, they analyzed the maximum sustainable yield of another population under the condition of limiting the harvest of one population. In [18], the authors discussed the case of H1=h1 and H2=r2h2v for system (1), namely constant yield harvesting for prey and constant effort harvesting for the predator population. At the same time, the relevant dynamic characteristics of the system were analyzed and the effects of different management strategies on the stability of the system were compared. Moreover, in [19, 20], the authors considered the constant yield harvesting of the prey population in system (1), that is, H1=h1 and H2=0. They obtained the conditions of Bogdanov-Takens bifurcation and supercritical/subcritical Hopf bifurcation of codimension 1.Based on the research of the appeal literature, it is found that only the impact of linear capture is taken into account. However, from the perspective of biology and economics, the linear harvesting term cannot accurately reflect the real social activities (harvesting population and developing biological resources) of humans and their impact on the predator-prey system. Thus, the nonlinear harvesting terms have been introduced to model the dynamics of predator-prey systems in recent years (see [21–23]). In [24], a realistic harvesting functional form was proposed as follows:(2)Ht=qEUt1+aE,where q represents the catchability coefficient and E stands for the harvesting effort. The conditions for Hopf bifurcation were derived and a nonlinear state feedback controller was designed to control the Hopf bifurcation.Consider the logistic equation(3)dNdt=rN1−NC.Here,r is the intrinsic growth rate and C represents the carrying capacity. Meanwhile, we can see from (3) that the predicted relation of the specific growth rate, dN/Ndt, to mass, N, is a straight line. However, when considering a measure of the portion of available limiting factors not yet utilized by the population, the linear average growth rate cannot accurately describe the growth trend of the population. Research shows the growth and development of organisms depend on the availability and utilization of food in the living environment, which implies that different degrees of food supply rates will affect the stable age composition of the population and then have an influence on the average growth rate of the population. Thus, when growth limitations are based on the proportion of available resources not utilized, Smith [25] established a food-limited growth function(4)dMdt=rMK−MK+r/cM,where c is positive constant and represents the rate of replacement of mass in the population at saturation. It follows from (4) that the predicted relation of the specific growth rate, dM/Mdt, to mass, M, is not a straight line but a concave curve.Considering the growth limit is based on the proportion of unused available resources and settingH1=qEu/1+fE, H2=0, then system (1) becomes(5)dudt=r1uK−uK+αu−Fu,v−qEu1+fE,dvdt=r2v1−vp+bu,where Fu,v represents the functional response function between populations. p represents that when the density of prey population is low, the predator population can switch to other prey for capture. Then, we will analyse the cases on Fu,v=βuv and Fu,v=e1uv/e2+u in the discussion that follows, respectively.In nature, the survival and development of a population often depend on the amount of food and the living space available. Importantly, the greater the population density, the higher the requirements for the living environment. At the same time, the acquisition of food mainly depends on the living environment of the population, which shows that the change of the population’s living space affects the survival and development of the population to a great extent. Therefore, the population will instinctively migrate and diffuse in space to seek a more suitable environment for survival and development. It is necessary to consider the influence of the diffusion effect on population dynamics in predator-prey systems. Mathematically, the nonlinear system with diffusion will show complex dynamic properties [26–29]. In the reaction-diffusion system proposed by Turing [30, 31], the spatial heterogeneity caused by the internal reaction-diffusion characteristics of the system results in the loss of system symmetry and makes the system self-organize to produce some spatial patterns. The process of pattern formation is called Turing instability (Turing bifurcation). The symmetry of the system is broken, leading to the formation of Turing patterns. Therefore, we call this phenomenon “Turing instability caused by diffusive reaction” [32].On the other hand, time delay has become a factor that cannot be ignored in many biological dynamic systems. A large number of studies have revealed that time delay has an important impact on the dynamic characteristics of biological systems and it is common in predator-prey systems, mainly including mature time delay, capture time delay, and pregnancy time delay [33–35]. Local stability of the system means that if the initial state is adjacent to the equilibrium state, the system will not vibrate, and its state trajectory will eventually fall to the equilibrium state. In particular, Hopf bifurcation is a dynamic bifurcation phenomenon, which shows that when the parameters change near the critical value, the stability of the equilibrium point will change and periodic solutions will be generated in its small neighborhood. Meanwhile, it is found that time delay as a bifurcation parameter can induce Hopf bifurcation [29, 36, 37]. Therefore, in this paper, we consider the pregnancy delay of the predator population and analyse the dynamic characteristics of ecological competition systems.Based on the discussion above, we introduce the diffusion effect of the population [26–29] and the pregnancy delay of the predator population [33–35] into system (5) to explore its impacts on the dynamic characteristics of the ecological competition system, which can be described by(6)∂ut,x∂t=D1Δu+r1uK−uK+αu−Fu,v−qEu1+fE,∂vt,x∂t=D2Δv+r2v1−vp+but−τ,with Neumann boundary conditions and initial conditions(7)ut,x=u0t,x≥0,vt,x=v0t,x≥0,t,x∈−τ,0×Ω¯,∂ut,x∂η=∂vt,x∂η=0,t>0,x∈∂Ω,Here,ut,x,vt,x stand for the population density of prey and predator at a spatial location x and time t, respectively. D1,D2>0 represent the diffusion coefficients associated to u and v, respectively. Δ denotes the Laplacian operator in Rn. Suppose Ω=0,π is a bounded domain with a smooth boundary ∂Ω.Lemma 1.
For any solution of system (6) without delay,(8)limsupt⟶∞ut,x≤r1K1+fE−qEKqEα+r11+fE,limsupt⟶∞vt,x≤p+br1K1+fE−qEKqEα+r11+fE.Proof.
Letut,x,vt,x be a nonnegative solution of system (6) without delay. Note that the functional response Fu,v>0. Then,(9)∂u∂t≤D1Δu+r1uK−uK+αu−qEu1+fE,x∈0,π,t>0,uxt,0=uxt,π=0,t>0,∂v∂t≤D2Δv+r2v1−vp+bu,x∈0,π,t>0,vxt,0=vxt,π=0,t>0.
We can estimate the upper limits ofut,x and vt,x due to the standard comparison principle:(10)ut,x≤r1K1+fE−qEKqEα+r11+fE,vt,x≤p+br1K1+fE−qEKqEα+r11+fE.
In other words, for arbitraryε1>0,ε2>0, there exists positive constants t1,t2 such that(11)ut,x≤r1K1+fE−qEKqEα+r11+fE+ε1,for t≥t1, x∈Ω, and(12)vt,x≤p+br1K1+fE−qEKqEα+r11+fE+ε2,for t≥t2, x∈Ω. So, the conclusion follows immediately.
Compared with the models proposed in [1–5], model (6) is a new measure of the population with a nonlinear average growth rate based on food restriction. At the same time, we consider the reaction-diffusion factors for system (6) to study their influence on the dynamic behaviors of systems. In addition, on the basis of references [16–20], we introduce the nonlinear harvesting term of the harvesting-effort coefficient E into system (6) and study the stability and related dynamic characteristics of the system under Holling I and Holling II functional response functions. Importantly, the addition of pregnancy delay can more accurately reflect the evolution of the population and make the system show more complex dynamic characteristics than the model without delay, which is also a widely concerned direction in the research of biological systems.
The main contributions of this paper can be stated as follows. In Section2, the effects of diffusion on the dynamic behavior of the systems without time delay are investigated and some conditions for system stability and Turing instability are determined. It is found that the appropriate diffusion coefficients will lead to Turing instability. In Section 3, we analyse the stability of equilibrium and Hopf bifurcation in the predator-prey system with time delay as the bifurcation parameter. The condition for Hopf bifurcation is constructed and the expression of the bifurcation threshold is given. In Section 4, we calculated the direction of the bifurcations to get more information about the bifurcations. Section 5 uses some numerical examples to verify the correctness of the previous derivation. Section 6 gives the conclusion of this paper.


## 2. Equilibrium Stability and Turing Instability Analysis

Assume that the predator-prey relationship in system (6) satisfies Holling type I functional response, that is Fu,v=βuv, then we make the following nondimensional scaling transformation [22]:(13)u∼=uK,v∼=βvr1,t∼=r1t,D1∼=D1r1,D2∼=D2r2,h=qr1,θ=r2r1,b∼=r2bβK,p∼=pbK,Dropping the tildes, system (6) can be rewritten by(14)∂ut,x∂t=D1Δu+u1−u1+αu−uv−hEu1+fE,∂vt,x∂t=D2Δv+vθ−bvp+ut−τ.Consider the case of no time delay in system (9), namely, τ=0, then system (9) becomes(15)∂ut,x∂t=D1Δu+u1−u1+αu−uv−hEu1+fE,∂vt,x∂t=D2Δv+vθ−bvp+u.Considering the practical significance of the ecosystem, we are interested in the coexisting equilibrium. In order to obtain the positive equilibrium of system (10), let(16)u1−u1+αu−uv−hEu1+fE=0,vθ−bvp+u=0.Then system (10) has a positive equilibrium E∗=u∗,v∗, where v∗=θp+u∗/b and u∗ satisfies the following quadratic equation:(17)A0u2+A1u+A2=0,where(18)A0=αn>0,A1=1+αnp+n+αS>0,A2=S+np−1,n=θb,S=hE1+fE.Lemma 2.
For equation (12), we come to the following results:(1)
IfA2>0, then equation (12) has no positive roots.(2)
IfA2<0, then equation (12) has a unique positive root u∗=−A1+A12−4A0A2/2A0.Hence, whenA2<0, system (10) has a unique positive equilibrium E∗=u∗,v∗, where v∗=θp+u∗/b.Consider system (10) without diffusion(19)dutdt=u1−u1+αu−uv−hEu1+fE,dvtdt=vθ−bvp+u.The Jacobian matrix of (19) at E∗ is a11a12a21a22, where(20)a11=1−2u∗−αu∗21+αu∗2−v∗−hE1+fE,a12=−u∗<0,a21=bv∗2p+u∗2>0,a22=θ−2bv∗p+u∗=−θ<0.Lemma 3.
WhenA2<0, then we can deduce a11<0.Proof.
Notice that(21)a111+αu∗2=1−2u∗−2u∗2−v∗1+αu∗2−S1+αu∗2=1−2u∗−2u∗2−nc+u∗1+2αu∗+α2u∗2−S1+2αu∗+α2u∗2=−αA0u∗3+A1u∗2+nu∗2+1α2A1−nu∗−A2α=−α−A2u∗+A0u∗2α+2A1u∗α−nu∗α−A2α=−α−A2u∗+A0u∗2+A1u∗α+A1u∗α−nu∗α−A2α=−α−A2u∗+A1u∗α−nu∗α−2A2α=−α1+1αu∗−2A2α.
Thus, we can obtaina111+αu∗2<0 when A2<0. As 1+αu∗2>0, then we can deduce a11<0.
Then from the Jacobian matrix of (13), we get the characteristic equation(22)λ2−B1λ+B2=0,where(23)B1=a11+a22,B2=a11a22−a12a21.Theorem 1.
IfA2<0, then E∗ of system (13) is locally asymptotically stable.Proof.
By Lemma3, we know a11<0 when A2<0. Together with a12<0,a21>0 and a22<0, we can get B1<0 and B2>0. Thus, by the Routh-Hurwitz criterion, E∗ of system (13) is locally asymptotically stable.
It follows from [38] that the Laplacian operator −Δ has the eigenvalue k2k∈0,1,2,… under the homogeneous Neumann boundary condition. And the corresponding eigenfunctions are βk1=ιk,0T,βk2=0,ιkT, where ιk=coskx and βk1,βk2k=0∞ construct a basis of the phase space X and X is defined by(24)X=v∈W2,20,π:v⋅0=v⋅π=0,with the inner product ⋅,⋅. Thus, for system (10), the characteristic equation at E∗ is(25)λ2+C1kλ+C2k=0,where(26)C1k=D1+D2k2−a11−a22,C2k=D1D2k4−a11D2+a22D1k2+a11a22−a12a21.
Hence,(27)λ1+λ2=−C1k=−D1+D2k2+a11+a22,λ1λ2=C2k=D1D2k4−a11D2+a22D1k2+a11a22−a12a21.Theorem 2.
IfA2<0, then E∗ of system (10) is asymptotically stable.Proof.
It follows from Lemma3 that a11<0 when A2<0. Notice that a12<0,a21>0 and a22<0. For equation (14), we obtain(28)λ1+λ2=−C1k=−D1+D2k2+a11+a22<0,λ1λ2=C2k=D1D2k4−a11D2+a22D1k2+a11a22−a12a21>0,k∈N0.
Thus, all roots of (25) have negative real parts for k∈N0, which implies that E∗ of system (10) is asymptotically stable.Remark 1.
Assume that there are no time delays andFu,v=βuv for system (6). It can be seen from Theorems 1 and 2 that when A2<0, the introduction of diffusion terms does not change the stability of E∗, which means that Turing instability does not occur.
Next, we consider the predator-prey relationship among populations in system (6) satisfying Holling type II functional response, that is Fu,v=e1uv/e2+u, where e1 represents the maximum per capita reduction rate of the prey population and e2 stands for the average saturation rate [3]. Then we make the following nondimensional scaling transformation [22]: (29)u∼=uK,v∼=e1vr1,t∼=r1t,D1∼=D1r1,D2∼=D2r2,h=qr1,θ=r2r1,b∼=r2be1K,p∼=pbK,e1∼=1K,e2∼=e2K.
Dropping the tildes, then system (6) without time delay turns into(30)∂ut,x∂t=D1Δu+u1−u1+αu−e1uve2+u−hEu1+fE,∂vt,x∂t=D2Δv+vθ−bvp+u.
System (15) has a positive equilibrium E∗=u∗,v∗, where v∗=θp+u∗/b and u∗ satisfies the following quadratic equation:(31)Y0u2+Y1u+Y2=0,where(32)Y0=αm+αS+1>0,Y1=αmp+m+αe2S+S+e2−1,Y2=e2S+mp−e2,m=e1θb.Lemma 4.
For equation (16), we have the following results:(1)
IfY1<0,Y12−4Y0Y2>0 and Y2>0, then equation (16) has positive roots(33)u1,2=−Y1±Y12−4Y0Y22Y0.(2)
IfY1<0 and Y12−4Y0Y2=0, then equation (16) has a unique positive root u0=−Y1/2Y0.(3)
IfY2<0, then equation (16) has a unique positive root u∗=−Y1+Y12−4Y0Y2/2Y0.Theorem 3.
For system (15), we come to the following results:(1)
IfY1<0,Y12−4Y0Y2>0 and Y2>0, then system (15) has two positive equilibrium points E1u1,v1 and E2u2,v2, where v1,2=θp+u1,2/b.(2)
IfY1<0 and Y12−4Y0Y2=0, then system (15) has a unique positive equilibrium point E0=u0,v0, where v0=θp+u0/b.(3)
IfY2<0, then system (15) has a unique positive equilibrium point E∗=u∗,v∗, where v∗=θp+u∗/b.The linearization of system (15) at E∗ can be represented by(34)∂u∂t∂v∂t=DΔut,xvt,x+Jut,xvt,x,where D=diagD1,D2,J=b11b12b21b22, and(35)b11=1−2u∗−αu∗21+αu∗2−e1e2v∗e2+u∗2−hE1+fE,b12=−e1u∗e2+u∗,b21=bv∗2p+u∗2,b22=θ−2bv∗p+u∗.Thus, the characteristic equation is(36)λ2+C3kλ+C4k=0,where(37)C3k=D1+D2k2−b11−b22,C4k=D1D2k4−b11D2+b22D1k2+b11b22−b12b21.Obviously, for (36), we have(38)λ1+λ2=−C3k=−D1+D2k2+b11+b22,λ1λ2=C4k=D1D2k4−b11D2+b22D1k2+b11b22−b12b21.We make the following assumptions:(39)H1b11+b22<0,H2b11b22−b12b21>0,H3b11D2+b22D1<0.Theorem 4.
IfH1−H3 hold and Y2<0, then system (15) is locally asymptotically stable at E∗.Proof.
IfH1−H3 and Y2<0 hold, for equation (17), we have(40)λ1+λ2=−C3k=−D1+D2k2+b11+b22<0,λ1λ2=C4k=D1D2k4−b11D2+b22D1k2+b11b22−b12b21>0.
Thus, all roots of (36) have negative real parts. By the Routh-Hurwitz criterion, system (15) is locally asymptotically stable at E∗.
Then we can describe system (15) without diffusion term by the following equations:(41)dut∂t=u1−u1+αu−e1uve2+u−hEu1+fE,dvt∂t=vθ−bvp+u.
The characteristic equation of system (18) at E∗ is(42)λ2−C3λ+C4=0,where(43)C3=b11+b22,C4=b11b22−b12b21.Lemma 5.
IfH1,H2, and Y2<0 are satisfied, then E∗ of system (18) is locally asymptotically stable.Let(44)Φ1ε=D1D2ε4−b11D2+b22D1ε2+b11b22−b12b21.We make the following assumptions:(45)H4b11D2+b22D1>0,H5b11D2+b22D1−2D1D2b11b22−b12b21>0.Lemma 6.
IfH2,H4, and H5 hold, then Φ1ε=0 has two positive roots ε1,ε2 .Theorem 5.
Assume thatY2<0 and H1,H2,H4, and H5 hold. Then, diffusion-driven instability occurs for system (15) at E∗ if there exists a k∈N0 such that C4k<0 for 0<ε1<k<ε2.Proof.
It follows from Lemmas5 and 6 that system (18) is stable at E∗ and Φ1ε=0 has two positive roots ε1,ε2 when Y2<0 and H1,H2,H4, and H5 hold. Then for the spatial system (15), the corresponding characteristic equation (17) has an eigenvalue with positive real part if there exists a k∈N0 such that C4k<0 for 0<ε1<k<ε2. This implies the spatial system (15) is unstable at E∗=u∗,v∗, that is, the diffusion-driven instability occurs.Remark 2.
Assume that there are no time delays andFu,v=e1uv/e2+u for system (6). It can be seen from Theorem 5 that the introduction of diffusion terms may change the stability of E∗, resulting in Turing instability.


## 3. Hopf Bifurcation Analysis

In this section, we consider the effect of time delay on the dynamics of the system and get the conditions for Hopf bifurcation of system (9).Linearizing system (9) at E∗=u∗,v∗ [34–36, 39], we have(46)∂u∂t∂v∂t=DΔut,xvt,x+J1ut,xvt,x+J2ut−τ,xvt−τ,x,where D=diagD1,D2, J1=a11a120a22, J2=00a210. Thus, the characteristic equation is(47)detλI−Dk−J1−J2e−λτ=0,where I stands for 2×2 identity matrix and Dk=−k2diagD1,D2,k∈N0. (47) is equivalent to(48)λ2+m1kλ+m2k+m3e−λτ=0,where(49)m1k=D1+D2k2−a11−a22,m2k=D1D2k4−a11D2+a22D1k2+a11a22,m3=−a12a21.Remark 3.
It should be pointed out that this paper adopts the linearization method [29, 36, 37, 39,40] to deal with the dynamics analysis of system (9), including the local stability, Turing instability, and Hopf bifurcation. It is common knowledge that Lyapunov’s second method is important to the stability theory of dynamical systems and control theory. However, this method is not suitable for investigating the dynamics of the ecological competitive system with delay and diffusion proposed in this paper. The Lyapunov stability criterion can only give a sufficient condition for the stability of a system. In this paper, not only the condition of the local stability is established, but also the boundary of stability (the onset of Hopf bifurcation) is determined.Lemma 7.
IfA2<0, then all roots of equation (21) have negative real parts when τ=0.Whenτ>0, assume λ=±iωω>0 is a pair of pure imaginary roots of (48), ω satisfies the following equation:(50)−ω2+im1kω+m2k+m3e−iωτ=0.By separating real and imaginary parts, we obtain(51)−ω2+m2k=−m3cosωτ,m1kω=m3sinωτ.Add the squares of both sides of equation (23) to get(52)ω4+Pkω2+Qk=0,where(53)Pk=m1k2−2m2k=a11−D1k22+a22−D2k22>0,Qk=m2k2−m32=m2k+m3m2k−m3.Letz=ω2, (52) can be converted to(54)z2+Pkz+Qk=0.Then we make the following assumptions:(55)H6a11a22+a12a21>0,H7a11a22+a12a21<0.Theorem 6.
IfH6 and A2<0 hold, then system (9) is locally asymptotically stable at E∗ for all τ≥0.Proof.
By Theorem2, we know m2k+m3=C2k>0 when A2<0. If H6 is satisfied, we obtain(56)m2k−m3=D1D2k4−a11D2+a22D1k2+a11a22+a12a21>0.
Thus, (54) has no positive roots, which implies (48) has no purely imaginary roots. By Lemma 7, we deduce that all roots of (48) have negative real parts for all τ≥0, that is, system (9) is locally asymptotically stable at E∗=u∗,v∗ for all τ≥0.
Let(57)Φ2ε=D1D2ε4−a11D2+a22D1ε2+a11a22+a12a21.
It follows that whenH7 and A2<0 hold, Φ2ε=0 has a unique positive root ε0 having the following form:(58)ε0=a11D2+a22D1+a11D2+a22D12−4D1D2a11a22+a12a212D1D2.
Denote(59)k∗=ε0,ε0∉N,ε0−1,ε0∈N.Lemma 8.
IfH7 and A2<0 are satisfied, equation (21) has a pair of purely imaginary roots ±iωk0≤k≤k∗ at τ=τkj, and(60)τkj=τk0+2jπωk,j∈N0,τk0=1ωkarccosωk2−m2km3,ωk=2m2k−m1k2+m1k2−2m2k2−4m2k2−m322.Theorem 7.
IfH7 and A2<0 hold, we have(61)τk∗j≥τk∗−1j≥τk∗−2j⋅⋅⋅τ1j≥τ0j.Proof.
From equation (28), we obtain(62)ωk2=2m2k−m1k2+m1k2−2m2k2−4m2k2−m322=2m1k2−2m2k/m32−m2k2+m1k2−2m2k2/m32−m2k22+4/m32−m2k2.
By Lemma3, we know that a11<0 when A2<0. Together with a12<0,a21>0 and a22<0, we can get that a11D2+a22D1<0, m2k2−2m2k is increasing and m32−m2k2 is decreasing in k∈0,k∗ when H7 hold. Thus,(63)ωk∗≤ωk∗−1≤ωk∗−2≤⋅⋅⋅ω1≤ω0.
Notice thatm3=−a12a21>0 and m2k is increasing in k∈0,k∗ with a11D2+a22D1<0. We can deduce that ωk2−m2k/m3 is decreasing in k∈0,k∗. Hence, τkj=1/ωkarccosωk2−m2k/m3+2jπ/ωk is increasing in k∈0,k∗.
Asτkj+1>τkj, we obtain(64)τ00=minτkj,0≤k≤k∗,j∈N0.
Letλτ=ρτ+iστ be the root of (48) satisfying ρτkj=0 and στkj=ωk.Theorem 8.
Fork∈0,1,⋅⋅⋅k∗, we have(65)Redλdτ|τ=τkj>0.Proof.
By derivation of equation (21) with respect to the delay τ, we have(66)2λ+m1k−m3τe−λτdλdτ=λm3e−λτ,then(67)dλdτ−1=2λ+m1keλτ−m3τm3λ.
It follows from (50) and (51) that(68)Redλdτ−1|τ=τkj=Re2iωk+m1keiωkτkj−m3τkjiωkm3=2ωkcosωkτkj+m1ksinωkτkjωkm3=2ωk2+m1k2−2m2km32.
Noticingm1k2−2m2k=Pk>0, we obtain(69)Redλdτ−1|τ=τkj=2ωk2+m1k2−2m2km32=2ωk2+Pkm32>0.
Therefore, the transversality condition holds.Theorem 9.
Suppose thatH7 and A2<0 hold. The following statements are valid:(1)
System (9) is locally asymptotically stable at E∗ when τ∈0,τ00 and unstable for τ>τ00.(2)
System (9) undergoes a Hopf bifurcations at E∗ when τ=τ00.


## 4. Direction and Stability of Hopf Bifurcation

The previous analysis has shown that system (9) admits a series of periodic solutions bifurcating from the trivial uniform steady state E∗ at some critical values. Then, in this section, we are concerned with the direction of Hopf bifurcations and the stability of bifurcating periodic solution.Denoteτkjj∈N0,0≤k≤k∗ by τ˜ and introduce the new parameter σ=τ−τ˜. Normalizing the delay τ by the timing-scaling t⟶t/τ, then system (9) can be transformed in the phase space ϑ∗=C−1,0,X(70)V˙=τ˜DΔV+Lτ˜Vt+GVt,σ,where V=u⋅,tv⋅,tT,L:ϑ∗⟶R2 and G:ϑ∗×R+⟶R2 are given, respectively, by(71)Lτ˜ζ=τ˜a11ζ10+a12ζ20a21ζ1−1+a22ζ20,(72)Gζ,σ=σDΔζ+Lσζ+gζ,σ,with gζ,σ=τ˜+σG1ζ,σG2ζ,σT, and(73)G1ζ,σ=ζ10+u∗1−ζ10−u∗1+αζ10+u∗−ζ10+u∗ζ20+v∗−hζ10+u∗1+fE−a11ζ10−a12ζ20,G2ζ,σ=θζ20+v∗−bζ20+v∗2p+ζ1−1+u∗−a21ζ1−1−a22ζ20,for ζ=ζ1ζ2T∈ϑ∗.From the previous discussions, we know thatΛk=iωkτ˜,−iωkτ˜ are a pair of simple purely imaginary eigenvalues of the linear equation(74)V˙t=−τ˜DΔVt+Lτ˜Vt.By the Riesz representation theorem, there exists a2×2 matrix function ηδ,τ,k−1≤δ≤0 such that(75)−τ˜Dk2ζ0+Lτ˜ζ=∫−10dηδ,τ˜,kζδ.Actually, we can take(76)ηδ,τ˜,k=τ˜a11−D1k2a120a22−D2k2,δ=0,0,δ∈−1,0,−τ˜00a210,δ=−1.LetQτ˜ denote the infinitesimal generator of the semigroup and Q∗ be the formal adjoint of Qτ˜ under the bilinear pairing(77)ςδ,ξs=ς0ξ0−∫−10∫ε=0sςε−sdηsξεdε=ς0ξ0+τ˜∫−10ςs+100a210ξsds,for ξ∈C,ς∈C∗=C0,1,R2. Then, we can give the definition of Qτ˜ is(78)Qτ˜=τ˜a11−D1k2a12a21e−λa22−D2k2.From the above analysis, we get thatQτ˜ has a pair of simple purely imaginary eigenvalues ±iωkτ˜. Supposing that p1δ=eiωkτ˜δ1MT,q1s=1Ne−iωkτ˜s, then we can calculate the eigenfunction p1δ of Qτ˜ corresponding to iωkτ˜ and the eigenfunction q1s of Q∗ corresponding to −iωkτ˜. From (78), we derive(79)λI−Qτ˜p10=0|λ=iωkτ˜,(80)λI−Q∗q10T=0|λ=−iωkτ˜.Then, we obtain(81)M=iωk+D1k2−a11a12,N=−a12a22+iωk−D2k2.Letς=ς1,ς2,ξ=ξ1,ξ2T, where(82)ς1δ=p1δ+p1δ¯2=Reeiωkτ˜δReMeiωkτ˜δ,ς2δ=p1δ−p1δ¯2i=Imeiωkτ˜δImMeiωkτ˜δ,for δ∈−1,0, and(83)ξ1s=q1s+q1s¯2=Ree−iωkτ∗sReNe−iωkδτs,ξ2s=q1s−q1s¯2i=Ime−iωkτ∗sImNe−iωkτ˜s,for s∈0,1. Let a1k=a11−D1k2,a2k=a22−D2k2, then by (77), we can calculate(84)ξ1,ς1=1+a1ka2ka2k2+ωk2+τ˜B3a21sinωkτ˜2−τ˜B4a21cosωkτ˜2+sinωkτ˜2ωkτ˜,ξ1,ς2=−a2kωka2k2+ωk2+τ˜B3a21cosωkτ˜2−sinωkτ∗2ωkτ∗+τ˜B4a21sinωkτ˜2,ξ1,ς1=−a1kωka2k2+ωk2+τ˜B3a21cosωkτ˜2+sinωkτ˜2ωkτ˜+τ˜B4a21sinωkτ˜2,ξ1,ς2=ωk2a2k2+ωk2−τ˜B3a21sinωkτ˜2+τ˜B4a21cosωkτ˜2−sinωkτ˜2ωkτ˜,where(85)B3=a12ωka2k2+ωk2,B4=a12a2ka2k2+ωk2.Next, we defineξ,ς=ξi,ςji,j=1,2 and construct a new basis(86)Φ=Φ1,Φ2T=ξ,ς−1ξ.Then, we obtainΦ,ς=I2. Denote fk=βk1,βk2 and a⋅fk=a1βk1+a2βk2, where a=a1,a2T∈ϑ∗. Thus, the center space of the (74) is given by PCNϑ∗ having the following form:(87)PCNψ=ςΦ,<ψ,fk>⋅fk,ψ∈ϑ∗,where(88)ψ,fk=ψ,βk1,ψ,βk2T,a,b=∫0πa1b1¯+a2b2¯dxπ,for a=a1a2T,b=b1b2T.LetQτ˜ be the infinitesimal generator induced by the solution of (74), then system (9) can be rewritten by(89)V˙t=Qτ˜Vt+Y0GVt,σ,where(90)Y0δ=0,−1≤δ<0,I,δ=0.From the decompositionϑ∗=PCNϑ∗⊕PSϑ∗ and (38), where PSϑ∗ denotes the complement subspace of PCNϑ∗ in ϑ∗, the solutions of system (39) can be written as(91)Vt=ςγ1γ2⋅fk+hγ1,γ2,σ,where(92)γ1γ2=Φ,<Vt,fk>.The bifurcation direction of this part is based onσ=0, so the solution of system (31) on the center manifold is given by(93)Vt∗=ςγ1γ2⋅fk+hγ1,γ2,0.Letx=γ1−iγ2 and notice p1=ς1+iς2, then (93) becomes(94)Vt∗=12ςx+x¯ix−x¯⋅fk+Wx,x¯=12p1x+p¯1x¯⋅fk+Wx,x¯,where(95)Wx,x¯=hx+x¯2,−x−x¯2i,0.By [40], it is easy to know x satisfies(96)x.=iωkτ˜x+gx,x¯with(97)gx,x¯=Φ10−iΦ20<GVt∗,0,fk>.By Taylor expansion ofWx,x¯ and gx,x¯, we get(98)Wx,x¯=w20x22+w11xx¯+w02x¯22+⋅⋅⋅,gx,x¯=g20x22+g11xx¯+g02x¯22+g21x2x¯2+⋅⋅⋅.From (94) and (98), we obtain(99)ut0=12x,x¯ιk+w2010x22+w1110xx¯+w0210x¯22+⋅⋅⋅,vt0=12Mx+M¯x¯ιk+w2020x22+w1120xx¯+w0220x¯22+⋅⋅⋅,ut−1=12xe−iωkτ˜+x¯e−iωkτ˜ιk+w201−1x22+w111−1xx¯+w021−1x¯22+⋅⋅⋅,vt−1=12Mxe−iωkτ˜+Mx¯e−iωkτ˜ιk+w202−1x22+w112−1xx¯+w022−1x¯22+⋅⋅⋅,(100)G¯1Vt∗,0=G1τ˜=fuu2ut20+fuvut0vt0+fvv2vt20+fuuu6ut30+fuuv2ut20vt0+fuvv2vt20ut0+fvvv6vt30+⋅⋅⋅,G¯2Vt∗,0=G2τ˜=guu2ut2−1+guvut−1vt0+gvv2vt20+guuu6ut3−1+guuv2ut2−1vt0+guvv2vt20ut−1+gvvv6vt30+⋅⋅⋅,where(101)fuu=−2+2α1+αu∗3,fuuu=6α+6α21+αu∗4,fuv=−1,fvv=fvvv=fuvv=fuuv=0,guu=−2bv∗2p+u∗3,guv=2bv∗p+u∗2,gvv=−2bp+u∗,guuu=6bv∗2p+u∗4,guuv=−4bv∗p+u∗3,guvv=2bp+u∗2,gvvv=0.Thus, we can get(102)G¯1Vt∗,0=x22Ψ20+xx¯Ψ11+x¯22Ψ¯20ιk2+x2x¯2Ψ1ιk+Ψ2ιk3…,G¯2Vt∗,0=x22ν20+xx¯ν11+x¯22ν¯20ιk2+x2x¯2ν1ιk+ν2ιk3…,(103)GVt∗,0,fk=τ˜G¯Vt∗,0,fk=τ˜G¯1Vt∗,0βk1+G¯2Vt∗,0βk2=x22τ˜Ψ20v20Γ+xx¯τ˜Ψ11v11Γ+x¯22τ˜Ψ¯20v¯20Γ+x2x¯2τ˜μ1μ2+⋅⋅⋅,with(104)Γ=1π∫0πιk3xdx,μ1=1π∫0πΨ1ιk2x+Ψ2ιk4xdx,μ2=1π∫0πν1ιk2x+ν2ιk4xdx,where(105)Ψ20=14fuu+2Mfuv+M2fvv,Ψ11=14fuu+Mfuv+M¯fuv+MM¯fvv,Ψ1=w2010fuu+M¯fuv2+w1110fuu+Mfuv+w2020fuv+M¯fvv2+w1120fuv+Mfvv,Ψ2=18fuuu+2M+M¯fuuv+M2+2MM¯fuvv+M2M¯fvvv,and,(106)ν20=14guue−2iωkτ˜+2Mguve−iωkτ˜+M2gvv,ν11=14guu+Meiωkτ˜+M¯e−iωkτ˜guv+MM¯gvv,ν1=w201−1guueiωkτ˜+M¯guv2+w111−1guue−iωkτ˜+Mguv+w2020eiωkτ˜guv+M¯gvv2+w1120guve−iωkτ˜+Mgvv,ν2=18guuue−iωkτ˜+182M+M¯e−2iωkτ˜guuv+18M2M¯gvvv+18M2eiωkτ˜+2MM¯e−iωkτ˜guvv.Denote(107)Φ10−iΦ20=χ1χ2,and notice that(108)1π∫0πιk3xdx=1π∫0πcos3kxdx=0,k=1,2,3⋅⋅⋅.Then we deduce(109)gx,x¯=Φ10−iΦ20<GVt∗,0,fk>=x22Ψ20χ1+ν20χ2Γτ˜+xx¯Ψ11χ1+ν11χ2Γτ˜+x¯22Ψ¯20χ1+ν¯20χ2Γτ˜+x2x¯2μ1χ1+μ2χ2τ˜+⋅⋅⋅.Thus, by equations (43), (45), and (46), we have g20=g02=g11=0 for k=1,2,…. If k=0, we get(110)g20=Ψ20χ1+ν20χ2τ˜,g11=Ψ11χ1+ν11χ2τ˜,g02=Ψ¯20χ1+ν¯20χ2τ˜.Also, fork=N0, we can get(111)g21=μ1χ1+μ2χ2τ˜.So far, we have given the expressions ofg20,g11,g02.It follows from (98) that(112)Wx,x¯=w20xx⋅+w11z⋅x¯+xx¯⋅+w02x¯x¯⋅+⋅⋅⋅,and,(113)Qτ˜Wx,x¯=Qτ˜w20x22+Qτ˜w11xx¯+Qτ˜w02x¯22+⋅⋅⋅.By [40], we know Wx,x¯ satisfies(114)W⋅=Qτ˜W+Zx,x¯,where(115)Zx,x¯=Z20x22+Z11xx¯+Z02x¯22+⋅⋅⋅=Y0δGVt∗,0−ςΦ,<Y0δGVt∗,0,fk>⋅fk.Substituting (96) into the derivative of (98) and comparing the coefficients with (114) and (115), we have(116)2iωkτ˜−Qτ˜w20δ=Z20δ,(117)Qτ˜w11δ=−Z11δ.From (115), we have that for θ∈−1,0,(118)Zx,x¯,δ=−ςδΦ0<GVt∗,0,fk>⋅fk=−p1δ+p1δ¯2,p1δ−p1δ¯2iΦ10,Φ20T⋅<GVt∗,0,fk>⋅fk=−12g20p1δ+g¯02p1δ¯x22⋅fk−12g11p1δ+g¯11p1δ¯xx¯⋅fk+⋯.Hence, forδ∈−1,0,(119)Z20δ=0,k∈N,−12g20p1δ+g¯02p1δ¯⋅fk,k=0,Z20δ=0,k∈N,−12g20p1δ+g¯11p1δ¯⋅fk,k=0,(120)Zx,x¯0=GVt∗,0−ςΦ,<GVt∗,0,fk>⋅fk,where(121)Z200=τ˜Ψ20ν20cos2kxk∈N,τ˜Ψ20ν20−12g20p10+g¯02p10¯⋅fkk=0,Z110=τ˜Ψ11ν11cos2kxk∈N,τ˜Ψ11ν11−12g11p10+g¯11p10¯⋅fkk=0.Thus, we can get from (78) and (116) that(122)w20⋅δ=2iωkτ˜w20δ+12g20p1δ+g¯02p1δ¯⋅fk,w11⋅δ=−i2ωkτ˜g11p1δ+g¯11p1δ¯fk,for −1≤δ<0. That is,(123)w20δ=i2ωkτ˜g20p1δ+g¯023p1δ¯⋅fk+E1e2iωkτ˜δ,w11δ=i2ωkτ˜g¯11p1δ¯−g11p1δ⋅fk+E2,where −1≤δ<0 and E1,E2∈R2 are constant vectors. From the above discussion and (116), we have(124)2iωkE1−Qτ˜−Lτ˜E1e2iωkτ˜=τ˜Ψ20ν20cos2kx,−Qτ˜E2−Lτ˜E2=τ˜Ψ11ν11cos2kx.Then, we obtain(125)E1=τ˜2iωk+D1k2−a11−a12−a21e−2iωkτ˜2iωk+D2k2−a22−1Ψ20ν20ιk2,E2=τ˜D1k2−a11−a12−a21D2k2−a22−1Ψ11ν11ιk2.Therefore, eachg20,g11,g02,g21 can be determined and we can evaluate the following values determining the direction and stability of bifurcating periodic orbits:(126)C10=i2ωkτ˜g11g20−2g112−13g022+12g21,κ2=−ReC10Reλ′τ˜,β2=2ReC10,T2=−ImC10+γ2Imλ′τ˜ωkτ˜.Theorem 10.
For any critical valueτkj, we have(1)
Whenκ2>0κ2<0, the Hopf bifurcation is supercritical (subcritical), that is, the bifurcating periodic solutions exist for τ>τkjτ<τkj.(2)
Whenβ2>0β2<0, the bifurcating periodic solutions on the center manifold are orbitally asymptotically unstable (stable).(3)
WhenT2>0T2<0, the period of bifurcating periodic solutions increases (decreases).


## 5. Numerical Simulations

In the following, we will carry out numerical simulations of three impacts of diffusion, time delay, and harvesting effort to illustrate our theoretical findings in the previous sections.


### 5.1. The Impact of Diffusion

We confirm that system (13) has a unique equilibrium E∗=0.1877,0.5158 with the parameters α=0.6,b=0.4,h=0.375,f=0.5,E=0.8,p=0.5,θ=0.3. By a simple verification, we can see that A2<0 is satisfied. By Theorem 1, we deduce that system (13) is locally asymptotically stable at E∗. Meanwhile, it is revealed from Theorem 2 that system (10) is still stable at E∗, which implies Turing instability will not occur in system (10) (see Figures 1–3).Figure 1
The graph of each node and spacex in system (10) when time t=100 and D1,D2>0.Figure 2
Spatiotemporal evolution behaviors of prey and predator of system (10) with D1=0.2,D2=2 and the initial condition u0,x=0.187+0.1cosx,v0,x=0.515+0.3cosx.Figure 3
Spatiotemporal evolution behaviors of prey and predator of system (10) with D1=0.5,D2=2 and the initial condition u0,x=0.187+0.1cosx,v0,x=0.515+0.4cosx.Next, we consider system (15) and set α=0.5,e1=0.2,e2=0.32,h=0.3,f=0.6,E=0.6,p=0.52,b=0.15,θ=0.36. When D1=D2=0, system (15) can be transformed into the nondiffusion system (18). By calculation, we get E∗=0.1196,1.5350, and H1, H2 and Y2<0 are satisfied. Thus, from Lemma 5, system (18) is asymptotically stable at E∗u∗,v∗=0.1196,1.5350. Then we set D1=0.1,D1=0.8 and other variables are the same as above. Clearly, H1−H3 and Y2<0 are reached. It follows from Theorem 4 that system (15) is still stable at E∗u∗,v∗=0.1196,1.5350 (see Figure 4). When D1=0.02,D2=5 and other variables are the same as above. By a simple calculation, we can see that Y2<0, H1, H2, H4, and H5 are satisfied. By Theorem 5, the diffusion-driven instability occurs for system (15) at E∗u∗,v∗, which is shown in Figures 5 and 6.Figure 4
Spatiotemporal evolution behaviors of prey and predator of system (15) with D1=0.1,D2=0.8 and the initial condition u0,x=0.119+0.1cosx,v0,x=1.535+0.5cosx.Figure 5
Spatiotemporal evolution behaviors of prey and predator of system (15) with D1=0.02,D2=5 and the initial condition u0,x=0.119+0.1cosx,v0,x=1.535+0.5cosx.Figure 6
Turing patterns of prey and predator when Turing instability occurs.


### 5.2. The Impact of Time Delay

For system (9), we take α=0.7,b=0.4,h=0.25,f=0.5,E=0.8,p=0.45,θ=0.6,D1=0.1,D2=0.1. By calculation, we get E∗u∗,v∗=0.05813,0.7622. In the meantime, A2<0 and H6 are satisfied. By Theorem 6, system (9) is asymptotically stable at E∗u∗,v∗=0.05813,0.7622 for all τ≥0, as displayed in Figures 7 and 8.Figure 7
The evolution process of nodeut and vt in system (9) when D1=0.1,D2=0.1, τ=20.Figure 8
The evolution process of nodeut and vt in system (9) when D1=0.1,D2=0.1, τ=50.Next, we setα=0.6,b=0.4,h=0.134,f=0.8,E=1.5,p=0.3,θ=0.6,D1=0.1,D2=0.1 for system (9), which has a unique positive equilibrium E∗u∗,v∗=0.1549,0.6823. It is verified that A2<0 and H7 are satisfied. Meanwhile, we get ω0=0.0975,τ00=26.0037, λ′τ00=0.021−0.0307i, implying that the cross-section condition holds. By Theorem 9, system (9) is locally asymptotically stable at E∗u∗,v∗=0.1549,0.6823 when τ∈0,τ00 (see Figures 9 and 10) and the Hopf bifurcation occurs when τ crosses the critical value τ00 (see Figures 11 and 12).Figure 9
Phase diagram of system (9) at E∗u∗,v∗=0.1549,0.6823 when D1=0.1,D2=0.1 and τ=23<τ00=26.0037.Figure 10
The evolution process of nodeut and vt in system (9) when D1=0.1,D2=0.1 and τ=23<τ00=26.0037.Figure 11
Phase diagram of system (9) at E∗u∗,v∗=0.1549,0.6823 when D1=0.1,D2=0.1 and τ=26.2>τ00=26.0037.Figure 12
The evolution process of nodeut and vt in system (9) when D1=0.1,D2=0.1 and τ=26.2>τ00=26.0037.Furthermore, from formula (125), we obtain C10=−8.36984−2.3355i,κ2=398.56381>0,β2=−16.73968<0,T2=2.3552>0. It follows from Theorem 10 that the direction of the bifurcation in system (9) at τ00 is forward, the bifurcating period solutions are stable, and the period of bifurcating periodic solutions is increasing.


### 5.3. The Impact of Harvesting Effort

Considering the impact of harvesting efforts, we fixα=0.6,b=0.4,h=0.134,f=0.8,p=0.3,θ=0.6,D1=0.1,D2=0.1 and let E vary in 0.5,5. The stability and instability regions for system (9) are depicted by mapping the nonlinear harvesting E to the critical value τ00 in Figure 13. Obviously, the critical value τ00 of bifurcation increases with the increase of E, that is, the stable region of system (9) is expanded with the increase of E.Figure 13
Stability region exploring the dynamics of system (9) in (E, τ00) parameter space.


## 5.1. The Impact of Diffusion

We confirm that system (13) has a unique equilibrium E∗=0.1877,0.5158 with the parameters α=0.6,b=0.4,h=0.375,f=0.5,E=0.8,p=0.5,θ=0.3. By a simple verification, we can see that A2<0 is satisfied. By Theorem 1, we deduce that system (13) is locally asymptotically stable at E∗. Meanwhile, it is revealed from Theorem 2 that system (10) is still stable at E∗, which implies Turing instability will not occur in system (10) (see Figures 1–3).Figure 1
The graph of each node and spacex in system (10) when time t=100 and D1,D2>0.Figure 2
Spatiotemporal evolution behaviors of prey and predator of system (10) with D1=0.2,D2=2 and the initial condition u0,x=0.187+0.1cosx,v0,x=0.515+0.3cosx.Figure 3
Spatiotemporal evolution behaviors of prey and predator of system (10) with D1=0.5,D2=2 and the initial condition u0,x=0.187+0.1cosx,v0,x=0.515+0.4cosx.Next, we consider system (15) and set α=0.5,e1=0.2,e2=0.32,h=0.3,f=0.6,E=0.6,p=0.52,b=0.15,θ=0.36. When D1=D2=0, system (15) can be transformed into the nondiffusion system (18). By calculation, we get E∗=0.1196,1.5350, and H1, H2 and Y2<0 are satisfied. Thus, from Lemma 5, system (18) is asymptotically stable at E∗u∗,v∗=0.1196,1.5350. Then we set D1=0.1,D1=0.8 and other variables are the same as above. Clearly, H1−H3 and Y2<0 are reached. It follows from Theorem 4 that system (15) is still stable at E∗u∗,v∗=0.1196,1.5350 (see Figure 4). When D1=0.02,D2=5 and other variables are the same as above. By a simple calculation, we can see that Y2<0, H1, H2, H4, and H5 are satisfied. By Theorem 5, the diffusion-driven instability occurs for system (15) at E∗u∗,v∗, which is shown in Figures 5 and 6.Figure 4
Spatiotemporal evolution behaviors of prey and predator of system (15) with D1=0.1,D2=0.8 and the initial condition u0,x=0.119+0.1cosx,v0,x=1.535+0.5cosx.Figure 5
Spatiotemporal evolution behaviors of prey and predator of system (15) with D1=0.02,D2=5 and the initial condition u0,x=0.119+0.1cosx,v0,x=1.535+0.5cosx.Figure 6
Turing patterns of prey and predator when Turing instability occurs.


## 5.2. The Impact of Time Delay

For system (9), we take α=0.7,b=0.4,h=0.25,f=0.5,E=0.8,p=0.45,θ=0.6,D1=0.1,D2=0.1. By calculation, we get E∗u∗,v∗=0.05813,0.7622. In the meantime, A2<0 and H6 are satisfied. By Theorem 6, system (9) is asymptotically stable at E∗u∗,v∗=0.05813,0.7622 for all τ≥0, as displayed in Figures 7 and 8.Figure 7
The evolution process of nodeut and vt in system (9) when D1=0.1,D2=0.1, τ=20.Figure 8
The evolution process of nodeut and vt in system (9) when D1=0.1,D2=0.1, τ=50.Next, we setα=0.6,b=0.4,h=0.134,f=0.8,E=1.5,p=0.3,θ=0.6,D1=0.1,D2=0.1 for system (9), which has a unique positive equilibrium E∗u∗,v∗=0.1549,0.6823. It is verified that A2<0 and H7 are satisfied. Meanwhile, we get ω0=0.0975,τ00=26.0037, λ′τ00=0.021−0.0307i, implying that the cross-section condition holds. By Theorem 9, system (9) is locally asymptotically stable at E∗u∗,v∗=0.1549,0.6823 when τ∈0,τ00 (see Figures 9 and 10) and the Hopf bifurcation occurs when τ crosses the critical value τ00 (see Figures 11 and 12).Figure 9
Phase diagram of system (9) at E∗u∗,v∗=0.1549,0.6823 when D1=0.1,D2=0.1 and τ=23<τ00=26.0037.Figure 10
The evolution process of nodeut and vt in system (9) when D1=0.1,D2=0.1 and τ=23<τ00=26.0037.Figure 11
Phase diagram of system (9) at E∗u∗,v∗=0.1549,0.6823 when D1=0.1,D2=0.1 and τ=26.2>τ00=26.0037.Figure 12
The evolution process of nodeut and vt in system (9) when D1=0.1,D2=0.1 and τ=26.2>τ00=26.0037.Furthermore, from formula (125), we obtain C10=−8.36984−2.3355i,κ2=398.56381>0,β2=−16.73968<0,T2=2.3552>0. It follows from Theorem 10 that the direction of the bifurcation in system (9) at τ00 is forward, the bifurcating period solutions are stable, and the period of bifurcating periodic solutions is increasing.


## 5.3. The Impact of Harvesting Effort

Considering the impact of harvesting efforts, we fixα=0.6,b=0.4,h=0.134,f=0.8,p=0.3,θ=0.6,D1=0.1,D2=0.1 and let E vary in 0.5,5. The stability and instability regions for system (9) are depicted by mapping the nonlinear harvesting E to the critical value τ00 in Figure 13. Obviously, the critical value τ00 of bifurcation increases with the increase of E, that is, the stable region of system (9) is expanded with the increase of E.Figure 13
Stability region exploring the dynamics of system (9) in (E, τ00) parameter space.


## 6. Conclusion

We investigate the spatiotemporal dynamic evolution of a time-delay ecological competition system with food restriction and diffusion terms under Neumann boundary conditions. The conditions of asymptotic stability and Turing instability at the positive equilibrium of delay-free systems are obtained under different functional response functions. Compared with the classical population growth model described by the logistic equation, the model with food restriction studied in this paper is more in line with actual biological competition systems. The results show that the diffusion phenomenon caused by the change of population position in space will seriously affect the stability of biological competition systems, eventually resulting in the appearance of Turing instability. Then, by selecting the time delay as the bifurcation parameter, we reveal that the delay can cause very complex dynamic phenomena. When the delay is less than the bifurcation critical value, the system maintains asymptotic stability at the positive equilibrium point, while the system becomes unstable and produces a Hopf bifurcation when the delay is greater than the bifurcation critical value. Meanwhile, by using the central manifold method, we derive the conditions for determining the bifurcation direction and the stability of the bifurcation periodic solution. The harvesting effort has a major influence on the stability and Hopf bifurcation. As the harvesting effort increases in the appropriate range, the bifurcation critical value increases; that is, the stable region of the system is expanded. Therefore, we can conclude that the appropriate harvesting of biological populations and the rational development of biological resources can not only meet human economic needs but also play a beneficial role in biological systems.



---

*Source: 2823303-2022-08-03.xml*