# Stability and Bifurcation for a Delayed Diffusive Two-Zooplankton One-Phytoplankton Model with Two Different Functions

Authors: Xin-You Meng; Li Xiao
Journal: Complexity
 (2021)

---

## Abstract

In this paper, a diffusion two-phytoplankton one-zooplankton model with time delay, Beddington–DeAnglis functional response, and Holling II functional response is proposed. First, the existence and local stability of all equilibria of such model are studied. Then, the existence of Hopf bifurcation of the corresponding model without diffusion is given by taking time delay as the bifurcation parameter. Next, the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions are investigated by using the normal form theory and center manifold theorem. Furthermore, due to the local bifurcation theory of partial functional differential equations, Hopf bifurcation of the model is investigated by considering time delay as the bifurcation parameter. The explicit formulas to determine the properties of Hopf bifurcation are given by the method of the normal form theory and center manifold theorem for partial functional differential equations. Finally, some numerical simulations are performed to check out our theoretical results.


---

## Body

## 1. Introduction

Generally speaking, plankton is composed of phytoplankton and zooplankton. Phytoplankton is a main producer in the marine ecosystem and the basic link of the food chain. Zooplankton is multiplying, large, and widely distributed. The zooplankton controls the number of phytoplankton by predation and plays an important regulatory role in marine ecosystems.Plankton systems play an important role in the marine ecosystem, which have always attracted attention of ecologists and mathematicians [1–3]. Many researchers have established a lot of mathematical models to study the dynamic behaviors of plankton systems [4–10]. Phytoplankton can be divided into nontoxic phytoplankton and toxic phytoplankton (TPP). Chattopadhyay et al. [11, 12] showed that toxic phytoplankton and the toxins it released affect the growth rate of zooplankton populations and the interaction between phytoplankton and zooplankton. Chattopadhyay et al. [13] proposed a delayed differential equation for toxins and pointed out that toxic phytoplankton or toxic substances may be the key reason for the termination of plankton blooms. Similar results can be found in references [14–19]. In their work [15], Saha and Bandyopadhyay considered the following TPP-zooplankton model:(1)dPtdt=rPt1−PtK−μPtZtα+Pt,dZtdt=β1PtZtα+Pt−δZt−ρPt−τZtα+Pt−τ,where Pt and Zt represent the density of TPP and zooplankton, respectively, r and K denote the intrinsic growth rate and environmental capacity of TPP, μ is the maximum predation rate of zooplankton, β1 is the ratio of biomass consumed per zooplankton for the production of new zooplankton 0<β1<1, δ is the natural death rate of zooplankton, α denotes the half saturation constant, ρ is the ratio of toxic substances produced per unit biomass of TPP, and τ represents the discrete time required for maturation of TPP to release toxic substances. They [15] found that system (1) is locally asymptotically stable at the positive equilibrium when τ is less than its threshold; otherwise, system (1) undergoes Hopf bifurcation.There are many types of zooplankton in marine ecology and freshwater biological systems. Each zooplankton has its own living habit. There is an inseparable relationship among them. Therefore, the impact of different species on the ecosystem is the concern for most scholars. Some researchers considered two types of zooplankton [20–24]. Lv et al. [22] investigated the complex dynamics of the two zooplankton-phytoplankton model with delay:(2)dPtdt=rPt1−PtK−μ1PtZ1tα1+Pt−μ2PtZ2tα2+Pt,dZ1tdt=β1PtZ1tα1+Pt−ρ1Pt−τZ1tα1+Pt−τ−d1Z1t−g1Z12t,dZ2tdt=β2PtZ2tα2+Pt−ρ2Pt−τZ2tα2+Pt−τ−d2Z2t−g2Z22t,where Pt, Z1t, and Z2t represent the density of TPP, zooplankton 1, and zooplankton 2 at time t, respectively. The parameter μi is the maximum uptake rate of zooplankton ii=1and2. βi denotes the ratio of biomass conversion of zooplankton ii=1and2. ρi is the ratio of toxic substances produced by per unit biomass of TPP for the zooplankton ii=1and2. αi represents the half saturation constant of zooplankton ii=1and2. gi is the interspecific competition coefficient of zooplankton ii=1and2. di denotes the natural mortality of corresponding zooplankton ii=1and2. The biological significance of other parameters is the same as system (1). All parameters are positive constants. System (2) shows that the intrinsic characteristics of toxins, such as the release rate of toxicity and the delay of toxin production, does not irreversibly change the stability of this system. It is finally concluded that toxin-producing phytoplankton may be used as a biocontrol agent for the harmful algal bloom problems.In nature, species interact with each other to produce different biological effects. Each effect can vividly describe the relationship between species. Beddington [25] and DeAngelis et al. [26] proposed Beddington–DeAnglis functional response (B-D function), which is similar to the functional response of the famous Holling type II. The B-D function can better reflect the prey-dependent and the mutual interference between predators [27–33]. Cantrell and Cosner [27] investigated the predator-prey model with the B-D function and showed that the degree of mutual interference between predators affects the position and stability of the equilibria. Meng and Wang [31] established a delayed diffusive model with the B-D function and discussed the existence of Hopf bifurcation and its properties.From a biological point of view, individual organisms are distributed in space, usually interacting with the physical environment and other organisms in their spatial neighborhood. Generally speaking, self-diffusion refers to the transmission from an area of higher concentration to a low concentration area through random motion. The diffusion of individuals may be related to other things, such as finding food and escaping from high risk of infection. Hardy [34] pointed out that the spatial distribution of marine plankton is not uniform. The cooperativity of plankton is weak, but plankton can move freely under the influence of ocean currents and monsoons. This spatial diffusion is subject to Fick’s law. Therefore, the influence of spatial diffusion on the phytoplankton-plankton model has been paid more attention by many scholars [35–45]. In the work by Jia et al. [46], a three-component plankton model with spatial diffusion and time delay is proposed, which describes the relationship between a zooplankton and two phytoplankton. Rao [47] studied the complex dynamics of a spatial toxic-phytoplankton-zooplankton model with the Holling II function and showed that the interaction between toxic phytoplankton and zooplankton in the marine environment may be partially driven by diffusivity or environmental carrying capacity.Encouraged by the above work, we will consider the Beddington–DeAngelis functional response and the effect of self-diffusion for phytoplankton and zooplankton into system (2) in this paper. Thus, a diffusion toxic phytoplankton-zooplankton model with time delay and two kinds of functional responses is given as follows:(3)∂Pt,x∂t=d1ΔPt,x+rPt,x1−Pt,xK−μ1Pt,xZ1t,xα1+Pt,x+sZ1t,x−μ2Pt,xZ2t,xα2+Pt,x,t>0andx∈0,π,∂Z1t,x∂t=d2ΔZ1t,x+β1Pt,xZ1t,xα1+Pt,x+sZ1t,x−ρ1Pt−τ,xZ1t,xα1+Pt−τ,x+sZ1t,x−g1Z1t,x,t>0andx∈0,π,∂Z2t,x∂t=d3ΔZ2t,x+β2Pt,xZ2t,xα2+Pt,x−ρ2Pt−τ,xZ2t,xα2+Pt−τ,x−g2Z2t,x,t>0andx∈0,π,∂Pt,x∂x=∂Z1t,x∂x=∂Z2t,x∂x=0,t>0,x=0andπ,andPt,x=φ1t,x≥0,Z1t,x=φ2t,x≥0,Z2t,x=φ3t,x≥0,t∈−τ,0andx∈0,π,where Pt,x,Z1t,x, and Z2t,x represent the densities of the toxic-phytoplankton and two kinds of zooplankton at location x and time t, respectively. Δ is the Laplace operator, and the homogeneous Neumann boundary condition means that no plankton species enters or leaves this area. d1,d2, and d3 are the diffusion rates of the three species, gi is expressed as natural mortality of zooplankton ii=1and2, and s is the measure of the degree of mutual interference among zooplankton Z1. The remaining parameters have the same biological meanings as model (2). All parameters are positive constants.The remain parts of this paper are organized as follows. The existence and stability of equilibria are discussed in Section2. The existence and properties of Hopf bifurcation of model (2) with time delay by the normal form theory and the center manifold theorem are given in Section 3. In Section 4, we analyze the Hopf bifurcation of model (3) at the coexistence equilibrium and give the explicit formulas to determine the properties of Hopf bifurcation by using the normal form theory and center manifold theorem for partial functional differential equations (PFDEs). In Section 5, numerical simulations are given to verify the theoretical results. A discussion is given to conclude this work in Section 6.


## 2. Existence and Stability of Equilibria of ODE

In this section, we investigate the dynamic behavior of system (3) without delay and diffusion. The corresponding ordinary differential model of system (3) is as follows:(4)dPtdt=rPt1−PtK−μ1PtZ1tα1+Pt+sZ1t−μ2PtZ2tα2+Pt,dZ1tdt=β1PtZ1tα1+Pt+sZ1t−ρ1PtZ1tα1+Pt+sZ1t−g1Z1t,dZ2tdt=β2PtZ2tα2+Pt−ρ2PtZ2tα2+Pt−g2Z2t.The Jacobian matrix of system (4) is(5)J=Kr−2rPK−μ1Z1α1+sZ1α1+P+sZ12−μ2Z2α2α2+P2−μ1Pα1+Pα1+P+sZ12−μ2Pα2+Pβ1−ρ1Z1α1+sZ1α1+P+sZ12β1−ρ1Pα1+Pα1+P+sZ12−g10β2−ρ2Z2α2α2+P20β2−ρ2Pα2+P−g2.


### 2.1. The Trivial Equilibriums

It is obvious that system (4) has the trivial equilibrium E00,0,0 and the semitrivial equilibrium E1K,0,0.The characteristic equation of system (4) at equilibrium E0 is(6)λ−rλ+g1λ+g2=0.Therefore, we can obtain three eigenvalues:(7)λ1=r>0,λ2=−g1<0,λ3=−g2<0.Thus, the trivial equilibriumE0 is unstable.The local stability of system (4) at the semitrivial equilibrium E1K,0,0 is as follows.We can obtain that the characteristic equation of system (4) at the semitrivial equilibrium E1 is as follows:(8)λ+rλ−β1−ρ1Kα1+K−g1λ−β2−ρ2Kα2+K−g2=0.From equation (8), we can get three eigenvalues:(9)λ1=−r<0,λ2=β1−ρ1Kα1+K−g1,λ3=β2−ρ2Kα2+K−g2.Soλ2<0 and λ3<0 if and only if K<ming1α1/β1−ρ1−g1,g2α2/β2−ρ2−g2. Thus, the semitrivial equilibrium E1 is locally asymptotically stable; otherwise, E1 is unstable.Theorem 1.
System (4) is unstable at the trivial equilibrium E0. System (4) is locally asymptotically stable at the semitrivial equilibrium E1 if K<ming1α1/β1−ρ1−g1,g2α2/β2−ρ2−g2, but system (4) is unstable if K≥ming1α1/β1−ρ1−g1,g2α2/β2−ρ2−g2.


### 2.2. The Boundary Equilibriums

System (4) has two boundary equilibriums. From the second equation of system (4), we can obtain that Z1=1/g1sβ1−ρ1−g1P−g1α1. Then, by substituting it into the first equation of system (4), we can obtain that the following quadratic equation(10)M1P2+M2P+M3=0,where M1=rsβ1−ρ1,M2=Kμ1β1−ρ1−g1−rsKβ1−ρ1,andM3=−μ1Kg1α1.If the assumptionβ1−ρ1−g1>0, then M1>0 and Δ1=M22−4M1M3>0; thus, equation (10) has at least one positive root, defined as P¯=−M2+M22−4M1M3/2M1. From the biological meaning, Z¯1=1/g1sβ1−ρ1−g1P¯−g1α1 is nonnegative when β1−ρ1−g1P¯>g1α1. If the assumption H1:β1−ρ1−g1>0,β1−ρ1−g1P¯>g1α1 holds, then system (4) has the boundary equilibrium E10P¯,Z¯1,0.Now, under the assumptionH1, we prove the stability of the boundary equilibrium E10. Linearizing system (4) at E10, we get the characteristic equation:(11)λ−A33λ2−A11+A22λ+A11A22−A12A21=0,where(12)A11=μ1g12Z¯1β1−ρ12P¯−rP¯K,A22=g1g1α1+P¯g1+ρ1−β1β1−ρ1P¯,A13=−μ2P¯α2+P¯2,A21=β1−ρ1Z¯1α1+sZ¯1α1+P¯+sZ¯12,A12=−μ1P¯α1+P¯α1+P¯+sZ¯12,A33=β2−ρ2−g2P¯−g2α2α2+P¯.Furtherλ1=A33, λ2 and λ3 are the roots of the quadratic equation:(13)λ2+M4λ+M5=0,where M4=−A11+A22andM5=A11A22−A12A21=0. Let Δ2=M42−4M5; we assume that H2: Δ2=M42−4M5≥0, μ1Z¯1/α1+P¯+sZ¯12<r/K, and P¯<g2α2/β2−ρ2−g2. Under the assumption H2, A33<0,M4>0, and M5>0. Thus, equation (13) has two negative roots. Then, we have the following conclusion.Theorem 2.
The boundary equilibriumE10 exists if H1 holds. Furthermore, E10 of system (4) is locally asymptotically stable if assumptions H1 and H2 hold.Similarly, the other boundary equilibriumE01P˜,0,Z˜2 exists when the assumption H3: β2−ρ2−g2>0andKβ2−ρ2−g2>g2α2 holds, where P˜=g2α2/β2−ρ2−g2 and Z˜2=rα2Kβ2−ρ2−g2−g2α2β2−ρ2/μ2Kβ2−ρ2−g22.Next, we investigate stability of the boundary equilibriumE01. The method is the same to that of the boundary equilibrium E10.The characteristic equation of system (4) at the boundary equilibrium E01 is(14)λ−B22λ2−B11λ−B13B31=0,where(15)B11=μ2P˜Z˜2α2+P˜2−rK,B22=β1−ρ1−g1P˜−g1α1α1+P˜,B13=−μ2P˜α2+P˜<0,B31=β2−ρ2α2Z˜2α1+P˜2>0.Ifλ1=B22, then λ2 and λ3 are the roots of the following equation:(16)λ2−B11λ−B13B31=0.LetΔ3=B112+4B13B31; we assume that H4: Δ3>0 and β2−ρ2−g2β2−ρ2−g2K+g22α2>0. Under the assumption H4, λ1=B22<0, B11<0, and B13B31<0; thus, equation (16) has two real negative roots. Then, we have the following result.Theorem 3.
If assumptionH3 is true, the boundary equilibrium E01 of system (4) exists. Further, if assumptions H3 and H4 are true, then the boundary equilibrium E01 of system (4) is locally asymptotically stable.


### 2.3. The Coexistence Equilibrium

In biology, we are interested in the existence and stability of the coexistence equilibrium. Next, we will analyze the existence of the coexistence equilibrium of system (4). System (4) has a coexistence equilibrium E∗P∗,Z1∗,Z2∗, where P∗, Z1∗, and Z2∗ are the positive solution of the following equations:(17)r1−PK−μ1Z1α1+P+sZ1−μ2Z2α2+P=0,β1−ρ1Pα1+P+sZ1−g1=0,β2−ρ2Pα2+P−g2=0.From the third equation of (17), we can get that P∗=g2α2/β2−ρ2−g2 if β2−ρ2−g2>0 holds. By substituting P∗ into the second equation of (17), we can get that Z1∗=K1/g1sβ2−ρ2−g2 if K1>0, where K1=g2α2β1−ρ1−g1−g1α1β2−ρ2−g2. By substituting P∗ and Z1∗ into the first equation of (17), we can get that Z2∗=α1+P∗K2/g2μ2α2Kβ1−ρ1β2−ρ2−g2 if r−μ1Kβ2−ρ2−g2−rg2α2>0 holds, where K2=g2α2β1−ρ1r−μ1Kβ2−ρ2−g2−rg2α2+μ1g1Kβ2−ρ2−g2g2α2+α1β2−ρ2−g2. Thus, system (4) has one positive equilibrium E∗P∗,Z1∗,Z2∗ when the condition H5: β2>ρ2+g2, K1>0, and r−μ1Kβ2−ρ2−g2−rg2α2>0 is true.From a biological point of view, the stable coexistence of multiple species is very meaningful. Now we discuss the stability of the positive equilibriumE∗. The characteristic equation of system (4) at E∗ is(18)λ3+M6λ2+M7λ+M8=0,where(19)C11=μ1Z1∗P∗α1+P∗+sZ1∗2+μ2P∗Z2∗α2+P∗2−rP∗K,C12=−μ1P∗α1+P∗α1+P∗+sZ1∗2,C13=−μ2P∗α2+P∗,C21=β1−ρ1Z1∗α1+sZ1∗α1+P∗+sZ1∗2,C22=−sβ1−ρ1P∗Z1∗α1+P∗+sZ1∗2,C31=β2−ρ2Z2∗α2α2+P∗2,M6=−C11+C22,M7=C11C22−C12C21−C13C31,M8=C13C31C22.By the Routh–Hurwitz criterion [48], E∗ is locally asymptotically stable if M6>0, M8>0, and M6M7−M8>0. Hence, we assume that H6: μ1β1−ρ12g22α22K1/g13sβ2−ρ2−g23+α22β2−ρ22K2/Kg2β1−ρ1β2−ρ2−g23<r/K. Then, we have the following conclusion.Theorem 4.
Under the assumptionH5, the coexistence equilibrium E∗ of system (4) is locally asymptotically stable if the assumption H6 holds.Now, we consider the parameters as a bifurcation parameter to study the existence of Hopf bifurcation. According to Liu [49], we can get that the following result.Theorem 5.
If the characteristic equation of the coexistence equilibrium is given by(20)λ3+M6sλ2+M7sλ+M8s=0,where M6s,Ω=M6sM7s−M8s, and M8s are smooth function of s in an open interval about s∗∈R such that(1)
M6s∗>0, Ωs∗=0, and M8s∗>0(2)
dΩ/ds|s=s∗≠0; then, a simple Hopf bifurcation occurs at s=s∗Proof. In fact, M6s=M6, Ωs=M6M7−M8, and M8s=M8. It is easy to see that there exists s∗ such that M6s∗>0, Ωs∗=0, M8s∗>0, dP∗/ds|s=s∗=0, dZ1∗/ds|s=s∗=−β1−ρ1−g1P∗−g1α1/g1s2≠0, and dZ2∗/ds|s=s∗=−g1μ1β2−ρ2/μ2g2β1−ρ1⋅dZ1∗/ds≠0. So, we can get dΩ/ds|s=s∗≠0. Thus, the Hopf bifurcation occurs at s=s∗.□


## 2.1. The Trivial Equilibriums

It is obvious that system (4) has the trivial equilibrium E00,0,0 and the semitrivial equilibrium E1K,0,0.The characteristic equation of system (4) at equilibrium E0 is(6)λ−rλ+g1λ+g2=0.Therefore, we can obtain three eigenvalues:(7)λ1=r>0,λ2=−g1<0,λ3=−g2<0.Thus, the trivial equilibriumE0 is unstable.The local stability of system (4) at the semitrivial equilibrium E1K,0,0 is as follows.We can obtain that the characteristic equation of system (4) at the semitrivial equilibrium E1 is as follows:(8)λ+rλ−β1−ρ1Kα1+K−g1λ−β2−ρ2Kα2+K−g2=0.From equation (8), we can get three eigenvalues:(9)λ1=−r<0,λ2=β1−ρ1Kα1+K−g1,λ3=β2−ρ2Kα2+K−g2.Soλ2<0 and λ3<0 if and only if K<ming1α1/β1−ρ1−g1,g2α2/β2−ρ2−g2. Thus, the semitrivial equilibrium E1 is locally asymptotically stable; otherwise, E1 is unstable.Theorem 1.
System (4) is unstable at the trivial equilibrium E0. System (4) is locally asymptotically stable at the semitrivial equilibrium E1 if K<ming1α1/β1−ρ1−g1,g2α2/β2−ρ2−g2, but system (4) is unstable if K≥ming1α1/β1−ρ1−g1,g2α2/β2−ρ2−g2.


## 2.2. The Boundary Equilibriums

System (4) has two boundary equilibriums. From the second equation of system (4), we can obtain that Z1=1/g1sβ1−ρ1−g1P−g1α1. Then, by substituting it into the first equation of system (4), we can obtain that the following quadratic equation(10)M1P2+M2P+M3=0,where M1=rsβ1−ρ1,M2=Kμ1β1−ρ1−g1−rsKβ1−ρ1,andM3=−μ1Kg1α1.If the assumptionβ1−ρ1−g1>0, then M1>0 and Δ1=M22−4M1M3>0; thus, equation (10) has at least one positive root, defined as P¯=−M2+M22−4M1M3/2M1. From the biological meaning, Z¯1=1/g1sβ1−ρ1−g1P¯−g1α1 is nonnegative when β1−ρ1−g1P¯>g1α1. If the assumption H1:β1−ρ1−g1>0,β1−ρ1−g1P¯>g1α1 holds, then system (4) has the boundary equilibrium E10P¯,Z¯1,0.Now, under the assumptionH1, we prove the stability of the boundary equilibrium E10. Linearizing system (4) at E10, we get the characteristic equation:(11)λ−A33λ2−A11+A22λ+A11A22−A12A21=0,where(12)A11=μ1g12Z¯1β1−ρ12P¯−rP¯K,A22=g1g1α1+P¯g1+ρ1−β1β1−ρ1P¯,A13=−μ2P¯α2+P¯2,A21=β1−ρ1Z¯1α1+sZ¯1α1+P¯+sZ¯12,A12=−μ1P¯α1+P¯α1+P¯+sZ¯12,A33=β2−ρ2−g2P¯−g2α2α2+P¯.Furtherλ1=A33, λ2 and λ3 are the roots of the quadratic equation:(13)λ2+M4λ+M5=0,where M4=−A11+A22andM5=A11A22−A12A21=0. Let Δ2=M42−4M5; we assume that H2: Δ2=M42−4M5≥0, μ1Z¯1/α1+P¯+sZ¯12<r/K, and P¯<g2α2/β2−ρ2−g2. Under the assumption H2, A33<0,M4>0, and M5>0. Thus, equation (13) has two negative roots. Then, we have the following conclusion.Theorem 2.
The boundary equilibriumE10 exists if H1 holds. Furthermore, E10 of system (4) is locally asymptotically stable if assumptions H1 and H2 hold.Similarly, the other boundary equilibriumE01P˜,0,Z˜2 exists when the assumption H3: β2−ρ2−g2>0andKβ2−ρ2−g2>g2α2 holds, where P˜=g2α2/β2−ρ2−g2 and Z˜2=rα2Kβ2−ρ2−g2−g2α2β2−ρ2/μ2Kβ2−ρ2−g22.Next, we investigate stability of the boundary equilibriumE01. The method is the same to that of the boundary equilibrium E10.The characteristic equation of system (4) at the boundary equilibrium E01 is(14)λ−B22λ2−B11λ−B13B31=0,where(15)B11=μ2P˜Z˜2α2+P˜2−rK,B22=β1−ρ1−g1P˜−g1α1α1+P˜,B13=−μ2P˜α2+P˜<0,B31=β2−ρ2α2Z˜2α1+P˜2>0.Ifλ1=B22, then λ2 and λ3 are the roots of the following equation:(16)λ2−B11λ−B13B31=0.LetΔ3=B112+4B13B31; we assume that H4: Δ3>0 and β2−ρ2−g2β2−ρ2−g2K+g22α2>0. Under the assumption H4, λ1=B22<0, B11<0, and B13B31<0; thus, equation (16) has two real negative roots. Then, we have the following result.Theorem 3.
If assumptionH3 is true, the boundary equilibrium E01 of system (4) exists. Further, if assumptions H3 and H4 are true, then the boundary equilibrium E01 of system (4) is locally asymptotically stable.


## 2.3. The Coexistence Equilibrium

In biology, we are interested in the existence and stability of the coexistence equilibrium. Next, we will analyze the existence of the coexistence equilibrium of system (4). System (4) has a coexistence equilibrium E∗P∗,Z1∗,Z2∗, where P∗, Z1∗, and Z2∗ are the positive solution of the following equations:(17)r1−PK−μ1Z1α1+P+sZ1−μ2Z2α2+P=0,β1−ρ1Pα1+P+sZ1−g1=0,β2−ρ2Pα2+P−g2=0.From the third equation of (17), we can get that P∗=g2α2/β2−ρ2−g2 if β2−ρ2−g2>0 holds. By substituting P∗ into the second equation of (17), we can get that Z1∗=K1/g1sβ2−ρ2−g2 if K1>0, where K1=g2α2β1−ρ1−g1−g1α1β2−ρ2−g2. By substituting P∗ and Z1∗ into the first equation of (17), we can get that Z2∗=α1+P∗K2/g2μ2α2Kβ1−ρ1β2−ρ2−g2 if r−μ1Kβ2−ρ2−g2−rg2α2>0 holds, where K2=g2α2β1−ρ1r−μ1Kβ2−ρ2−g2−rg2α2+μ1g1Kβ2−ρ2−g2g2α2+α1β2−ρ2−g2. Thus, system (4) has one positive equilibrium E∗P∗,Z1∗,Z2∗ when the condition H5: β2>ρ2+g2, K1>0, and r−μ1Kβ2−ρ2−g2−rg2α2>0 is true.From a biological point of view, the stable coexistence of multiple species is very meaningful. Now we discuss the stability of the positive equilibriumE∗. The characteristic equation of system (4) at E∗ is(18)λ3+M6λ2+M7λ+M8=0,where(19)C11=μ1Z1∗P∗α1+P∗+sZ1∗2+μ2P∗Z2∗α2+P∗2−rP∗K,C12=−μ1P∗α1+P∗α1+P∗+sZ1∗2,C13=−μ2P∗α2+P∗,C21=β1−ρ1Z1∗α1+sZ1∗α1+P∗+sZ1∗2,C22=−sβ1−ρ1P∗Z1∗α1+P∗+sZ1∗2,C31=β2−ρ2Z2∗α2α2+P∗2,M6=−C11+C22,M7=C11C22−C12C21−C13C31,M8=C13C31C22.By the Routh–Hurwitz criterion [48], E∗ is locally asymptotically stable if M6>0, M8>0, and M6M7−M8>0. Hence, we assume that H6: μ1β1−ρ12g22α22K1/g13sβ2−ρ2−g23+α22β2−ρ22K2/Kg2β1−ρ1β2−ρ2−g23<r/K. Then, we have the following conclusion.Theorem 4.
Under the assumptionH5, the coexistence equilibrium E∗ of system (4) is locally asymptotically stable if the assumption H6 holds.Now, we consider the parameters as a bifurcation parameter to study the existence of Hopf bifurcation. According to Liu [49], we can get that the following result.Theorem 5.
If the characteristic equation of the coexistence equilibrium is given by(20)λ3+M6sλ2+M7sλ+M8s=0,where M6s,Ω=M6sM7s−M8s, and M8s are smooth function of s in an open interval about s∗∈R such that(1)
M6s∗>0, Ωs∗=0, and M8s∗>0(2)
dΩ/ds|s=s∗≠0; then, a simple Hopf bifurcation occurs at s=s∗Proof. In fact, M6s=M6, Ωs=M6M7−M8, and M8s=M8. It is easy to see that there exists s∗ such that M6s∗>0, Ωs∗=0, M8s∗>0, dP∗/ds|s=s∗=0, dZ1∗/ds|s=s∗=−β1−ρ1−g1P∗−g1α1/g1s2≠0, and dZ2∗/ds|s=s∗=−g1μ1β2−ρ2/μ2g2β1−ρ1⋅dZ1∗/ds≠0. So, we can get dΩ/ds|s=s∗≠0. Thus, the Hopf bifurcation occurs at s=s∗.□


## 3. Hopf Bifurcation of System (3) without Diffusion

### 3.1. The Existence of Bifurcation

In this section, under the assumptionH5, we regard time delay τ as the bifurcation parameter to study its influence on the stability of coexistence equilibrium E∗P∗,Z1∗,Z2∗. System (4) with time delay is as follows:(21)dPdt=rP1−PK−μ1PZ1α1+P+sZ1−μ2PZ2α2+P,dZ1dt=β1PZ1α1+P+sZ1−ρ1Pt−τZ1α1+Pt−τ+sZ1−g1Z1,dZ2dt=β2PZ2α2+P−ρ2Pt−τZ2α2+Pt−τ−g2Z2.Letu1t=Pt−P∗,u2t=Z1t−Z1∗,u3t=Z2t−Z2∗,andUt=u1t,u2t,u3tT. Substituting Ut into system (21), we can get the following linearized system:(22)dUtdt=L1Ut+L2Ut−τ,where(23)L1=a11a12a13a21a220a3100,L2=000b2100b3100,where a1k=C1kk=1,2,3 and(24)a21=β1Z1∗α1+sZ1∗α1+P∗+sZ1∗2,a22=β1−ρ1P∗α1+P∗α1+P∗+sZ1∗2−g1,a31=β2Z2∗α2α2+P∗2,b21=−ρ1Z1∗α1+sZ1∗α1+P¯+sZ1∗2,b31=−ρ2Z2∗α2α2+P∗2.Therefore, the characteristic equation of system (22) is(25)λ3+D2λ2+D1λ+D0+e−λτF1λ+F0=0,where(26)D2=−a11+a22,D1=a11a22−a12a21−a13a31,D0=a13a22a31,F1=−a13b31+a12b21,F0=a13a22b31.Whenτ=0, the distribution of roots of equation (25) has been discussed above. Next, we will discuss the distribution of characteristic roots of equation (25) when τ>0. First, we suppose that iω1ω1>0 is a root of equation (25). Then, we obtain that(27)−iω13−D2ω12+iD1ω1+D0+cosω1τ−isinω1τiF1ω1−F0=0.Separating the real and imaginary parts of equation (27), we have(28)ω13−D1ω1=F1ω1cosω1τ−F0sinω1τ,D2ω12−D0=F0cosω1τ+F1ω1sinω1τ.From equation (28), we get that(29)ω16+D22−2D1ω14+D12−2D0D2−F12ω12+D02−F02=0.Letω12=y, M9=D22−2D1, M10=D12−2D0D2−F12, and M11=D02−F02. Then, equation (29) can be rewritten as(30)fy≜y3+M9y2+M10y+M11=0.To obtain the existence of Hopf bifurcation, equation (30) has at least one positive root. Due to f0=M11=a132a222a31+b31a31+b31>0, limy⟶∞fy=∞, and f′y=3y2+2M9y+M10. Hence, we assume that H7: Δ4=M92−3M10≤0; H8: Δ4=M92−3M10>0, y+=−M9+Δ4/3>0, and fy+≤0.From Lemmas 2.2 and 4.2 in [50], we obtain the following lemma.Lemma 1.
For polynomial equation (30), the following conclusions are true.(1)
If the assumptionH7 holds, equation (30) has no positive root(2)
If the assumptionH8 holds, equation (30) has at least one positive rootWithout losing of generality, we assume that equation (30) has three positive roots, denoted by y1, y2, and y3. Hence, equation (29) has three positive roots ω1k=ykk=1,2,and3. By equation (28), we can get that the critical values of time delay τ are(31)τkj=1ω1karccosF1ω1k4+F0D2−F1D1ω1k2−F0D0F12ω1k2+F02+2jπω1k,for k=1,2,and3,j∈N0. Therefore, ±iω1k is a pair of purely imaginary roots of equation (25) when τ=τkjk=1,2,and3;j=1,2,…. Denote τ0=τ00=mink=1,2,3,j∈N0τkj and ω10=ω1k|τ=τ0.Letλτ=ατ±ντ be the root of equation (25) near τ=τkj with ατkj=0 and ντkj=ω1k,k=1,2,and3andj∈N0. Then, we have the following transversality condition.Lemma 2.
Suppose thatf′yk≠0; then, the following transversality condition is satisfied: signdReλτ/dτλ=iω10≠0.Proof.
Differentiating the two sides of equation (25) with respect to τ, we have(32)dλdτ−1=3λ2+2D2λ+D1eλτ+F1λF1λ+F0−τλ.Then, we obtain(33)dReλτdτλ=iω10−1=Re3λ2+2D2λ+D1eλτλF1λ+F0+F1λF1λ+F0λ=iω10.By some calculation, we get(34)dReλτdτλ=iω10−1=3ω106+2D22−2D1ω104+D12−2D0D2−F12ω102Λ1,=ykΛ13yk2+2D22−2D1yk+D12−2D0D2−F12,=ykΛ1f′yk≠0,where Λ1=F12ω12+F02.Hence,(35)signdReλτdτλ=iω10=signdReλτdτλ=iω10−1≠0.By the theory of Hopf bifurcation, we have the following result.Theorem 6.
Assume thatH5 and H6 hold, then the following statements are true.(1)
If the hypothesisH7 holds, then the coexistence equilibrium E∗ of model (21) is locally asymptotically stable for τ≥0.(2)
If the hypothesisH8 and f′yk≠0k=1,2,and3 hold, then the coexistence equilibrium E∗ of model (21) is locally asymptotically stable for τ∈0,τ0 and unstable for τ>τ0. Moreover, model (21) undergoes Hopf bifurcation at the coexistence equilibrium E∗ when τ=τ0.


### 3.2. The Properties of Hopf Bifurcation

In the previous parts, we have discussed the existence of Hopf bifurcation of model (21) at the coexistence equilibrium E∗. Next, we will study the direction of Hopf bifurcation and the stability of bifurcated periodic solution by using the normal form theory and center manifold theory developed by Hassard et al. [51]. Let μ=τ−τ0 with μ∈R; then, μ=0 is a Hopf bifurcation value of model (21). Let u1t=Pt−P∗, u2t=Z1t−Z1∗, and u3t=Z2t−Z2∗; then, model (21) can be rewritten in the phase space C=C−τ0,0,R3:(36)U˙t=LUt+Fμ,Ut,where Ut=u1t,u2t,u3tT∈R3,L:C⟶R3,F:R×C⟶R3,φθ=φ1θ,φ2θ,φ3θ∈R3 has the following:(37)Lφ=L1φ0+L2φ−τ0−μ,Fμ,Ut=a14φ120+a15φ220+a16φ10φ20+a17φ10φ30a24φ120+a25φ220+a16φ10φ20+b24φ12−τ0−μ+b25φ1−τ0−μφ20a34φ120+a35φ10φ30+b34φ12−τ0−μ+b35φ1−τ0−μφ30,where(38)a14=−1K+μ1Z1∗α1+sZ1∗α1+P∗+sZ1∗3+μ2α2Z2∗α2+P∗3,a15=μ1sP∗α1+P∗α1+P∗+sZ1∗3,a17=−μ2α2α2+P∗2,a16=−μ1α1α1+P∗+sZ1∗+2μ1sP∗Z1∗α1+P∗+sZ1∗3,a24=−β1Z1∗α1+sZ1∗α1+P∗+sZ1∗3,a25=β1−ρ1P∗α1+P∗α1+P∗+sZ1∗3,a26=β1α1α1+P∗+sZ1∗+2β1sP∗Z1∗α1+P∗+sZ1∗3,a34=−β2α2Z2∗α2+P∗2,a35=β2α2α2+P∗2,b25=−ρ1α1α1+P∗+sZ1∗+2ρ1sP∗Z1∗α1+P∗+sZ1∗3,b24=ρ1Z1∗α1+sZ1∗α1+P∗+sZ1∗3,b35=−ρ2α2α2+P∗2,b34=ρ2α2Z2∗α2+P∗2,and L1 and L2 are given in (23).According to the Riesz representation theorem [52], there exists a function whose elements are ηθ,μ−τ0<θ<0 of bounded variance such that(39)Lφ=∫−τ00dηθ,μφθ,φ∈C−τ0,0,R3.Thus, we chose(40)ηθ,μ=L1,θ=0,0,θ∈−τ0,0,−L2,θ=−τ0.Forφ∈C1−τ0,R3, we define the operator:(41)Aμφθ=dφθdθ,θ∈−τ0,0,∫−τ00dηθ,μφθ,θ=0,Rμφθ=0,θ∈−τ0,0,Fφ,μ,θ=0.Then, system (36) is equivalent to(42)U˙t=AμUt+RμUt.Forφ∈C1−τ0,0,R3, A∗ is defined as the adjoint operator of A as(43)A∗ψξ=−dψξdξ,ξ∈0,τ0,∫−τ00dηTζ,μψ−ζ,ξ=0.The bilinear product is given by(44)ψ,φ=ψ¯0φ0−∫−τ00∫0θψ¯ξ−θdηθφξdξ,where ηθ=ηθ,0.By the discussion in Section3.1, ±iω10 are the eigenvalues of A0 and A∗0, which satisfy A0qθ=iω10qθ and A∗0=−iω10q∗ξ. Suppose that qθ=1,q2,q3Teiω10θ is the eigenvector of A0 corresponding to the eigenvalue iω10 and q∗ξ=D1,q2∗,q3∗e−iω10ξ is the eigenvector of A∗0 corresponding to the eigenvalue −iω10. By computation, we have(45)q1=a21iω10−a22,q2=iω10−a11iω10−a22−a12a21a13iω10−a22,q1∗=−a12iω10+a22,q2∗=−iω10+a11iω10+a22−a12a21+b21eiω10τ0a31+b31eiω10τ0iω10+a22.From equation (44) and q∗,q=1, we can get(46)q∗,q=q¯∗0q0−∫−τ00∫0θq¯∗ξ−θdηθqξdξ,=D¯1,q¯1∗,q¯2∗1,q1,q2T−∫−τ00∫0θ1,q¯1∗,q¯2∗1,q1,q2Teiω10θdηθdξ,=D¯1+q¯1∗q1+q¯2∗q2+τ0e−iω10τ0b21q¯1∗+b31q¯2∗,=1.Thus,(47)D¯=1+q¯1∗q1+q¯2∗q2+τ0e−iω10τ0b21q¯1∗+b31q¯2∗−1.In the remainder of this section, through using the method in [51], we can get the following relevant parameters which can determine the direction of Hopf bifurcation and the stability of the bifurcation periodic solutions:(48)g20=2D¯k11+q¯1∗k21+q¯2∗k31,g11=D¯k12+q¯1∗k22+q¯2∗k32,g02=2D¯k13+q¯1∗k23+q¯2∗k33,g21=2D¯k14+q¯1∗k24+q¯2∗k34,where(49)k11=a14+q12a15+q1a16+q2a17,k13=a14+q¯12a15+q¯1a16+q¯2a17,k12=2a14+2q1q¯1a15+q1+q¯1a16+q2+q¯2a17,k21=a24+q12a25+q1a26+e−2iω10τ0b24+q1eiω10τ0b25,k22=2a24+2q1q¯1a25+q1+q¯1a26+2b24+q1eiω10τ0+q¯1e−iω10τ0b25,k23=a24+q¯12a25+q¯1a26+e2iω10τ0b24+q¯1e−iω10τ0b25,k31=a34+q2a35+e−2iω10τ0b34+q2e−iω10τ0b35,k32=a34+q2+q¯2a35+2b34+q2eiω10τ0+q¯1e−iω10τ0b35,k33=a34+q¯2a35+e2iω10τ0b34+q¯2eiω10τ0b35,k14=a14W2010+2W1110+a15q¯1W2020+2q1W1120+a1612W2020+W1120+q1W1110+12q¯1W2010+a1712W2030+12q¯2W2010+W1130+q2W1110,k24=a24W2010+2W1110+a25q¯1W2020+2q1W1120+a2612W2020+W1120+q1W1110+12q¯1W2010+b24eiω10τ0W201−τ0+2e−iω10τ0W112−τ0+b2512q¯1W201−τ0+q1W111−τ0+e−iω10τ0W1120+12eiω10τ0W202,k34=a34W2010+2W1110+a3512W2030+12q¯2W2010+W1130+q2W1110+b34eiω10τ0W201−τ0+2e−iω10τ0W112−τ0+b35q2W111−τ0+12q¯2W201−τ0+12eiω10τ0W2030+e−iω10τ0W1130,W20θ=ig20q0ω10eiω10θ+ig¯02q¯03ω10e−iω10θ+E1e2iω10θ,W11θ=−ig11q0ω10eiω10θ+ig¯11q¯0ω10e−iω10θ+E2,where E1=E11,E12,E13T∈R3 and E2=E21,E22,E23T∈R3 can be determined by the following, respectively:(50)E1=22iω10−a11−a12−a13−a21−b21e−2iω10τ02iω10−a220−a31−b31e−2iω10τ002iω10−1k11k21k31,E2=−a11−a12−a13−a21−b21−a220−a31−b3100−1k12k22k32.So we can getg21 and compute the following values:(51)c10=i2ω10τ0g20g11−2g112−g0223+g212,ν1=−Rec10Reλ′τ0,ϱ1=2Rec10,T1=−Imc10+ν1Imλ′τ0ω10τ0,which determine the direction of Hopf bifurcation and the stability of bifurcating periodic solutions when τ=τ0. From the conclusions of Hassard et al. [51], we have the following results.Theorem 7.
SupposeH5, H6, H8, and f′yii=1,2,and3 are satisfied. For the properties of Hopf bifurcation of model (21), we have that(1)
Ifν1>0<0, then the direction of the Hopf bifurcation is forward (backward) and the bifurcating periodic solutions exist for τ>τ0τ<τ0(2)
Ifϱ1<0>0, then the bifurcating periodic solutions are orbitally stable (unstable)(3)
IfT1>0<0, then the period of bifurcation periodic solutions increase (decrease)


## 3.1. The Existence of Bifurcation

In this section, under the assumptionH5, we regard time delay τ as the bifurcation parameter to study its influence on the stability of coexistence equilibrium E∗P∗,Z1∗,Z2∗. System (4) with time delay is as follows:(21)dPdt=rP1−PK−μ1PZ1α1+P+sZ1−μ2PZ2α2+P,dZ1dt=β1PZ1α1+P+sZ1−ρ1Pt−τZ1α1+Pt−τ+sZ1−g1Z1,dZ2dt=β2PZ2α2+P−ρ2Pt−τZ2α2+Pt−τ−g2Z2.Letu1t=Pt−P∗,u2t=Z1t−Z1∗,u3t=Z2t−Z2∗,andUt=u1t,u2t,u3tT. Substituting Ut into system (21), we can get the following linearized system:(22)dUtdt=L1Ut+L2Ut−τ,where(23)L1=a11a12a13a21a220a3100,L2=000b2100b3100,where a1k=C1kk=1,2,3 and(24)a21=β1Z1∗α1+sZ1∗α1+P∗+sZ1∗2,a22=β1−ρ1P∗α1+P∗α1+P∗+sZ1∗2−g1,a31=β2Z2∗α2α2+P∗2,b21=−ρ1Z1∗α1+sZ1∗α1+P¯+sZ1∗2,b31=−ρ2Z2∗α2α2+P∗2.Therefore, the characteristic equation of system (22) is(25)λ3+D2λ2+D1λ+D0+e−λτF1λ+F0=0,where(26)D2=−a11+a22,D1=a11a22−a12a21−a13a31,D0=a13a22a31,F1=−a13b31+a12b21,F0=a13a22b31.Whenτ=0, the distribution of roots of equation (25) has been discussed above. Next, we will discuss the distribution of characteristic roots of equation (25) when τ>0. First, we suppose that iω1ω1>0 is a root of equation (25). Then, we obtain that(27)−iω13−D2ω12+iD1ω1+D0+cosω1τ−isinω1τiF1ω1−F0=0.Separating the real and imaginary parts of equation (27), we have(28)ω13−D1ω1=F1ω1cosω1τ−F0sinω1τ,D2ω12−D0=F0cosω1τ+F1ω1sinω1τ.From equation (28), we get that(29)ω16+D22−2D1ω14+D12−2D0D2−F12ω12+D02−F02=0.Letω12=y, M9=D22−2D1, M10=D12−2D0D2−F12, and M11=D02−F02. Then, equation (29) can be rewritten as(30)fy≜y3+M9y2+M10y+M11=0.To obtain the existence of Hopf bifurcation, equation (30) has at least one positive root. Due to f0=M11=a132a222a31+b31a31+b31>0, limy⟶∞fy=∞, and f′y=3y2+2M9y+M10. Hence, we assume that H7: Δ4=M92−3M10≤0; H8: Δ4=M92−3M10>0, y+=−M9+Δ4/3>0, and fy+≤0.From Lemmas 2.2 and 4.2 in [50], we obtain the following lemma.Lemma 1.
For polynomial equation (30), the following conclusions are true.(1)
If the assumptionH7 holds, equation (30) has no positive root(2)
If the assumptionH8 holds, equation (30) has at least one positive rootWithout losing of generality, we assume that equation (30) has three positive roots, denoted by y1, y2, and y3. Hence, equation (29) has three positive roots ω1k=ykk=1,2,and3. By equation (28), we can get that the critical values of time delay τ are(31)τkj=1ω1karccosF1ω1k4+F0D2−F1D1ω1k2−F0D0F12ω1k2+F02+2jπω1k,for k=1,2,and3,j∈N0. Therefore, ±iω1k is a pair of purely imaginary roots of equation (25) when τ=τkjk=1,2,and3;j=1,2,…. Denote τ0=τ00=mink=1,2,3,j∈N0τkj and ω10=ω1k|τ=τ0.Letλτ=ατ±ντ be the root of equation (25) near τ=τkj with ατkj=0 and ντkj=ω1k,k=1,2,and3andj∈N0. Then, we have the following transversality condition.Lemma 2.
Suppose thatf′yk≠0; then, the following transversality condition is satisfied: signdReλτ/dτλ=iω10≠0.Proof.
Differentiating the two sides of equation (25) with respect to τ, we have(32)dλdτ−1=3λ2+2D2λ+D1eλτ+F1λF1λ+F0−τλ.Then, we obtain(33)dReλτdτλ=iω10−1=Re3λ2+2D2λ+D1eλτλF1λ+F0+F1λF1λ+F0λ=iω10.By some calculation, we get(34)dReλτdτλ=iω10−1=3ω106+2D22−2D1ω104+D12−2D0D2−F12ω102Λ1,=ykΛ13yk2+2D22−2D1yk+D12−2D0D2−F12,=ykΛ1f′yk≠0,where Λ1=F12ω12+F02.Hence,(35)signdReλτdτλ=iω10=signdReλτdτλ=iω10−1≠0.By the theory of Hopf bifurcation, we have the following result.Theorem 6.
Assume thatH5 and H6 hold, then the following statements are true.(1)
If the hypothesisH7 holds, then the coexistence equilibrium E∗ of model (21) is locally asymptotically stable for τ≥0.(2)
If the hypothesisH8 and f′yk≠0k=1,2,and3 hold, then the coexistence equilibrium E∗ of model (21) is locally asymptotically stable for τ∈0,τ0 and unstable for τ>τ0. Moreover, model (21) undergoes Hopf bifurcation at the coexistence equilibrium E∗ when τ=τ0.


## 3.2. The Properties of Hopf Bifurcation

In the previous parts, we have discussed the existence of Hopf bifurcation of model (21) at the coexistence equilibrium E∗. Next, we will study the direction of Hopf bifurcation and the stability of bifurcated periodic solution by using the normal form theory and center manifold theory developed by Hassard et al. [51]. Let μ=τ−τ0 with μ∈R; then, μ=0 is a Hopf bifurcation value of model (21). Let u1t=Pt−P∗, u2t=Z1t−Z1∗, and u3t=Z2t−Z2∗; then, model (21) can be rewritten in the phase space C=C−τ0,0,R3:(36)U˙t=LUt+Fμ,Ut,where Ut=u1t,u2t,u3tT∈R3,L:C⟶R3,F:R×C⟶R3,φθ=φ1θ,φ2θ,φ3θ∈R3 has the following:(37)Lφ=L1φ0+L2φ−τ0−μ,Fμ,Ut=a14φ120+a15φ220+a16φ10φ20+a17φ10φ30a24φ120+a25φ220+a16φ10φ20+b24φ12−τ0−μ+b25φ1−τ0−μφ20a34φ120+a35φ10φ30+b34φ12−τ0−μ+b35φ1−τ0−μφ30,where(38)a14=−1K+μ1Z1∗α1+sZ1∗α1+P∗+sZ1∗3+μ2α2Z2∗α2+P∗3,a15=μ1sP∗α1+P∗α1+P∗+sZ1∗3,a17=−μ2α2α2+P∗2,a16=−μ1α1α1+P∗+sZ1∗+2μ1sP∗Z1∗α1+P∗+sZ1∗3,a24=−β1Z1∗α1+sZ1∗α1+P∗+sZ1∗3,a25=β1−ρ1P∗α1+P∗α1+P∗+sZ1∗3,a26=β1α1α1+P∗+sZ1∗+2β1sP∗Z1∗α1+P∗+sZ1∗3,a34=−β2α2Z2∗α2+P∗2,a35=β2α2α2+P∗2,b25=−ρ1α1α1+P∗+sZ1∗+2ρ1sP∗Z1∗α1+P∗+sZ1∗3,b24=ρ1Z1∗α1+sZ1∗α1+P∗+sZ1∗3,b35=−ρ2α2α2+P∗2,b34=ρ2α2Z2∗α2+P∗2,and L1 and L2 are given in (23).According to the Riesz representation theorem [52], there exists a function whose elements are ηθ,μ−τ0<θ<0 of bounded variance such that(39)Lφ=∫−τ00dηθ,μφθ,φ∈C−τ0,0,R3.Thus, we chose(40)ηθ,μ=L1,θ=0,0,θ∈−τ0,0,−L2,θ=−τ0.Forφ∈C1−τ0,R3, we define the operator:(41)Aμφθ=dφθdθ,θ∈−τ0,0,∫−τ00dηθ,μφθ,θ=0,Rμφθ=0,θ∈−τ0,0,Fφ,μ,θ=0.Then, system (36) is equivalent to(42)U˙t=AμUt+RμUt.Forφ∈C1−τ0,0,R3, A∗ is defined as the adjoint operator of A as(43)A∗ψξ=−dψξdξ,ξ∈0,τ0,∫−τ00dηTζ,μψ−ζ,ξ=0.The bilinear product is given by(44)ψ,φ=ψ¯0φ0−∫−τ00∫0θψ¯ξ−θdηθφξdξ,where ηθ=ηθ,0.By the discussion in Section3.1, ±iω10 are the eigenvalues of A0 and A∗0, which satisfy A0qθ=iω10qθ and A∗0=−iω10q∗ξ. Suppose that qθ=1,q2,q3Teiω10θ is the eigenvector of A0 corresponding to the eigenvalue iω10 and q∗ξ=D1,q2∗,q3∗e−iω10ξ is the eigenvector of A∗0 corresponding to the eigenvalue −iω10. By computation, we have(45)q1=a21iω10−a22,q2=iω10−a11iω10−a22−a12a21a13iω10−a22,q1∗=−a12iω10+a22,q2∗=−iω10+a11iω10+a22−a12a21+b21eiω10τ0a31+b31eiω10τ0iω10+a22.From equation (44) and q∗,q=1, we can get(46)q∗,q=q¯∗0q0−∫−τ00∫0θq¯∗ξ−θdηθqξdξ,=D¯1,q¯1∗,q¯2∗1,q1,q2T−∫−τ00∫0θ1,q¯1∗,q¯2∗1,q1,q2Teiω10θdηθdξ,=D¯1+q¯1∗q1+q¯2∗q2+τ0e−iω10τ0b21q¯1∗+b31q¯2∗,=1.Thus,(47)D¯=1+q¯1∗q1+q¯2∗q2+τ0e−iω10τ0b21q¯1∗+b31q¯2∗−1.In the remainder of this section, through using the method in [51], we can get the following relevant parameters which can determine the direction of Hopf bifurcation and the stability of the bifurcation periodic solutions:(48)g20=2D¯k11+q¯1∗k21+q¯2∗k31,g11=D¯k12+q¯1∗k22+q¯2∗k32,g02=2D¯k13+q¯1∗k23+q¯2∗k33,g21=2D¯k14+q¯1∗k24+q¯2∗k34,where(49)k11=a14+q12a15+q1a16+q2a17,k13=a14+q¯12a15+q¯1a16+q¯2a17,k12=2a14+2q1q¯1a15+q1+q¯1a16+q2+q¯2a17,k21=a24+q12a25+q1a26+e−2iω10τ0b24+q1eiω10τ0b25,k22=2a24+2q1q¯1a25+q1+q¯1a26+2b24+q1eiω10τ0+q¯1e−iω10τ0b25,k23=a24+q¯12a25+q¯1a26+e2iω10τ0b24+q¯1e−iω10τ0b25,k31=a34+q2a35+e−2iω10τ0b34+q2e−iω10τ0b35,k32=a34+q2+q¯2a35+2b34+q2eiω10τ0+q¯1e−iω10τ0b35,k33=a34+q¯2a35+e2iω10τ0b34+q¯2eiω10τ0b35,k14=a14W2010+2W1110+a15q¯1W2020+2q1W1120+a1612W2020+W1120+q1W1110+12q¯1W2010+a1712W2030+12q¯2W2010+W1130+q2W1110,k24=a24W2010+2W1110+a25q¯1W2020+2q1W1120+a2612W2020+W1120+q1W1110+12q¯1W2010+b24eiω10τ0W201−τ0+2e−iω10τ0W112−τ0+b2512q¯1W201−τ0+q1W111−τ0+e−iω10τ0W1120+12eiω10τ0W202,k34=a34W2010+2W1110+a3512W2030+12q¯2W2010+W1130+q2W1110+b34eiω10τ0W201−τ0+2e−iω10τ0W112−τ0+b35q2W111−τ0+12q¯2W201−τ0+12eiω10τ0W2030+e−iω10τ0W1130,W20θ=ig20q0ω10eiω10θ+ig¯02q¯03ω10e−iω10θ+E1e2iω10θ,W11θ=−ig11q0ω10eiω10θ+ig¯11q¯0ω10e−iω10θ+E2,where E1=E11,E12,E13T∈R3 and E2=E21,E22,E23T∈R3 can be determined by the following, respectively:(50)E1=22iω10−a11−a12−a13−a21−b21e−2iω10τ02iω10−a220−a31−b31e−2iω10τ002iω10−1k11k21k31,E2=−a11−a12−a13−a21−b21−a220−a31−b3100−1k12k22k32.So we can getg21 and compute the following values:(51)c10=i2ω10τ0g20g11−2g112−g0223+g212,ν1=−Rec10Reλ′τ0,ϱ1=2Rec10,T1=−Imc10+ν1Imλ′τ0ω10τ0,which determine the direction of Hopf bifurcation and the stability of bifurcating periodic solutions when τ=τ0. From the conclusions of Hassard et al. [51], we have the following results.Theorem 7.
SupposeH5, H6, H8, and f′yii=1,2,and3 are satisfied. For the properties of Hopf bifurcation of model (21), we have that(1)
Ifν1>0<0, then the direction of the Hopf bifurcation is forward (backward) and the bifurcating periodic solutions exist for τ>τ0τ<τ0(2)
Ifϱ1<0>0, then the bifurcating periodic solutions are orbitally stable (unstable)(3)
IfT1>0<0, then the period of bifurcation periodic solutions increase (decrease)


## 4. Hopf Bifurcation of System (3)

### 4.1. The Existence of Hopf Bifurcation of System (3)

In this section, we study the dynamic behavior of system (3) at the coexistence equilibrium. According to the properties of the Laplacian operator defined on the bounded domain with homogeneous Neumann boundary conditions [53], we get that the operator Δ has the eigenvalues −n2/l2n∈N0≜0,1,2,… with the corresponding eigenfunctions on X denoted by ηn1=cosnx,0,0T,ηn2=0,cosnx,0T,andηn3=0,0,cosnxTn∈N0. Here, we assume that l=1.Under the conditionH4, the linearization of system (3) at E∗ in the phase space C=C−τ,0,X can be rewritten as follows:(52)dUtdt=DΔUt+L1Ut+L2Ut−τ,where D=diagd1,d2,d3, L1 and L2 are defined in Section 3.1, and Ut=u1t,u2t,u3tT=u1t,⋅,u2t,⋅,u3t,⋅T.The characteristic equation of model (52) is(53)Pλ,τ=λ3+h2nλ2+h1nλ+h0n+e−λτc1nλ+c0n=0,where(54)h2n=d1n2+d2n2+d3n2−a11−a22,c1n=−a12b21+a13b31,h1n=d1d2+d1d3+d2d3n4−a11+a22d2+d3n2+a11a22−a12a21−a13a31,h0n=d1d2d3n6−d3d1a22+d2a11n4+d3n2a11a22−a12a21−a13a31d2n2+a13a31a22,c0n=a13a22b31−d3n2a12b21−d2n2a13b31.Next, we seek for condition determining the stability of the coexistence equilibriumE∗ of system (3) with τ=0. For this purpose, substituting τ=0 into equation (53), we obtain that(55)Pλ,0=λ3+h2nλ2+h1n+c1nλ+h0n+c0n=0.Simple computation yields that(56)h1n+c1n=d1d2+d1d3+d2d3n4−a11d2+d3n2−a22d1+d3n2+a11a22−a12a21−a13a31−a12b21−a13b31,h0n+c0n=d1d2d3n6−d3d1a22+d2a11n4+d3n2a11a22−a12a21−a12b21−d2n2a13b31+a13a31+a13a31a22+a13a22b31,h2nh1n+c1n−h0n−c0n=d1nn2h1n+c1n−d2d3n4+d2nn2h1n+c1n+a13b31+a13a31+d3n2h1n+c1n−a11a22+a12a21+a12b21−a11h1n+c1n−a22h1n+c1n+a13a31+a13b31.Under the conditionH6, h1n+c1n>0,h0n+c0n>0, and h2nh1n+c1n−h0n−c0n>0, for all n∈N0. Through the Routh–Hurwitz criterion [48], all the roots of equation (55) have negative real parts. From Theorems 5.1.1 and 5.1.3 of Henry [54], system (3) at E∗ is locally asymptotically stable when τ=0. To sum up, diffusion does not change the stability of system (3) at equilibrium E∗. Then, we can get the following theorem.Theorem 8.
Assume the conditionsH5 and H6 hold, then the coexistence equilibrium E∗P¯,Z¯1,Z¯2 of system (3) with τ=0 is locally asymptotically stable.Now, we study the effect of the delayτ≠0 on the stability of the positive equilibrium E∗ of system (3). First, we seek critical values of τ such that equation (53) has a pair of purely imaginary roots.Letλ=iω2ω2>0 be a root of equation (53); then, substituting λ=iω2 into equation (53) yields(57)−iω23−h2nω22+ih1nω2+h0n+e−iω2τic1nω2+c0n=0,which we get(58)−ω23+h1nω2=c0nsinω2τ−c1nω2cosω2τ,−ω22h2n+h0n=−c1nω2sinω2τ−c0ncosω2τ.Two equations of equation (58) are squared on both sides and then added, which yields(59)ω26+h2n2−2h1nω24+h1n2−2h0nh2n−c1n2ω22+h0n2−c0n2=0.Letω22=y1; then, equation (59) becomes(60)gy1=y13+pny12+qny1+en=0,where(61)pn=h2n2−2h1n=d12+d22+d32n4−2d1n2a11−2d2n2a22+a112+a222+2a12a21+2a13a31,qn=h1n2−2h0nh2n−c1n2=ℱn+a11a22−a12a21−a13a312−2a11a22−a12a21−a13a31a13a31a22−a12b21+a13b312,en=h0n2−c0n2,h0n−c0n=d1d2d3n6−d3d1a22+d2a11n4+d3a11a22−a12a21+a12b21n2+d2a13b31−a31n2+a22a13a31−b31,ℱn=d1d2d3d1d2+d1d3+d2d3n10+d1d2+d1d3+d2d32+d3d1d2+d1d3+d2d3d1a22+d2a11+2a11d2+d3d1d2d3+2a22d1+d3d1d2d3n8+2−a11d2+d3d1d2+d1d3+d2d3−a22d1+d3d1d2+d1d3+d2d3−d1d2+d1d3+d2d3d3a11a22−a12a21+a13a31d2d1d2+d1d3+d2d3−a11d3d2+d3d1a22+d2a11−a22d1+d3d1a22+d2a11−a11a22−a12a21−a13a31d1d2d3n6+a112d2+d32+a222d2+d32+2a11a22−a12a21−a13a31d1d2+d1d3+d2d3+2a11a22d2+d3d2+d3−2a13a31a22d1d2+d1d3+d2d3+2a11d3d2+d3d1a22+d2a11−2a11a13a31d2d2+d3+2a22d3d1+d3a11a22+a12a21−2a22a13a31d2d1+d3+2a11a22−a12a21−a13a31d3d1a22+d2a11n4+−2a11d3d2+d3a11a22−a12a21+a13a31−2a22d1+d3a11a22−a12a21+a13a31+2a11a22a13a31d2+d3+2a222a13a31d1+d2−2a11a22−a12a21+a13a31d3a11a22−a12a21+2a11a22−a12a21+a13a31a13a31d2n2.Under the conditionH6, for all n∈N0,g0=en=h0n2−c0n2>0,limy1⟶∞gy1=∞, and g′y1=3y12+2pny1+qn.From Lemmas 2.2 and 4.2 in [50], we obtain the following lemma.Lemma 3.
For polynomial equation (60), we can get the following conclusions about the distribution of its positive roots.(1)
If there existsn1∈N0, such that pn1>0, qn1>0, and en1>0, then equation (60) has no positive roots(2)
If there existsn2∈N0, such that en2<0, then equation (60) has at least one positive root(3)
Ifen≥0 and Δ5=pn2−3qn>0, then equation (60) has positive roots if and only if y1+=−pn+Δ5/3>0 and gy1+≤0Next, we assume thatH9: Δ6=p02−3q0>0, y1+=−p0+Δ6/3>0, and gy1+≤0. Thus, we have the following result.Theorem 9.
If the conditionsH6 and H9 hold, then equation (60) has at least one positive root for n=0.Proof. Under the conditionH6, then h00+c00>0 and h00−c00>0; hence, e0>0. If H9 holds, according to Lemma 3, then the proof is completed.□Assume that equation (60) has at least one positive root for n=0. Without loss of generality, we suppose that it has three positive roots, denoted by y11,0,y12,0, and y13,0. Hence, equation (59) has three positive roots ω2k,0=y1k,0, k=1,2,and3. Through equation (58), we can get(62)τk,j0=1ω2k,0arccosc10ω2k,04+h20c00−h10c10ω2k,02−h00c00c102ω2k,02+c002+2jπ,for k=1,2,and3,j∈N0. So ±iω2k,0 is a pair of purely imaginary roots of equation (55) with τ=τk,j0. Define τ00=mink=1,2,3,j∈N0τk,j0, and ω20=ω2k,0|τ=τ00.Letλτ=στ±ετi be the roots of equation (53) near τ=τk,j0; then, στk,j0=0 and ετk,j0=ω2k,0j∈N0. So we get the following transversality condition.Lemma 4.
Assume thatg′y1k,0≠0; then, the following transversality condition holds: signdReλτ/dτλ=iω20≠0.Proof.
Differentiating the two sides of equation (53) with respect to τ yields(63)dλdτ−1=3λ2+2h20λ+h10λc10λ+c00+c10λc10λ+c00−τλ.
Thus, we have(64)dReλτdτλ=iω20−1=Re3λ2+2h20λ+h10λc10λ+c00λ=iω20+Rec10λc10λ+c00λ=iω20.
Via calculations, we get(65)dReλτdτλ=iω20−1=3ω206+2h202−2h10ω204+h102−2h20h00−c102ω202Λ2,=y1k,0Λ23y1k,02+2h202−2h10y1k,0+h102−2h20h00−c102=y1k,0Λ2g′y1k,0≠0,where Λ2=c10ω204+c00ω202. Then, we have(66)signdReλτdτλ=iω20=signdReλτdτλ=iω20−1≠0.
This completes the proof of the lemma.By the theory of Hopf bifurcation in Wu [55], we have the results on the existence of spatially homogeneous Hopf bifurcation.Theorem 10.
Suppose that the conditionsH5 and H6 hold, then the following statements about system (3) are true.(i)
If the conditionH8 holds, then the coexistence equilibrium E∗ of system (3) is locally asymptotically stable for τ∈0,τ00; E∗ is unstable when τ>τ00(ii)
Ifg′y1k,0≠0,k=1,2,and3, holds, then system (3) undergoes a spatially homogeneous Hopf bifurcation at the coexistence equilibrium E∗ when τ=τ00To study the existence of spatially inhomogeneous Hopf bifurcation, we need to seek the conditions which satisfy that equation (60) has at least one positive root. From Lemma 3, we know that there exists N1∈N0N1≥1, such that pN1>0, qN1>0, and eN1>0. When, n>N1; equation (60) has no positive roots. When n<N1, Δ5=pn2−3qn>0, y1+=−pn+Δ5/3>0, and gy1+≤0 are satisfied. We know that equation (60) has at least one positive root for 1≤n≤N1. Without loss of generality, we assume that equation (60) has three positive roots, which are denoted by y1k,n0, k=1,2,and3. Then, equation (59) has three positive roots ω2k,n0=y1k,n0.From equation (58), we can get(67)τk,jn0=1ω2k,n0c1n0ω2k,n04+h2n0c0n0−c1n0h1n0ω2k,n02+c0n0h0n0c0n0+c1n02ω2k,n02+2jπω2k,n0,k=1,2,and3andj∈N0.Then,±iω2k,n0 is a pair of purely imaginary roots of equation (53) with τ=τk,jn0. Since g′y1k,n0≠0, using the same calculation method in Lemma 4, we can get signdReλτ/dττ=τk,jn0≠0.By the theory of Hopf bifurcation in Hassard et al. [51], we have the results on the existence of spatially inhomogeneous Hopf bifurcation.Theorem 11.
Suppose that the conditionsH5 and H6 are true, and there exists n0∈N1 such that pn0, qn0, and en0. Moreover, if g′yk,n0≠0k=1,2,and3 holds, then system (3) undergoes a spatially inhomogeneous Hopf bifurcation at E∗ when τ=τk,jn0, k=1,2,and3andj∈N0.To sum up, system (3) without diffusion undergoes a spatially homogeneous Hopf bifurcation at the coexistence equilibrium E∗. System (3) undergoes a spatially inhomogeneous Hopf bifurcation at the coexistence equilibrium E∗.


### 4.2. Direction and Stability of Bifurcated Periodic Solutions

We have get some conditions, such that a spatially homogeneous or inhomogeneous periodic solutions bifurcate from the coexistence equilibriumE∗ of system (3) when τ crosses through the critical values τk,jnk=1,2,and3andj∈N0,n=0andn0. In this part, we will study the direction of Hopf bifurcation and the stability of bifurcated periodic solutions by utilizing the normal form theory and center manifold theorem for partial functional differential equations [55]. For convenience, we fix k=1,2,and3andj=0,1,2,…, and denote τk,jnn=0andn0 by τ˜ and ω2k,nn=0andn0 by ωn.First, if we introduce a new parameterτ=τ˜−υ,υ∈R, then υ=0 is the new Hopf bifurcation value. Let χ=t/τ and rewrite χ as t; system (3) can be rewritten in the phase space C1=C1−1,0,X as(68)U˙t=τ˜DΔUt+Γτ˜Ut+HUt,υ,where u1t=u1t,⋅,u2t=u2t,⋅,u3t=u3t,⋅,U=u1t,u2t,u3tT,D=diagd1,d2,d3, Γθ⋅:C1⟶X, and H:C1×R⟶X are given by(69)Γθφ=θa11φ10+a12φ20+a13φ30a21φ10+a22φ20+b21φ1−1a31φ10+b31φ1−1,Hφ,υ=υDΔφ0+Γυφ+hφ,υ,where(70)hφ,υ=τ˜+υa14φ120+a15φ220+a16φ10φ20+a17φ10φ30a24φ120+a25φ220+a26φ10φ20+b24φ12−1+b25φ1−1φ20a34φ120+a35φ10φ30+b34φ12−1+b35φ1−1φ30,for φ=φ1,φ2,φ3T∈C1.Linearizing model (68) at the origin leads to the following equation:(71)U˙t=τ˜DΔUt+Γτ˜Ut.From the discussion in Section4.1, we can get that the characteristic equation of equation (71) has a pair of purely imaginary roots Υn=iωnτ˜,−iωnτ˜ when τ=τ˜.Let℘:=C1−1,0,R3. Next, we consider the functional differential equation on ℘:(72)z˙=Γτ˜zt.Obviously,Γτ˜ is continuous linear function mapping C1−1,0,R3 into R3. According to the Riesz representation theorem [52], there exists a 3×3 matrix function η1θ,τ˜,n−1<θ<0, whose entries are bounded variation that(73)Γτ˜ϕ=∫−10dη1θ,τ˜,nϕθ,ϕ∈C1−1,0,R3.We can choose(74)η1θ,τ˜,n=τ˜a11a12a13a21a220a3100δθ−τ˜000b2100b3100δθ+1,where δ is a Dirac delta function.DefineC1∗=C−1,0,R3∗, where R3∗ is the three-dimensional vector space of row vectors. Aτ˜ denotes the infinitesimal generator of semigroup induced by the solutions of equation (72), and Aτ˜∗ is the adjoint matrix of Aτ˜ under the bilinear pairing:(75)φ,ϕ=φ0,ϕ0−∫−10∫ξ=0θφξ−θdη1θ,τ˜,nϕξdξ,where ϕ∈C1, φ∈C1−1,0,R3∗. Thus, Aτ˜ and Aτ˜∗ are a pair of adjoint operators, and their eigenvalues are Υn=±iωnτ˜. Let F and F∗ be the center space, that is, the generalized characteristic space of Aτ˜ and Aτ˜∗, respectively. Then, F∗ is adjoint space of F, and dim F = dim F∗. Thus, we have the following lemma.Lemma 5.
Let(76)γ1=a21+b21e−iωnτ˜iωn−a22,γ2=iωn−a11iωn−a22−a12a21+b21e−iωnτ˜a13iωn−a22,γ1∗=−a12iωn+a22,γ2∗=a12a21+b21eiωnτ˜−iωn+a11iωn+a22iωn+a22a31+b31eiωnτ˜,and then, p1θ=1,γ1,γ2Teiωnτ˜θ and p2θ=p1θ¯,−1<θ<0 are a basis of F associated with Υn, and p1∗ξ=1,γ1∗,γ2∗e−iωnτ˜θ and p2∗ξ=p1∗ξ¯, 0≤ξ≤1, are a basis of F∗ associated with Υn.DenotingΦ=Φ1,Φ2andΨ∗=Ψ1∗,Ψ2∗T, where(77)Φ1=p1θ+p2θ2,Φ2=p1θ−p2θ2i,θ∈−1,0,Ψ1∗=p1∗ξ+p2∗ξ2,Ψ1∗=p1∗ξ−p2∗ξ2i,ξ∈0,1,and we define that(78)Ψ∗,Φ=Ψk∗,Φl=Ψ1∗,Φ1Ψ1∗,Φ2Ψ2∗,Φ1Ψ2∗,Φ2,and construct a new basis Ψ for F∗ by Ψ=Ψ1,Ψ2T=Ψ∗,Φ−1Ψ∗, so Ψ,Φ=I.Letfn=ηn1,ηn2,ηn3, where ηn1=cosnx,0,0T,ηn2=0,cosnx,0T,andηn3=0,0,cosnxT. Define(79)v⋅fn=v1ηn1+v2ηn2+v3ηn3,v=v1,v2,v3T∈C1−1,0,X.From Wu [55], the center space of linear equation (71) is given by PCN℘, where(80)PCNφ=ΦΨ,φ,fn⋅fn,φ∈℘,and C1=PCN℘⊕PSC1; PSC1 denotes the complementary subspace of PCN℘.According to Hassard et al. [51] and Wu [55], the center space of the linear model of model (68) with μ2=0 is given by PCN℘, where(81)PCN℘=12p1θz+p2θz¯⋅fn,z∈C1.Hence, the solutions of model (68) are(82)Ut=12p1θz+p2θz¯⋅fn+Qzt,z¯tθ,where Qzt,z¯tθ=wz+z¯/2,iz−z¯/2,0, z=x1−ix2. Therefore, z satisfies(83)z˙=iωnτ˜z+g1z,z¯,where g1z,z¯=Ψ10−iΨ20hUt,0,fn, Ψ0=Ψ10,Ψ20T.Let(84)Qz,z¯=Q20z2z+Q11zz¯+Q02z¯2z+⋯,(85)g1z,z¯=g˜20z22+g˜11zz¯+g˜02z¯2z+g˜21z2z¯2+⋯,and κ1,κ2,κ3=Ψ10−iΨ20. From equations (70) and (85), we can get(86)g˜20=τ˜2a14+γ12a15+γ1a16+γ2a17κ1+a24+γ12a25+γ1a26+b24e−i2ωnτ˜+γ1b25e−iωnτ˜κ2+a34+γ2a35+b34e−i2ωnτ˜+γ2b35e−iωnτ˜κ3cos2nx,cosnx,g˜11=τ˜42a14+2γ1γ¯1a25+γ1+γ¯1a16+γ2+γ¯2a17κ1+2a24+2γ1γ¯1a25+γ1+γ¯1a26+2b24+γ1eiωnτ˜+γ¯1e−iωnτ˜b25κ2+2a34+γ2+γ¯2a35+2b34+γ2eiωnτ˜+γ¯2e−iωnτ˜b35κ3cos2nx,cosnx,g˜21=τ˜a14Q2010+2Q1110+a15Q2020+2Q1120+a16γ1Q1110+12γ¯1Q2010+Q1120+12Q2020+a17γ2Q1110+12γ¯2Q2010+Q1130+12Q2030cosnx,cosnxκ1+a24Q2010+2Q1110+a25Q2020+2Q1120+a26γ1Q1110+12γ¯1Q2010+Q1120+12Q2020+b242e−iωnτ˜Q111−1+eiωnτ˜Q201−1+b25γ1Q111−1+12γ¯1Q201−1+e−iωnτ˜Q1120+12eiωnτ˜Q2020cosnx,cosnxκ2+a34Q2010+2Q1110+a35γ2Q11102Q1110+12γ¯2Q2010+Q1130+12Q2030+b342e−iωnτ˜Q111−1+eiωnτ˜Q201−1+b35γ2Q111−1+12γ¯2Q201−1+e−iωnτ˜Q1130+12eiωnτ˜Q2030cosnx,cosnxκ3,where n∈N0. It is obvious that ∫0πcos3nxdx=0 and g˜02=g˜20¯. When n=1,2,…, we have g˜20=g˜11=g˜02=0. When n=0, we obtain(87)g˜20=τ˜2a14+γ12a15+γ1a16+γ2a17κ1+a24+γ12a25+γ1a26+b24e−i2ωnτ˜+γ1b25e−iωnτ˜κ2+a34+γ2a35+b34e−i2ωnτ˜+γ2b35e−iωnτ˜κ3,g˜11=τ˜42a14+2γ1γ¯1a15+γ1+γ¯1a16+γ2+γ¯2a17κ1+2a24+2γ1γ¯1a25+γ1+γ¯1a26+2b24+γ1eiωnτ˜+γ¯1e−iωnτ˜b25κ2+2a34+γ2+γ¯2a35+2b34+γ2eiωnτ˜+γ¯2e−iωnτ˜b35κ3,and Q200,Q20−1,Q110,andQ11−1 appear in g21. From equation (84), we have(88)Q˙z,z¯=Q20zz˙+Q11z˙z¯+Q11zz¯˙+Q02z¯z˙¯+⋯,(89)Aτ˜Q=Aτ˜Q20z22+Aτ˜Q11zz¯+Aτ˜Q02z¯22+⋯.Therefore, according to Wu [55], we know that Qzt,z¯t satisfies(90)Q˙=Aτ˜Q+Gz,z¯,where(91)Gz,z¯=G20z22+G11zz¯+G02z¯22+⋯=X0hUt,0−ΦΨ,X0hUt,0,fn⋅fn,where the function X0 defined by X00=IandX0θ=0,−1≤θ<0, and Gij∈PS℘i+j=2.According to (82) and (88)–(90), we have(92)2iωnτ˜−Aτ˜Q20=G20,−Aτ˜Q11=G11,−2iωnτ˜−Aτ˜Q02=G02.SinceAτ˜ has a pair of purely imaginary characteristic roots, which are ±iωnτ˜, equation (91) has the unique solution Qij.From equation (92), we obtain that, for θ∈−1,0,(93)Gz,z˜=−ΦθΨ0hUt,0,fn⋅fn,=−12τ˜p1θg˜20+g˜02¯p1θ¯⋅fnz22+p1θg˜11+g˜11¯p1θ¯⋅fnzz¯+p1θg˜02+g˜20¯p1θ¯⋅fnz¯22+⋯.Therefore, by comparing the coefficients of equations (91) and (93) and noticing that(94)Gz,z˜0=hUt,0−ΦΨ,hUt,0,fn⋅fn=hUt,0−ΦθΨ0hUt,0,fn⋅fn,we obtain that(95)G20θ=−12τ˜p1θg˜20+g˜20¯p1θ¯cosnx,θ∈−1,0,τ˜2a14+γ12a15+γ1a16+γ2a17a24+γ12a25+γ1a26+b24e−2iωnτ˜+γ1b25e−iωnτ˜a34+γ2a35+b34e−2iωnτ˜+γ2b35e−iωnτ˜cos2nx−12τ˜p1θg˜20+g˜20¯p1θ¯cosnx,θ=0,(96)G11θ=−12τ˜p1θg˜11+g˜11¯p1θ¯cosnx,θ∈−1,0,τ˜42a14+2γ1γ¯1a25+γ1+γ¯1a16+γ2+γ¯2a172a24+2γ1γ¯1a25+γ1+γ¯1a26+2b24+γ1eiωnτ˜+γ¯1e−iωnτ˜b252a34+γ2+γ¯2a35+2b34+γ2eiωnτ˜+γ¯2e−iωnτ˜b35cos2nx−12τ˜p1θg˜11+g˜11¯p1θ¯cosnx,θ=0.Next, we computeQ200, Q20−1, Q110, and Q11−1 by the variation-of-constants method. According to equations (92)–(96), we have(97)Q20θ=−12iωnτ˜p10g˜20eiωnτ˜θ+p10¯g˜02¯3e−iωnτ˜θ+e2iωnτ˜θE3,Q11θ=12iωnτ˜p10g˜11eiωnτ˜θ−p10¯g˜11¯e−iωnτ˜θ+E4.Furthermore, forQij∈C1, we have(98)Q200=−12iωnτ˜p10g˜20+p10¯g˜02¯3+E3,Q110=12iωnτ˜p10g˜11−p10¯g˜11¯+E4,where E3=E31,E32,E33T∈R3 and E4=E41,E42,E43T∈R3 are given as follows, respectively:(99)E3=122iωn−a11−a12−a13−a21−b21e−2iωnτ˜2iωn−a220−a31−b31e−2iωnτ˜02iωn−1×a14+γ12a15+γ1a16+γ2a17a24+γ12a25+γ1a26+b24e−2iωnτ˜+γ1b25e−iωnτ˜a34+γ2a35+b34e−2iωnτ˜+γ2b35e−iωnτ˜cos2nx,E4=14−a11−a12−a13−a21−b21−a220−a31−a3100−1×2a14+2γ1γ¯1a25+γ1+γ¯1a16+γ2+γ¯2a172a24+2γ1γ¯1a25+γ1+γ¯1a26+2b24+γ1eiωnτ˜+γ¯1e−iωnτ˜b252a34+γ2+γ¯2a35+2b34+γ2eiωnτ˜+γ¯2e−iωnτ˜b35cos2nx.In the following, we can compute the following quantities that determine the properties of bifurcating periodic solutions:(100)c20=i2ωnτ˜g˜20g˜11−2g˜112−g˜0223+g˜212,ν2=−Rec20Reλ′τ˜,ϱ2=2Rec20,T2=−Imc20+ν2Imλ′τ˜ωnτ˜.Then, by the Hopf bifurcation theory [51], we can obtain the following result.Theorem 12.
Assume that Theorem9 or Theorem 10 holds, the following results on equation (100) are true.(1)
ν2 determines the Hopf bifurcation direction. If ν2>0<0, then the Hopf bifurcation is forward (backward) and the bifurcating periodic solutions exist for τ>τ˜.(2)
ϱ2 determines the stability of bifurcating periodic solutions. If ϱ2<0>0, then the periodic solutions are stable (unstable).(3)
T2 determines the period of the bifurcating periodic solutions. If T2>0<0, then the period increases (decreases).


## 4.1. The Existence of Hopf Bifurcation of System (3)

In this section, we study the dynamic behavior of system (3) at the coexistence equilibrium. According to the properties of the Laplacian operator defined on the bounded domain with homogeneous Neumann boundary conditions [53], we get that the operator Δ has the eigenvalues −n2/l2n∈N0≜0,1,2,… with the corresponding eigenfunctions on X denoted by ηn1=cosnx,0,0T,ηn2=0,cosnx,0T,andηn3=0,0,cosnxTn∈N0. Here, we assume that l=1.Under the conditionH4, the linearization of system (3) at E∗ in the phase space C=C−τ,0,X can be rewritten as follows:(52)dUtdt=DΔUt+L1Ut+L2Ut−τ,where D=diagd1,d2,d3, L1 and L2 are defined in Section 3.1, and Ut=u1t,u2t,u3tT=u1t,⋅,u2t,⋅,u3t,⋅T.The characteristic equation of model (52) is(53)Pλ,τ=λ3+h2nλ2+h1nλ+h0n+e−λτc1nλ+c0n=0,where(54)h2n=d1n2+d2n2+d3n2−a11−a22,c1n=−a12b21+a13b31,h1n=d1d2+d1d3+d2d3n4−a11+a22d2+d3n2+a11a22−a12a21−a13a31,h0n=d1d2d3n6−d3d1a22+d2a11n4+d3n2a11a22−a12a21−a13a31d2n2+a13a31a22,c0n=a13a22b31−d3n2a12b21−d2n2a13b31.Next, we seek for condition determining the stability of the coexistence equilibriumE∗ of system (3) with τ=0. For this purpose, substituting τ=0 into equation (53), we obtain that(55)Pλ,0=λ3+h2nλ2+h1n+c1nλ+h0n+c0n=0.Simple computation yields that(56)h1n+c1n=d1d2+d1d3+d2d3n4−a11d2+d3n2−a22d1+d3n2+a11a22−a12a21−a13a31−a12b21−a13b31,h0n+c0n=d1d2d3n6−d3d1a22+d2a11n4+d3n2a11a22−a12a21−a12b21−d2n2a13b31+a13a31+a13a31a22+a13a22b31,h2nh1n+c1n−h0n−c0n=d1nn2h1n+c1n−d2d3n4+d2nn2h1n+c1n+a13b31+a13a31+d3n2h1n+c1n−a11a22+a12a21+a12b21−a11h1n+c1n−a22h1n+c1n+a13a31+a13b31.Under the conditionH6, h1n+c1n>0,h0n+c0n>0, and h2nh1n+c1n−h0n−c0n>0, for all n∈N0. Through the Routh–Hurwitz criterion [48], all the roots of equation (55) have negative real parts. From Theorems 5.1.1 and 5.1.3 of Henry [54], system (3) at E∗ is locally asymptotically stable when τ=0. To sum up, diffusion does not change the stability of system (3) at equilibrium E∗. Then, we can get the following theorem.Theorem 8.
Assume the conditionsH5 and H6 hold, then the coexistence equilibrium E∗P¯,Z¯1,Z¯2 of system (3) with τ=0 is locally asymptotically stable.Now, we study the effect of the delayτ≠0 on the stability of the positive equilibrium E∗ of system (3). First, we seek critical values of τ such that equation (53) has a pair of purely imaginary roots.Letλ=iω2ω2>0 be a root of equation (53); then, substituting λ=iω2 into equation (53) yields(57)−iω23−h2nω22+ih1nω2+h0n+e−iω2τic1nω2+c0n=0,which we get(58)−ω23+h1nω2=c0nsinω2τ−c1nω2cosω2τ,−ω22h2n+h0n=−c1nω2sinω2τ−c0ncosω2τ.Two equations of equation (58) are squared on both sides and then added, which yields(59)ω26+h2n2−2h1nω24+h1n2−2h0nh2n−c1n2ω22+h0n2−c0n2=0.Letω22=y1; then, equation (59) becomes(60)gy1=y13+pny12+qny1+en=0,where(61)pn=h2n2−2h1n=d12+d22+d32n4−2d1n2a11−2d2n2a22+a112+a222+2a12a21+2a13a31,qn=h1n2−2h0nh2n−c1n2=ℱn+a11a22−a12a21−a13a312−2a11a22−a12a21−a13a31a13a31a22−a12b21+a13b312,en=h0n2−c0n2,h0n−c0n=d1d2d3n6−d3d1a22+d2a11n4+d3a11a22−a12a21+a12b21n2+d2a13b31−a31n2+a22a13a31−b31,ℱn=d1d2d3d1d2+d1d3+d2d3n10+d1d2+d1d3+d2d32+d3d1d2+d1d3+d2d3d1a22+d2a11+2a11d2+d3d1d2d3+2a22d1+d3d1d2d3n8+2−a11d2+d3d1d2+d1d3+d2d3−a22d1+d3d1d2+d1d3+d2d3−d1d2+d1d3+d2d3d3a11a22−a12a21+a13a31d2d1d2+d1d3+d2d3−a11d3d2+d3d1a22+d2a11−a22d1+d3d1a22+d2a11−a11a22−a12a21−a13a31d1d2d3n6+a112d2+d32+a222d2+d32+2a11a22−a12a21−a13a31d1d2+d1d3+d2d3+2a11a22d2+d3d2+d3−2a13a31a22d1d2+d1d3+d2d3+2a11d3d2+d3d1a22+d2a11−2a11a13a31d2d2+d3+2a22d3d1+d3a11a22+a12a21−2a22a13a31d2d1+d3+2a11a22−a12a21−a13a31d3d1a22+d2a11n4+−2a11d3d2+d3a11a22−a12a21+a13a31−2a22d1+d3a11a22−a12a21+a13a31+2a11a22a13a31d2+d3+2a222a13a31d1+d2−2a11a22−a12a21+a13a31d3a11a22−a12a21+2a11a22−a12a21+a13a31a13a31d2n2.Under the conditionH6, for all n∈N0,g0=en=h0n2−c0n2>0,limy1⟶∞gy1=∞, and g′y1=3y12+2pny1+qn.From Lemmas 2.2 and 4.2 in [50], we obtain the following lemma.Lemma 3.
For polynomial equation (60), we can get the following conclusions about the distribution of its positive roots.(1)
If there existsn1∈N0, such that pn1>0, qn1>0, and en1>0, then equation (60) has no positive roots(2)
If there existsn2∈N0, such that en2<0, then equation (60) has at least one positive root(3)
Ifen≥0 and Δ5=pn2−3qn>0, then equation (60) has positive roots if and only if y1+=−pn+Δ5/3>0 and gy1+≤0Next, we assume thatH9: Δ6=p02−3q0>0, y1+=−p0+Δ6/3>0, and gy1+≤0. Thus, we have the following result.Theorem 9.
If the conditionsH6 and H9 hold, then equation (60) has at least one positive root for n=0.Proof. Under the conditionH6, then h00+c00>0 and h00−c00>0; hence, e0>0. If H9 holds, according to Lemma 3, then the proof is completed.□Assume that equation (60) has at least one positive root for n=0. Without loss of generality, we suppose that it has three positive roots, denoted by y11,0,y12,0, and y13,0. Hence, equation (59) has three positive roots ω2k,0=y1k,0, k=1,2,and3. Through equation (58), we can get(62)τk,j0=1ω2k,0arccosc10ω2k,04+h20c00−h10c10ω2k,02−h00c00c102ω2k,02+c002+2jπ,for k=1,2,and3,j∈N0. So ±iω2k,0 is a pair of purely imaginary roots of equation (55) with τ=τk,j0. Define τ00=mink=1,2,3,j∈N0τk,j0, and ω20=ω2k,0|τ=τ00.Letλτ=στ±ετi be the roots of equation (53) near τ=τk,j0; then, στk,j0=0 and ετk,j0=ω2k,0j∈N0. So we get the following transversality condition.Lemma 4.
Assume thatg′y1k,0≠0; then, the following transversality condition holds: signdReλτ/dτλ=iω20≠0.Proof.
Differentiating the two sides of equation (53) with respect to τ yields(63)dλdτ−1=3λ2+2h20λ+h10λc10λ+c00+c10λc10λ+c00−τλ.
Thus, we have(64)dReλτdτλ=iω20−1=Re3λ2+2h20λ+h10λc10λ+c00λ=iω20+Rec10λc10λ+c00λ=iω20.
Via calculations, we get(65)dReλτdτλ=iω20−1=3ω206+2h202−2h10ω204+h102−2h20h00−c102ω202Λ2,=y1k,0Λ23y1k,02+2h202−2h10y1k,0+h102−2h20h00−c102=y1k,0Λ2g′y1k,0≠0,where Λ2=c10ω204+c00ω202. Then, we have(66)signdReλτdτλ=iω20=signdReλτdτλ=iω20−1≠0.
This completes the proof of the lemma.By the theory of Hopf bifurcation in Wu [55], we have the results on the existence of spatially homogeneous Hopf bifurcation.Theorem 10.
Suppose that the conditionsH5 and H6 hold, then the following statements about system (3) are true.(i)
If the conditionH8 holds, then the coexistence equilibrium E∗ of system (3) is locally asymptotically stable for τ∈0,τ00; E∗ is unstable when τ>τ00(ii)
Ifg′y1k,0≠0,k=1,2,and3, holds, then system (3) undergoes a spatially homogeneous Hopf bifurcation at the coexistence equilibrium E∗ when τ=τ00To study the existence of spatially inhomogeneous Hopf bifurcation, we need to seek the conditions which satisfy that equation (60) has at least one positive root. From Lemma 3, we know that there exists N1∈N0N1≥1, such that pN1>0, qN1>0, and eN1>0. When, n>N1; equation (60) has no positive roots. When n<N1, Δ5=pn2−3qn>0, y1+=−pn+Δ5/3>0, and gy1+≤0 are satisfied. We know that equation (60) has at least one positive root for 1≤n≤N1. Without loss of generality, we assume that equation (60) has three positive roots, which are denoted by y1k,n0, k=1,2,and3. Then, equation (59) has three positive roots ω2k,n0=y1k,n0.From equation (58), we can get(67)τk,jn0=1ω2k,n0c1n0ω2k,n04+h2n0c0n0−c1n0h1n0ω2k,n02+c0n0h0n0c0n0+c1n02ω2k,n02+2jπω2k,n0,k=1,2,and3andj∈N0.Then,±iω2k,n0 is a pair of purely imaginary roots of equation (53) with τ=τk,jn0. Since g′y1k,n0≠0, using the same calculation method in Lemma 4, we can get signdReλτ/dττ=τk,jn0≠0.By the theory of Hopf bifurcation in Hassard et al. [51], we have the results on the existence of spatially inhomogeneous Hopf bifurcation.Theorem 11.
Suppose that the conditionsH5 and H6 are true, and there exists n0∈N1 such that pn0, qn0, and en0. Moreover, if g′yk,n0≠0k=1,2,and3 holds, then system (3) undergoes a spatially inhomogeneous Hopf bifurcation at E∗ when τ=τk,jn0, k=1,2,and3andj∈N0.To sum up, system (3) without diffusion undergoes a spatially homogeneous Hopf bifurcation at the coexistence equilibrium E∗. System (3) undergoes a spatially inhomogeneous Hopf bifurcation at the coexistence equilibrium E∗.


## 4.2. Direction and Stability of Bifurcated Periodic Solutions

We have get some conditions, such that a spatially homogeneous or inhomogeneous periodic solutions bifurcate from the coexistence equilibriumE∗ of system (3) when τ crosses through the critical values τk,jnk=1,2,and3andj∈N0,n=0andn0. In this part, we will study the direction of Hopf bifurcation and the stability of bifurcated periodic solutions by utilizing the normal form theory and center manifold theorem for partial functional differential equations [55]. For convenience, we fix k=1,2,and3andj=0,1,2,…, and denote τk,jnn=0andn0 by τ˜ and ω2k,nn=0andn0 by ωn.First, if we introduce a new parameterτ=τ˜−υ,υ∈R, then υ=0 is the new Hopf bifurcation value. Let χ=t/τ and rewrite χ as t; system (3) can be rewritten in the phase space C1=C1−1,0,X as(68)U˙t=τ˜DΔUt+Γτ˜Ut+HUt,υ,where u1t=u1t,⋅,u2t=u2t,⋅,u3t=u3t,⋅,U=u1t,u2t,u3tT,D=diagd1,d2,d3, Γθ⋅:C1⟶X, and H:C1×R⟶X are given by(69)Γθφ=θa11φ10+a12φ20+a13φ30a21φ10+a22φ20+b21φ1−1a31φ10+b31φ1−1,Hφ,υ=υDΔφ0+Γυφ+hφ,υ,where(70)hφ,υ=τ˜+υa14φ120+a15φ220+a16φ10φ20+a17φ10φ30a24φ120+a25φ220+a26φ10φ20+b24φ12−1+b25φ1−1φ20a34φ120+a35φ10φ30+b34φ12−1+b35φ1−1φ30,for φ=φ1,φ2,φ3T∈C1.Linearizing model (68) at the origin leads to the following equation:(71)U˙t=τ˜DΔUt+Γτ˜Ut.From the discussion in Section4.1, we can get that the characteristic equation of equation (71) has a pair of purely imaginary roots Υn=iωnτ˜,−iωnτ˜ when τ=τ˜.Let℘:=C1−1,0,R3. Next, we consider the functional differential equation on ℘:(72)z˙=Γτ˜zt.Obviously,Γτ˜ is continuous linear function mapping C1−1,0,R3 into R3. According to the Riesz representation theorem [52], there exists a 3×3 matrix function η1θ,τ˜,n−1<θ<0, whose entries are bounded variation that(73)Γτ˜ϕ=∫−10dη1θ,τ˜,nϕθ,ϕ∈C1−1,0,R3.We can choose(74)η1θ,τ˜,n=τ˜a11a12a13a21a220a3100δθ−τ˜000b2100b3100δθ+1,where δ is a Dirac delta function.DefineC1∗=C−1,0,R3∗, where R3∗ is the three-dimensional vector space of row vectors. Aτ˜ denotes the infinitesimal generator of semigroup induced by the solutions of equation (72), and Aτ˜∗ is the adjoint matrix of Aτ˜ under the bilinear pairing:(75)φ,ϕ=φ0,ϕ0−∫−10∫ξ=0θφξ−θdη1θ,τ˜,nϕξdξ,where ϕ∈C1, φ∈C1−1,0,R3∗. Thus, Aτ˜ and Aτ˜∗ are a pair of adjoint operators, and their eigenvalues are Υn=±iωnτ˜. Let F and F∗ be the center space, that is, the generalized characteristic space of Aτ˜ and Aτ˜∗, respectively. Then, F∗ is adjoint space of F, and dim F = dim F∗. Thus, we have the following lemma.Lemma 5.
Let(76)γ1=a21+b21e−iωnτ˜iωn−a22,γ2=iωn−a11iωn−a22−a12a21+b21e−iωnτ˜a13iωn−a22,γ1∗=−a12iωn+a22,γ2∗=a12a21+b21eiωnτ˜−iωn+a11iωn+a22iωn+a22a31+b31eiωnτ˜,and then, p1θ=1,γ1,γ2Teiωnτ˜θ and p2θ=p1θ¯,−1<θ<0 are a basis of F associated with Υn, and p1∗ξ=1,γ1∗,γ2∗e−iωnτ˜θ and p2∗ξ=p1∗ξ¯, 0≤ξ≤1, are a basis of F∗ associated with Υn.DenotingΦ=Φ1,Φ2andΨ∗=Ψ1∗,Ψ2∗T, where(77)Φ1=p1θ+p2θ2,Φ2=p1θ−p2θ2i,θ∈−1,0,Ψ1∗=p1∗ξ+p2∗ξ2,Ψ1∗=p1∗ξ−p2∗ξ2i,ξ∈0,1,and we define that(78)Ψ∗,Φ=Ψk∗,Φl=Ψ1∗,Φ1Ψ1∗,Φ2Ψ2∗,Φ1Ψ2∗,Φ2,and construct a new basis Ψ for F∗ by Ψ=Ψ1,Ψ2T=Ψ∗,Φ−1Ψ∗, so Ψ,Φ=I.Letfn=ηn1,ηn2,ηn3, where ηn1=cosnx,0,0T,ηn2=0,cosnx,0T,andηn3=0,0,cosnxT. Define(79)v⋅fn=v1ηn1+v2ηn2+v3ηn3,v=v1,v2,v3T∈C1−1,0,X.From Wu [55], the center space of linear equation (71) is given by PCN℘, where(80)PCNφ=ΦΨ,φ,fn⋅fn,φ∈℘,and C1=PCN℘⊕PSC1; PSC1 denotes the complementary subspace of PCN℘.According to Hassard et al. [51] and Wu [55], the center space of the linear model of model (68) with μ2=0 is given by PCN℘, where(81)PCN℘=12p1θz+p2θz¯⋅fn,z∈C1.Hence, the solutions of model (68) are(82)Ut=12p1θz+p2θz¯⋅fn+Qzt,z¯tθ,where Qzt,z¯tθ=wz+z¯/2,iz−z¯/2,0, z=x1−ix2. Therefore, z satisfies(83)z˙=iωnτ˜z+g1z,z¯,where g1z,z¯=Ψ10−iΨ20hUt,0,fn, Ψ0=Ψ10,Ψ20T.Let(84)Qz,z¯=Q20z2z+Q11zz¯+Q02z¯2z+⋯,(85)g1z,z¯=g˜20z22+g˜11zz¯+g˜02z¯2z+g˜21z2z¯2+⋯,and κ1,κ2,κ3=Ψ10−iΨ20. From equations (70) and (85), we can get(86)g˜20=τ˜2a14+γ12a15+γ1a16+γ2a17κ1+a24+γ12a25+γ1a26+b24e−i2ωnτ˜+γ1b25e−iωnτ˜κ2+a34+γ2a35+b34e−i2ωnτ˜+γ2b35e−iωnτ˜κ3cos2nx,cosnx,g˜11=τ˜42a14+2γ1γ¯1a25+γ1+γ¯1a16+γ2+γ¯2a17κ1+2a24+2γ1γ¯1a25+γ1+γ¯1a26+2b24+γ1eiωnτ˜+γ¯1e−iωnτ˜b25κ2+2a34+γ2+γ¯2a35+2b34+γ2eiωnτ˜+γ¯2e−iωnτ˜b35κ3cos2nx,cosnx,g˜21=τ˜a14Q2010+2Q1110+a15Q2020+2Q1120+a16γ1Q1110+12γ¯1Q2010+Q1120+12Q2020+a17γ2Q1110+12γ¯2Q2010+Q1130+12Q2030cosnx,cosnxκ1+a24Q2010+2Q1110+a25Q2020+2Q1120+a26γ1Q1110+12γ¯1Q2010+Q1120+12Q2020+b242e−iωnτ˜Q111−1+eiωnτ˜Q201−1+b25γ1Q111−1+12γ¯1Q201−1+e−iωnτ˜Q1120+12eiωnτ˜Q2020cosnx,cosnxκ2+a34Q2010+2Q1110+a35γ2Q11102Q1110+12γ¯2Q2010+Q1130+12Q2030+b342e−iωnτ˜Q111−1+eiωnτ˜Q201−1+b35γ2Q111−1+12γ¯2Q201−1+e−iωnτ˜Q1130+12eiωnτ˜Q2030cosnx,cosnxκ3,where n∈N0. It is obvious that ∫0πcos3nxdx=0 and g˜02=g˜20¯. When n=1,2,…, we have g˜20=g˜11=g˜02=0. When n=0, we obtain(87)g˜20=τ˜2a14+γ12a15+γ1a16+γ2a17κ1+a24+γ12a25+γ1a26+b24e−i2ωnτ˜+γ1b25e−iωnτ˜κ2+a34+γ2a35+b34e−i2ωnτ˜+γ2b35e−iωnτ˜κ3,g˜11=τ˜42a14+2γ1γ¯1a15+γ1+γ¯1a16+γ2+γ¯2a17κ1+2a24+2γ1γ¯1a25+γ1+γ¯1a26+2b24+γ1eiωnτ˜+γ¯1e−iωnτ˜b25κ2+2a34+γ2+γ¯2a35+2b34+γ2eiωnτ˜+γ¯2e−iωnτ˜b35κ3,and Q200,Q20−1,Q110,andQ11−1 appear in g21. From equation (84), we have(88)Q˙z,z¯=Q20zz˙+Q11z˙z¯+Q11zz¯˙+Q02z¯z˙¯+⋯,(89)Aτ˜Q=Aτ˜Q20z22+Aτ˜Q11zz¯+Aτ˜Q02z¯22+⋯.Therefore, according to Wu [55], we know that Qzt,z¯t satisfies(90)Q˙=Aτ˜Q+Gz,z¯,where(91)Gz,z¯=G20z22+G11zz¯+G02z¯22+⋯=X0hUt,0−ΦΨ,X0hUt,0,fn⋅fn,where the function X0 defined by X00=IandX0θ=0,−1≤θ<0, and Gij∈PS℘i+j=2.According to (82) and (88)–(90), we have(92)2iωnτ˜−Aτ˜Q20=G20,−Aτ˜Q11=G11,−2iωnτ˜−Aτ˜Q02=G02.SinceAτ˜ has a pair of purely imaginary characteristic roots, which are ±iωnτ˜, equation (91) has the unique solution Qij.From equation (92), we obtain that, for θ∈−1,0,(93)Gz,z˜=−ΦθΨ0hUt,0,fn⋅fn,=−12τ˜p1θg˜20+g˜02¯p1θ¯⋅fnz22+p1θg˜11+g˜11¯p1θ¯⋅fnzz¯+p1θg˜02+g˜20¯p1θ¯⋅fnz¯22+⋯.Therefore, by comparing the coefficients of equations (91) and (93) and noticing that(94)Gz,z˜0=hUt,0−ΦΨ,hUt,0,fn⋅fn=hUt,0−ΦθΨ0hUt,0,fn⋅fn,we obtain that(95)G20θ=−12τ˜p1θg˜20+g˜20¯p1θ¯cosnx,θ∈−1,0,τ˜2a14+γ12a15+γ1a16+γ2a17a24+γ12a25+γ1a26+b24e−2iωnτ˜+γ1b25e−iωnτ˜a34+γ2a35+b34e−2iωnτ˜+γ2b35e−iωnτ˜cos2nx−12τ˜p1θg˜20+g˜20¯p1θ¯cosnx,θ=0,(96)G11θ=−12τ˜p1θg˜11+g˜11¯p1θ¯cosnx,θ∈−1,0,τ˜42a14+2γ1γ¯1a25+γ1+γ¯1a16+γ2+γ¯2a172a24+2γ1γ¯1a25+γ1+γ¯1a26+2b24+γ1eiωnτ˜+γ¯1e−iωnτ˜b252a34+γ2+γ¯2a35+2b34+γ2eiωnτ˜+γ¯2e−iωnτ˜b35cos2nx−12τ˜p1θg˜11+g˜11¯p1θ¯cosnx,θ=0.Next, we computeQ200, Q20−1, Q110, and Q11−1 by the variation-of-constants method. According to equations (92)–(96), we have(97)Q20θ=−12iωnτ˜p10g˜20eiωnτ˜θ+p10¯g˜02¯3e−iωnτ˜θ+e2iωnτ˜θE3,Q11θ=12iωnτ˜p10g˜11eiωnτ˜θ−p10¯g˜11¯e−iωnτ˜θ+E4.Furthermore, forQij∈C1, we have(98)Q200=−12iωnτ˜p10g˜20+p10¯g˜02¯3+E3,Q110=12iωnτ˜p10g˜11−p10¯g˜11¯+E4,where E3=E31,E32,E33T∈R3 and E4=E41,E42,E43T∈R3 are given as follows, respectively:(99)E3=122iωn−a11−a12−a13−a21−b21e−2iωnτ˜2iωn−a220−a31−b31e−2iωnτ˜02iωn−1×a14+γ12a15+γ1a16+γ2a17a24+γ12a25+γ1a26+b24e−2iωnτ˜+γ1b25e−iωnτ˜a34+γ2a35+b34e−2iωnτ˜+γ2b35e−iωnτ˜cos2nx,E4=14−a11−a12−a13−a21−b21−a220−a31−a3100−1×2a14+2γ1γ¯1a25+γ1+γ¯1a16+γ2+γ¯2a172a24+2γ1γ¯1a25+γ1+γ¯1a26+2b24+γ1eiωnτ˜+γ¯1e−iωnτ˜b252a34+γ2+γ¯2a35+2b34+γ2eiωnτ˜+γ¯2e−iωnτ˜b35cos2nx.In the following, we can compute the following quantities that determine the properties of bifurcating periodic solutions:(100)c20=i2ωnτ˜g˜20g˜11−2g˜112−g˜0223+g˜212,ν2=−Rec20Reλ′τ˜,ϱ2=2Rec20,T2=−Imc20+ν2Imλ′τ˜ωnτ˜.Then, by the Hopf bifurcation theory [51], we can obtain the following result.Theorem 12.
Assume that Theorem9 or Theorem 10 holds, the following results on equation (100) are true.(1)
ν2 determines the Hopf bifurcation direction. If ν2>0<0, then the Hopf bifurcation is forward (backward) and the bifurcating periodic solutions exist for τ>τ˜.(2)
ϱ2 determines the stability of bifurcating periodic solutions. If ϱ2<0>0, then the periodic solutions are stable (unstable).(3)
T2 determines the period of the bifurcating periodic solutions. If T2>0<0, then the period increases (decreases).


## 5. Numerical Simulation

In this section, we present some numerical results of system (3) to support and extend our theoretical results which have been established in the previous sections. The spatial domain is 0,π. Next, we simulate those results by using the software Matlab.First, we simulate system (4). We assume that P0=9,Z10=20,andZ20=1 and choose the following parameter values: r=11,K=25,μ1=10,μ2=9,α1=10,α2=20,β1=5,β2=5,ρ1=1,ρ2=1.400,g1=1,g2=1.110,ands=0.700 so that P¯>g1α1/β1−ρ1−g1, rK−P¯β1−ρ1P¯>Kg1μ1Z¯1, and μ1Z¯1/α1+P¯+sZ¯12+μ2P¯/α2+P¯2<P¯/K. According to Theorem 4, we can get that system (4) has the coexistence equilibrium E∗8.916,23.924,1.184 and E∗ is locally asymptotically stable, see Figure 1.Figure 1
The coexistence equilibriumE∗ of system (4) is locally asymptotically stable when τ=0. (a) Pt; (b) Z1t; (c) Z2t.
(a)(b)(c)For system (4), we consider s as the bifurcation parameter under P0=6,Z10=2,andZ20=1. We select a set of parameters: r=1,K=20,μ1=5,μ2=6,α1=6,α2=10,β1=4,β2=5,ρ1=1,ρ2=1.400,g1=1,andg2=1.110. By simple calculation, we can know that s∗=2.199. Thus, E∗4.578,1.944,0.121 is locally asymptotically stable when s=1.500<s∗, see Figure 2. When s=11>s∗, E∗ is unstable, see Figure 3.Figure 2
The coexistence equilibriumE∗4.458,1.944,0.121 of system (4) is locally asymptotically stable when s=1.500<s∗=2.199. (a) Dynamical response of population. (b) Phase portrait.
(a)(b)Figure 3
The coexistence equilibriumE∗ of system (4) is destabilized through a Hopf bifurcation when s=11>s∗=2.199. (a) Dynamical response of population. (b) Phase portrait.
(a)(b)Next, we choose the values of parameters:r=1,K=20,μ1=5,μ2=6,α1=6,α2=10,β1=4,β2=5,ρ1=1,ρ2=1.400,g1=1,g2=1.110,ands=4. We can get ω10=3.433 and τ0=1.830 when τ≠0. According to Theorem 6, we obtain that the coexistence equilibrium E∗4.458,0.729,1.216 is locally asymptotically stable when τ=2<τ0, see Figure 4. Then, we select the values of τ as τ=3.500>τ0 and get that c10≈−1.573−8.542i<0,ν1≈0.458>0,ϱ1≈−3.145<0,andT1≈0.711>0 by using formula (51). From Theorem 2, we know that system (21) undergoes Hopf bifurcation and the bifurcation direction is forward when τ passes through its critical value τ0. In addition, system (21) has a stable bifurcation periodic solution, and the period of the bifurcating periodic solutions is increasing. From Figure 5, we know that the coexistence equilibrium E∗4.458,0.729,1.216 is destabilized through a Hopf bifurcation.Figure 4
The coexistence equilibriumE∗4.458,0.729,1.216 of model (21) is locally asymptotically stable when τ=1<τ0=1.830. (a) Dynamical response of population. (b) Phase portrait.
(a)(b)Figure 5
Whenτ=3.500>τ0=1.830, Hopf bifurcation takes place at the coexistence equilibrium E∗4.458,0.729,1.216 of system (21). (a) Dynamical response of population. (b) Phase portrait.
(a)(b)Finally, we choose a new set of parameter values to simulate system (3). d1=0.500,d2=0.500,d3=0.400,r=1,K=20,μ1=5,μ2=6,α1=6,α2=10,β1=4,β2=5,ρ1=1,ρ2=1.400,g1=1,g2=1.110,ands=6. From Theorem 7, we can get that the coexistence equilibrium E∗4.458,0.486,1.435 is locally asymptotically stable when τ=0, see Figure 6.Figure 6
The coexistence equilibriumE∗4.458,0.486,1.435 of system (3) is locally asymptotically stable when τ=0. (a) Pt; (b) Z1t; (c) Z2t.
(a)(b)(c)Through simple calculation, we obtain thatω20=0.348, and τ00=8.640. According to Theorem 10, we know that E∗ is locally asymptotically stable when τ=8<τ00=8.640, see Figure 7. Then, we choose the value of τ=15>τ00=8.640 and obtain that c20≈−213−345i<0, ν2≈6120.800>0,ϱ2≈−4253.900<0,andT2≈660.823>0. From Theorem 11, we know that system (3) undergoes a forward Hopf bifurcation when τ across its critical value τ00=8.640. System (3) has stable bifurcating periodic solutions, and the period of the bifurcating periodic solutions is increasing. In addition, we see from Figure 8 that the coexistence equilibrium E∗4.458,0.486,1.435 is destabilized by a Hopf bifurcation.Figure 7
The coexistence equilibriumE∗4.458,0.486,1.435 of system (3) is locally asymptotically stable when τ=8<τ00=8.640. (a) Pt; (b) Z1t; (c) Z2t.
(a)(b)(c)Figure 8
The coexistence equilibriumE∗4.458,0.486,1.435 of system (3) is destabilized through a Hopf bifurcation when τ=15>τ00=8.640. (a) Pt; (b) Z1t; (c) Z2t.
(a)(b)(c)


## 6. Conclusions and Discussions

In this paper, we have proposed a reaction-diffusion model with time delay and two kinds of functional response incorporating three plankton populations, where the functional responses of two zooplankton to phytoplankton are different. Because there is interference between some zooplankton in freshwater lakes or oceans, so we assume that the functional response of zooplanktonZ1t,x is Beddington–DeAnglis functional response. We verify that the degree of mutual interference s can produce Hopf bifurcation using numerical simulation; thus, we can change the value of this parameter s to control the explosion of plankton. Then, we consider time delay τ as bifurcation parameters to study the dynamics of systems (21) and (3). We discuss the existence of the Hopf bifurcation when τ passes through the critical values τ0 or τ00. In addition, we get that properties of Hopf bifurcation by utilizing the normal form theory and the center manifold theorem. By analysis, we know that systems (3) and (21) are stable at the coexistence equilibrium E∗ when τ<τ00 and τ<τ0, respectively. But systems (3) and (21) are unstable at the coexistence equilibrium E∗ when τ>τ00 and τ>τ0, respectively. These results indicate that if time delay τ required for mature phytoplankton to produce toxins is small enough, then the densities of phytoplankton P and zooplankton Z1 and Z2 will tend to constant; but if τ is too large, systems (3) and (21) will have periodic solutions, which means that the densities of phytoplankton P and zooplankton Z1 and Z2 will change periodically. Through analysis, we find that the size of the diffusion coefficient does not change the stability of the system.In order to make the model more realistic, we will consider linear harvesting of zooplankton, such asq1EZ1 and q1EZ1, where q1 and q2 are the catchability coefficients of zooplankton Zii=1and2, respectively, and E represents the harvesting effort. The model is as follows:(101)∂Pt,x∂t=d1ΔPt,x+rPt,x1−Pt,xK−μ1Pt,xZ1t,xα1+Pt,x+sZ1t,x−μ2Pt,xZ2t,xα2+Pt,x,t>0andx∈0,π,∂Z1t,x∂t=d2ΔZ1t,x+β1Pt,xZ1t,xα1+Pt,x+sZ1t,x−g1Z1t,x−q1EZ1t,x−ρ1Pt−τ,xZ1t,xα1+Pt−τ,x+sZ1t,x,t>0andx∈0,π,∂Z2t,x∂t=d3ΔZ2t,x+β2Pt,xZ2t,xα2+Pt,x−ρ2Pt−τ,xZ2t,xα2+Pt−τ,x−g2Z2t,x−q2EZ2t,x,t>0andx∈0,π,∂Pt,x∂x=∂Z1t,x∂x=∂Z2t,x∂x=0,t>0andx=0andπ,Pt,x=φ1t,x≥0,Z1t,x=φ2t,x≥0,Z2t,x=φ3t,x≥0,t∈−τ,0andx∈0,π.We leave this work for the future.



---

*Source: 5560157-2021-06-17.xml*