# Stability and Hopf Bifurcation in a Three-Component Planktonic Model with Spatial Diffusion and Time Delay

Authors: Kejun Zhuang; Gao Jia; Dezhi Liu
Journal: Complexity
 (2019)

---

## Abstract

Due to the different roles that nontoxic phytoplankton and toxin-producing phytoplankton play in the whole aquatic system, a delayed reaction-diffusion planktonic model under homogeneous Neumann boundary condition is investigated theoretically and numerically. This model describes the interactions between the zooplankton and two kinds of phytoplanktons. The long-time behavior of the model and existence of positive constant equilibrium solution are first discussed. Then, the stability of constant equilibrium solution and occurrence of Hopf bifurcation are detailed and analyzed by using the bifurcation theory. Moreover, the formulas for determining the bifurcation direction and stability of spatially bifurcating solutions are derived. Finally, some numerical simulations are performed to verify the appearance of the spatially homogeneous and nonhomogeneous periodic solutions.


---

## Body

## 1. Introduction

Oceans have a major role in the global carbon cycling and so directly impact the pace and extent of climate change [1]. Marine organism can bring great economic and social values. It has the irreplaceable function in the global food processing, tourism, nutrient cycling, gas regulation, and so on. As the basis of the marine food chain system, plankton can supply food and oxygen to a myriad of marine life and can also absorb about half of the climate-warming carbon dioxide [2]. Besides, unlike fish or some intertidal creatures, plankton has rarely been commercially exploited. Moreover, this free-floating plankton can respond quickly to the temperature and change of oceanic system. As a consequence, mathematical modeling is a valuable tool for research fields of the marine ecology.In general, plankton can be broadly divided into autotrophic phytoplankton and herbivorous zooplankton. The interactive process between phytoplankton and zooplankton is very complicated. It depends not only on the foraging style and feeding behavior, but also on other processes in the system. To describe the interactive relationship between phytoplankton and zooplankton, the following general model can be utilized:(1)dPdt=GP-FP,Z-MP,dZdt=FP,Z-MZ.Here, P(t) and Z(t) denote the concentration of phytoplankton and zooplankton respectively, G(P) and M(P) represent the growth rate and death rate of phytoplankton, respectively, M(Z) is the death rate of zooplankton, and F(P,Z) is the predation rate of zooplankton on phytoplankton. Recently, various particular cases of model (1) have been well studied [3–5]. Some interesting results about the stability and Hopf bifurcation type periodic oscillations have been obtained.It should be noted that some phytoplankton can release toxic substances which will result in poisoning in both fish and shellfish. By accumulating in marine food webs, the toxins may have hazards on animal and human health. On the other hand, due to the rapid growth of plankton, harmful algal bloom may cause massive death of marine animals. However, most previous studies did not directly consider the toxin-producing phytoplankton population and the role of toxic phytoplankton can not be ignored. While in the establishment of planktonic mathematical models, it is necessary to include the nontoxic phytoplankton, toxin-producing phytoplankton, and zooplankton. In 2004, Chattopadhyay et al. [6] proposed the following kinetic model which has three interacting components and is with an additional factor that the release of toxic substance reduces the growth of zooplankton: (2)dP1dt=rP11-P1K-αP1Z,dP2dt=sP21-P2K-θP2Zγ+P2,dZdt=βP1Z-μZ-θ1P1Zγ+P2,where P1(t) denotes the concentration of the nontoxic phytoplankton at time t and P2(t) and Z(t) denote the concentration of toxic phytoplankton population and zooplankton at time t, respectively. It is assumed that the two phytoplankton populations share the same resource. In model (2), all the coefficients are positive constants, r and s are the growth rates of two phytoplankton populations, respectively, K is the carrying capacity, α and β are the maximum zooplankton ingestion rate and maximum zooplankton conversion rate, respectively, μ is the death rate of zooplankton, θ is the rate of toxin liberation by toxic phytoplankton, and θ1 is the specific predation rate of zooplankton population on toxic phytoplankton. This model shows that toxic substances released by phytoplankton have negative effects on the grazing pressure of zooplankton. It is finally concluded that toxin-producing phytoplankton may be used as a biocontrol agent for the harmful algal bloom problems. It is also mentioned that the role of time delay and environmental fluctuation in the planktonic dynamics may arouse some interesting results and needs further investigations.Motivated by (2), some modified models have been proposed recently. For instance, Sarkar et al. [7] established a new model made of two harmful phytoplankton populations and one zooplankton population. Roy et al. [8] investigated the model where the two phytoplankton populations compete with each other. Pal et al. [9] further considered the three-component model with both nonlinear predation functions by zooplankton. Further, the study was also extended from the perspectives of stochastic dynamics and plankton-nutrient interactions, respectively [10, 11].In view of the ocean current and monsoon, the plankton can freely drift and this spatial dispersal is subject to Fickian diffusion. So the effect of spatial diffusion has been investigated by many authors [12–17]. The results indicate that spatial diffusion has a vital role in the spatiotemporal dynamics of the planktonic model and spatial pattern may occur. Besides, the impact of time delay can not be ignored because it usually causes periodic oscillations, even chaotic behaviors, and time delay is ubiquitous in the real ecosystem [18–20].According to the above factors, we consider the following three-component planktonic model with spatial diffusion and time delay:(3)∂P∂t=d1ΔP+r1P1-PK-αPZ,x,t∈0,L×0,+∞,∂T∂t=d2ΔT+r2T1-TK-aTZγ+T,x,t∈0,L×0,+∞,∂Z∂t=d3ΔZ+βPZ-dZ2-bTx,t-τZγ+Tx,t-τ,x,t∈0,L×0,+∞,∂P∂x=∂T∂x=∂Z∂x=0,t≥0,x=0,L,Px,t=P0x,t≥0,Tx,t=T0x,t≥0,Zx,t=Z0x,t≥0,x,t∈0,L×-τ,0,where P(x,t), T(x,t), and Z(x,t) denote the densities of nontoxic phytoplankton, toxin-producing phytoplankton, and zooplankton at location x and time t, respectively, Δ is the usual Laplace operator, L denotes the depth of the water column, and the homogeneous Neumann boundary condition means that no plankton species is entering or leaving the column at the top or the bottom.All the parameters are positive constants,di, i=1,2,3, are the three species’ diffusion rates, respectively, d is the higher mortality of zooplankton, and τ is the time needed for zooplankton from ingesting toxic phytoplankton to dying. The other coefficients have the same meanings as in model (2). Note that the zooplankton may get eaten by higher predators, whose population is not being explicitly modelled [21–23]. So, we adopt the quadratic closure term dZ2 to describe the higher mortality of zooplankton.In this paper, we mainly investigate the spatiotemporal dynamics of delayed and diffusive system (3). The rest of the paper is organized as follows. In Section 2, the permanence and nonpersistence of system (3) are derived. In Section 3, the sufficient conditions for existence of positive constant equilibrium solution are obtained. In Section 4, the stability of equilibrium solution and delay-induced Hopf bifurcation are explored. In Section 5, the detailed formulae for determining the bifurcation properties are given by calculating the normal form on the center manifold. In Section 6, some numerical simulations are conducted to illustrate the theoretical results. Finally, some conclusions are given in Section 7.


## 2. Long-Time Behavior

In this section, we shall show that any nonnegative solutionP(x,t),T(x,t),Z(x,t) of system (3) lies in a bounded region as t→∞, for all x∈Ω.Theorem 1.
Ifdr1>αβK and dr1βγ-αβ2γK-bdr1>0 hold, then system (3) is permanent; that is, there exist positive constants m and M independent of solution such that (4)m≤liminft→+∞Px,t≤limsupt→+∞Px,t≤M,m≤liminft→+∞Tx,t≤limsupt→+∞Tx,t≤M,m≤liminft→+∞Zx,t≤limsupt→+∞Zx,t≤Mfor any nonnegative solution.Proof.
From the first equation of system (3), we have (5)∂P∂t-d1ΔP≤r1P1-PK.The standard comparison principle implies (6)limsupt→+∞Px,t≤K,and thus for every real number ε1>0, there exists a T1>0 such that P(x,t)≤K+ε1, for all t>T1.
Similarly, from the second equation of (3), we have (7)limsupt→+∞Tx,t≤K,And thus for every real number ε2>0, there exists a T2>0 such that T(x,t)≤K+ε2, for all t>T2.
The third equation of (3) can be reduced to (8)∂Z∂t-d3ΔZ≤βK+ε1Z-dZ2=ZβK+ε1-dZ,for all t>T1. Then, we have Z(x,t)≤β(K+ε1)/d, which indicates (9)limsupt→+∞Zx,t≤βKdas ε1→0. Therefore, for every real number ε3>0, there exists a T3>0 such that Z(x,t)≤βK/d+ε3, for all t>T3.
Again from the first equation of (3), we have (10)∂P∂t-d1ΔP≥r1P1-αr1βKd+ε3-PK,for all t>T3. This implies (11)Px,t≥1-α/r1βK/d+ε31/K;thus, (12)liminft→+∞Px,t≥Kdr1-αβKdr1>0.Then, for every real number ε4>0, there exists a T4>0 such that P(x,t)≥K(dr1-αβK)/dr1-ε4, for all t>T4.
From the second equation of (3), we have (13)Px,t≥Kdr1-αβKdr1-ε4=r2T1-1K+ar2γβKd+ε3T,which implies (14)liminft→+∞Tx,t≥dr2γKdr2γ+aβK2>0.
From the third equation of (3), we also have (15)∂Z∂t-d3ΔZ≥ZβKdr1-αβKdr1-ε4-dZ-bK+ε2γ,for all t>maxT2+τ,T4. Thus, (16)liminft→+∞Zx,t≥dr1βγK-αβ2γK2-bdr1Kd2r1γ>0.
Finally, if we set(17)m=minKdr1-αβKdr1,dr2γKdr2γ+aβK2,dr1βγK-αβ2γK2-bdr1Kd2r1γ,and (18)M=maxK,βKd,then the proof is complete.Definition 2.
System (3) is said to be not persistent if (19)minliminft→+∞Px,t,liminft→+∞Tx,t,liminft→+∞Zx,t=0for some of its nonnegative solutions.Next, we discuss the nonpersistence of system (3).Theorem 3.
Ifbdr2γ≥β(γ+K)(dr2γ+aβK2) holds, then system (3) is not persistent.Proof.
According to the process of Theorem1, for an arbitrary positive constant ε5, there exists a T5>0 such that (20)Tx,t≥dr2γKdr2γ+aβK2-ε5,for all t>T5. Further, we have (21)∂Z∂t-d3ΔZ≤ZβK+ε1-bdr2γK-bε5dr2γ+aβK2γ+K+ε2dr2γ+aβK2-dZ,for all t>maxT1,T2+τ,T5+τ. Then, by the arbitrariness of ε1, ε2, and ε5, we have (22)liminft→+∞Zx,t≤limsupt→+∞Zx,t≤βKd-br2γKγ+Kdr2γ+aβK2≤0as t→+∞. The proof is complete.


## 3. Steady State

In consideration of the biological significance of (3), we focus on the positive constant equilibrium solution. To determine the equilibrium solution of (3), we only need to solve the following algebraic equations: (23)r1P1-PK-αPZ=0,r2T1-TK-aTZγ+T=0,βPZ-dZ2-bTZγ+T=0.Simplifying the first equation of (23) and substituting it into the third equation, we have (24)Z=r1K-PαK,and (25)β+dr1αKP-dr1α-bTγ+T=0.From (24) and the second equation of (23), we have (26)aZγ+T=ar1αγ+T-ar1PαKγ+T=r2-r2KT,and (27)ar1PαKγ+T=ar1αγ+T-r2+r2KT;thus, (28)P=K-αr2KT+γar1+αr2TT+γar1,which is equivalent to (29)P=αr2ar1T-K-γ22+ar1Kαr2-γK-K-γ24.From (25) and (28), we have (30)β+dr1αKK-αr2KT+γar1+αr2TT+γar1-dr1α-bTγ+T=0,and (31)a3T3+a2T2+a1T+a0=0,where(32)a3=αr2ar1β+dr1αK>0,a2=αr2ar12γ-Kβ+dr1αK,a1=β+dr1αKK+αγr2ar1γ-2K-b-dr1α,a0=γKβ+dr1αK1-αγr2ar1-dr1γα.According to Descartes’ rule of signs, cubic equation (31) has at least one positive real root when a0<0, that is, (33)Kβ+dr1αK1-αγr2ar1<dr1α.For convenience, we denote any positive root of (31) by T∗. Combining (24) and (29), we can obtain the positive solution of (23) under the conditions 4ar1K/αr2>(K+γ)2 and K>P∗.From above analyses, we can establish the existence of positive constant equilibrium solution of (3).Theorem 4.
If the following assumption(H1)
Kβ+dr1/αK1-αγr2/ar1<dr1/α,  4ar1K/αr2>K+γ2and K>P∗
holds, then system (3) has positive constant equilibrium solution E∗=(P∗,T∗,Z∗).Based on the aim of this study, we always assume that condition (H1) is satisfied in the following sections.


## 4. Hopf Bifurcation Induced by Time Delay

Here, we will regard time delayτ as the bifurcation parameter to investigate its effect on the stability of coexistence equilibrium solution E∗=(P∗,T∗,Z∗).Linearizing system (3) at E∗ leads to the corresponding characteristic equation: (34)λ+d1kπL+r1KP∗0αP∗0λ+d2kπL+r2T∗K-aT∗Z∗γ+T∗2aT∗γ+T∗-βZ∗bγZ∗γ+T∗2e-λτλ+d3kπL+dZ∗=0,which can be simplified as the following transcendental equation: (35)λ3+Akλ2+Bkλ+Ck+Dkλ+Fke-λτ=0,where (36)Ak=d1+d2+d3kπL2+r1KP∗+r2KT∗+dZ∗-aT∗Z∗γ+T∗2,Bk=αβP∗Z∗+d1kπL2+r1KP∗d2kπL2+r2KT∗-aT∗Z∗γ+T∗2+d1kπL2+r1KP∗d3kπL2+dZ∗+d2kπL2+r2KT∗-aT∗Z∗γ+T∗2d3kπL2+dZ∗,Ck=d2kπL2+d3kπL2+r2KT∗-aT∗Z∗γ+T∗2+dZ∗·d1kπL2+r1KP∗+αβP∗Z∗d2kπL2+r2KT∗-aT∗Z∗γ+T∗2,Dk=-abγT∗Z∗γ+T∗3<0,Fk=-abγT∗Z∗γ+T∗3d1kπL2+r1KP∗<0.The special case of (35) with τ=0 is (37)λ3+Akλ2+Bk+Dkλ+Ck+Fk=0.If all the roots of (37) have negative real parts for every nonnegative integer k, then the positive equilibrium solution E∗ without time delay is asymptotically stable. With the help of Routh-Hurwitz criterion, E∗ is stable if and only if Ak>0, Bk+Dk>0, Ck+Fk>0, and Ak(Bk+Dk)-(Ck+Fk)>0.Assume that(H2)
r2/K>aZ∗/γ+P2∗2, and abγT∗/γ+T∗3<mind,αβP∗;then, we have the following stability conclusion.Theorem 5.
If assumption (H2) is satisfied, then the equilibrium solutionE∗ of (3) is asymptotically stable for τ=0.Next, we shall discuss the distribution of characteristic roots whenτ>0. Suppose iω(ω>0) is a root of (35). Then, for any nonnegative integer k, we have (38)-iω3-Akω2+iBkω+Ck+iDkω+Fkcosωτ-isinωτ=0,and(39)-iω3-Akω2+iBkω+Ck+iDkωcosωτ+Dkωsinωτ+Fkcosωτ-iFksinωτ=0.Separating the real and imaginary parts results in (40)ω3-Bkω=Dkωcosωτ-Fksinωτ,Akω2-Ck=Fkcosωτ+Dkωsinωτ;thus, (41)ω6+Ak2-2Bkω4+Bk2-2AkCk-Dk2ω2+Ck2-Fk2=0.For simplicity, we set σ=ω2; then, (41) can be rewritten in the form of (42)Rσ=σ3+Mkσ2+Nkσ+Pk=0,where (43)Mk=Ak2-2Bk,Nk=Bk2-2AkCk-Dk2,Pk=Ck2-Fk2.If the assumption(H3)
Mk>0 and Mk>0for  any  k>0holds, then we havePk>0 by combining with (H2). And the cubic equation (42) of σ has no positive root, so (35) has no purely imaginary root. It can be concluded that the equilibrium solution E∗ is always asymptotically stable for any τ>0 and that system (3) has no spatially nonhomogeneous periodic solution.On the other hand, whenk=0, based on Lemma 2.2 in [24], the following condition is needed to ensure the existence of positive root of (42):(H4)
σ1>0 and R(σ1)k=0≤0,  whereσ1=-M0+M02-3N0/3.Denote any positive root of (42) by σ∗; then, ±iω∗ is a pair of purely imaginary roots of (35), where ω∗=σ∗.Whenk=0, rewrite (35) as (44)λ3+A0λ2+B0λ+C0+D0λ+F0e-λτ=0,and we have (45)sinω∗τ=A0D0-F0ω∗3+B0F0-C0D0ω∗D02ω∗2+F02≜Fs,(46)cosω∗τ=D0ω∗4+A0F0-B0D0ω∗2-C0F0D02ω∗2+F02≜Fc.Define(47)τ0j=1ω∗2π-arccosFc+2jπ,Fs≥0,1ω∗arccosFc+2jπ,Fc<0,where j=0,1,2,… and τ0=τ00=minj=0,1,2,…τ0j. Then, (44) has a pair of purely imaginary roots ±iω∗ when τ=τ0j.We claim that if(H5)
R′(σ∗)≠0,then(48)dReλτdττ=τ0j,ω=ω∗≠0.In fact, differentiating both sides of (35) with respect to τ, it follows that (49)dλdτ-1=3λ2+2A0λ+B0eλτ+D0λD0λ+F0-τλ.Thus, we have (50)Redλdτ-1τ=τ0j,ω=ω∗=ω∗2R′σ∗D0ω∗22+F0ω∗2≠0.From what has been discussed above and the Hopf bifurcation Theorem by Hassard et al. [25], we can draw the conclusion on the existence of spatially homogeneous Hopf bifurcation.Theorem 6.
Suppose that conditions (H1)-(H5) are satisfied.
(i) The equilibrium solutionE∗ is locally asymptotically stable for τ∈[0,τ0).
(ii) The equilibrium solutionE∗ is unstable for τ>τ0 and τ0 is the Hopf bifurcation value.
(ii) System (3) undergoes spatially homogeneous periodic solutions at E∗ when τ=τ0.


## 5. Stability and Direction of the Bifurcation

In this section, we investigate the properties of spatially homogeneous periodic solutions, including bifurcation direction, stability of periodic solutions, monotonicity of periodic solutions. Here, we mainly apply the normal form theory and center manifold theorem for partial functional differential equations [25, 26].For fixedj=0,1,2,…, denote bifurcation value τ0j by τ∗ and introduce the new parameter μ=τ-τ∗; then, μ=0 is the new Hopf bifurcation value. Let ϑ=t/τ and rewrite ϑ as t; system (3) can be transformed into (51)dUtdt=τ∗DΔUt+Lτ∗Ut+FUt,μ,where L(μ)(φ):C→X and F(·,μ):C→X are given by(52)φ=φ1,φ2,φ3T∈C,D=diagd1,d2,d3,Δ=diag∂2∂x2,∂2∂x2,∂2∂x2,Utθ=Ut+θ,-τ≤θ≤0,Lμφ=μ-r1KP∗φ10-αP∗φ30aT∗Z∗γ+T∗2-r2KT∗φ20-aT∗γ+T∗φ30βZ∗φ10-bγZ∗γ+T∗2φ2-1-dZ∗φ30,Fφ,μ=μDΔφ0+Lμφ+fφ,μ,and (53)fφ,μ=τ∗+μ-2r1Kφ120-αφ10φ302aγZ∗γ+T∗3φ220-2r2Kφ220-aγγ+T∗2φ20φ30βφ10φ30-2dφ320-bγφ2-1φ30γ+T∗2+2bγZ∗γ+T∗3φ22-1+h.o.t.Then, the linearized system of (51) at origin is (54)dUtdt=τ∗DΔUt+Lτ∗Ut.From the discussion in Section 4, we can find that characteristic equation (35) has a pair of purely imaginary roots Λ0=iω∗τ∗,-iω∗τ∗ when τ=τ∗.LetC:=C([-1,0],R3). Then, we consider the functional differential equation on C: (55)z˙=Lτ∗zt.It is obvious that L(τ∗) is a continuous linear function mapping C([-1,0],R3) into R3. By the Riesz representation theorem, there exists a 3×3 matrix function η(θ,τ) (-1≤θ≤0) such that (56)Lτ∗φ=∫-10dηθ,τ∗φθ,forφ∈C;here we choose (57)ηθ,τ∗=τ∗-r1KP∗0-αP∗0aT∗Z∗γ+T∗2-r2KT∗-aT∗γ+T∗βZ∗0-dZ∗δθ-τ∗0000000-bγZ∗γ+T∗20δθ+1.LetA(τ∗) be the infinitesimal generator of the semigroup induced by the solutions of (55) and A∗ be the adjoint matrix of A(τ∗) under the bilinear pairing: (58)ψ,ϕ=ψ0,ϕ0-∫-10∫ξ=0θψξ-θdηθϕξdξ=ψ0,ϕ0+τ∗∫-10ψθ+10000000-bγZ∗γ+P2∗20ϕθdθ,where φ∈C, ψ∈C∗=C([0,1],R3). Then, A(τ∗) and A∗ are a pair of adjoint operators and they both have characteristic roots ±iω∗τ∗. Let P and P∗ be the generalized eigenspaces of A(τ∗) and A∗, respectively; then, P∗ is the adjoint space of P and dimP=dimP∗=2.By direct calculation, we have the following lemma.Lemma 7.
Let(59)ξ=aT∗ζγ+T∗aKT∗Z∗-r2T∗γ+T∗2-iω∗Kγ+T∗2,ζ=-iω∗K+r1P∗αKP∗,ξ∗=-aKT∗Z∗+bKγζ∗Z∗e-iω∗τ∗r2T∗γ+T∗2,ζ∗=iω∗K+r1P∗βKZ∗;then a basis of P with Λ0 is (60)p1θ=eiω∗τ∗θ1,ξ,ζT,p2θ=p1θ¯,-1≤θ≤0,and a basis a of P∗ with Λ0 is (61)q1s=e-iω∗τ∗s1,ξ∗,ζ∗T,q2s=q1s¯,0≤s≤1.LetΦ=(Φ1,Φ2), Ψ∗=(Ψ1∗,Ψ2∗)T, where (62)Φ1θ=p1θ+p2θ2,Φ2θ=p1θ-p2θ2i,θ∈-1,0,and (63)Ψ1∗s=q1s+q2s2,Ψ2∗s=q1s-q2s2i,s∈0,1.From (58), we can obtain (Ψ1∗,Φ2), (Ψ1∗,Φ1), and further (64)q1,p1=Ψ1∗,Φ1-Ψ2∗,Φ2+iΨ1∗,Φ2-Ψ2∗,Φ1;that is, (65)q1,p1=1+ξξ∗+ζζ∗-τ∗bZ∗γ+P2∗ξ∗e-iω∗τ∗≜D∗.Therefore, we have(66)Ψ1∗,Φ1-Ψ2∗,Φ2=ReD∗,Ψ1∗,Φ2-Ψ2∗,Φ1=ImD∗.Now, we define(Ψ∗,Φ)=(Ψl∗,Φm) (l,m=1,2) and construct a new basis for P∗ by Ψ=Ψ1,Ψ2T=Ψ∗,Φ-1Ψ∗. Moreover, we define f0=(ξ01,ξ02,ξ03), where (67)ξ01=100,ξ02=010,ξ03=001.Letc·f0 be defined by (68)c·f0=c1ξ01+c2ξ02+c3ξ03for c=(c1,c2,c3)T, cj∈R(j=1,2,3). Then, the center space of linear equation (54) is given by PCNC, where (69)PCNφ=ΦΨ,φ,f0·f0,φ∈C,and C=PCNC⊕PSC, PSC denotes the complementary subspace of PCNC.If(70)Aτ∗φθ=φ˙θ+X0θτ∗DΔφ0+Lτ∗φθ-φ˙0,φ∈BC,where X0:[-1,0]→B(X,X) and X0(θ)=0,-1≤θ≤0,I,θ=0. then Aτ∗ is the infinitesimal generator induced by the solutions of (51) and (54), which can be written as the following operator differential equation: (71)U˙t=Aτ∗Ut+X0FUt,μ,Therefore, the solution of (51) can be written in the form of (72)Ut=Φx1tx2t·f0+hx1,x2,μ,x1tx2t=Ψ,Ut,f0,where h(x1,x2,μ)∈PSC, h(0,0,0)=Dh(0,0,0)=0. Specifically, the solution of (51) on the center manifold is (73)Ut∗=Φx1tx2t·f0+hx1,x2,0.Let z=x1-ix2, and notice that p1=Φ1+iΦ2; then, we can rewrite (73) as (74)Ut∗=12Φz+z¯iz-z¯·f0+Wz,z¯=12p1z+p¯1z¯·f0+Wz,z¯,where W(z,z¯)=h(z+z¯/2,-z-z¯/2i,0). Furthermore, from [26], z also satisfies (75)z˙=iω∗τ∗z+gz,z¯,where (76)gz,z¯=Ψ10-iΨ20FUt∗,0,f0.Let (77)Wz,z¯=W20z22+W11zz¯+W02z¯22+⋯and (78)gz,z¯=g20z22+g11zz¯+g02z¯22+⋯.By (74), it is not difficult to compute that (79)
F
U
t
∗
,
0
,
f
0
=
τ
∗
z
2
4
-
r
1
K
-
α
ζ
a
γ
Z
∗
γ
+
T
∗
3
ξ
2
-
r
2
K
ξ
2
-
a
γ
γ
+
T
∗
2
ξ
ζ
β
ζ
-
2
d
ζ
2
-
2
b
γ
γ
+
T
∗
2
ξ
ζ
e
-
i
ω
∗
τ
∗
+
2
b
γ
Z
∗
γ
+
T
∗
3
ζ
2
e
-
2
i
ω
∗
τ
∗
+
τ
∗
z
z
¯
4
-
2
r
1
K
-
α
ζ
+
ζ
¯
2
a
γ
Z
∗
γ
+
T
∗
3
ξ
ξ
¯
-
2
r
2
K
ξ
ξ
¯
-
a
γ
γ
+
T
∗
2
ξ
ζ
¯
+
ξ
¯
ζ
β
ζ
+
ζ
¯
-
2
d
ζ
ζ
¯
-
b
γ
γ
+
T
∗
2
ξ
ζ
¯
e
-
i
ω
∗
τ
∗
+
ξ
¯
ζ
e
i
ω
∗
τ
∗
+
4
b
γ
Z
∗
γ
+
T
∗
3
ζ
ζ
¯
+
τ
∗
z
¯
2
4
-
r
1
K
-
α
ζ
¯
a
γ
Z
∗
γ
+
T
∗
3
ξ
¯
2
-
r
2
K
ξ
¯
2
-
a
γ
γ
+
T
∗
2
ξ
¯
ζ
¯
β
ζ
¯
-
2
d
ζ
¯
2
-
2
b
γ
γ
+
T
∗
2
ξ
¯
ζ
¯
e
i
ω
∗
τ
∗
+
2
b
γ
Z
∗
γ
+
T
∗
3
ζ
¯
2
+
τ
∗
z
2
z
¯
4
-
r
1
K
4
W
11
1
0
+
2
W
20
1
0
-
α
W
20
3
0
+
2
W
11
1
0
+
2
ζ
W
11
1
0
+
ζ
¯
W
20
1
0
,
1
a
γ
Z
∗
γ
+
T
∗
3
-
r
2
K
2
ξ
¯
W
20
2
0
+
4
ξ
W
11
2
0
-
a
γ
γ
+
T
∗
2
2
ξ
W
11
3
0
+
ξ
¯
W
20
3
0
+
2
ζ
W
11
2
0
+
ζ
¯
W
20
2
0
,
1
-
b
γ
γ
+
T
∗
2
2
ξ
e
-
i
ω
∗
τ
∗
W
11
3
0
+
ξ
¯
e
i
ω
∗
τ
∗
W
20
3
0
+
2
ζ
W
11
2
-
1
+
ζ
¯
W
20
2
-
1
-
2
d
2
ζ
¯
W
20
3
0
+
4
ζ
W
11
3
0
+
β
W
20
3
0
+
2
W
11
1
0
+
2
ζ
W
11
1
0
+
ζ
¯
W
20
1
0
+
2
b
γ
Z
∗
γ
+
T
∗
3
2
ξ
e
-
i
ω
∗
τ
∗
W
11
2
-
1
+
ξ
¯
e
i
ω
∗
τ
∗
W
20
2
-
1
,
1
+
⋯
,where Wij(n)(θ),1=1/π∫0π‍Wij(n)(θ)(x)dx, i+j=2, n=1,2,3. Let (ψ1,ψ2,ψ3)=Ψ1(0)-iΨ2(0). We can get the following expressions: (80)g20=τ∗2-r1K-αζψ1+aγZ∗γ+T∗3ξ2-r2Kξ2-aγγ+T∗2ξζψ2+βζ-2dζ2-2bγγ+T∗2ξζe-iω∗τ∗+2bγZ∗γ+T∗3ζ2e-2iω∗τ∗ψ3,g11=τ∗4-2r1K+αζ+ζ¯ψ1+2aγZ∗γ+T∗3ξξ¯-2r2Kξξ¯-aγγ+T∗2ξζ¯+ξ¯ζψ2+βζ+ζ¯-2dζζ¯-bγγ+T∗2ξζ¯e-iω∗τ∗+ξ¯ζeiω∗τ∗+4bγZ∗γ+T∗3ζζ¯ψ3,g02=g20¯,g21=τ∗2-r1K4W1110+2W2010-αW2030+21+ζW1110+ζ¯W2010,1ψ1+aγZ∗γ+T∗3-r2K2ξ¯W2020+4ξW1120-aγγ+T∗22ξW1130+ξ¯W2030+2ζW1120+ζ¯W2020,1ψ2+-bγγ+T∗22ξe-iω∗τ∗W1130+ξ¯eiω∗τ∗W2030+2ζW112-1+ζ¯W202-1-2d2ζ¯W2030+4ζW1130+βW2030+2W1110+2ζW1110+ζ¯W2010+2bγZ∗γ+T∗32ξe-iω∗τ∗W112-1+ξ¯eiω∗τ∗W202-1,1ψ3.Since W20(θ) and W11(θ) appear in g21, we should also establish them. According to (77), we have (81)W˙z,z¯=W20zz˙+W11z˙z¯+zz˙¯+W02z¯z˙¯+⋯,and (82)Aτ∗W=Aτ∗W20z22+Aτ∗W11zz¯+Aτ∗W02z¯22+⋯.In addition, by [26], W(z(t),z¯(t)) also satisfies (83)W˙=Aτ∗W+Hz,z¯,where (84)Hz,z¯=H20z22+H11zz¯+H02z¯22+⋯=X0FUt∗,0-ΦΨ,X0FUt∗,0,f0·f0with Hij∈PSC and i+j=2.Consequently, (74) and (81)-(83) lead to (85)2iω∗τ∗-Aτ∗W20=H20,-Aτ∗W11=H11.Since Aτ∗ has a pair of purely imaginary characteristic roots, (83) has the unique solution such that (86)W20=2iω∗τ∗-Aτ∗-1H20,W11=-Aτ∗-1H11.For θ∈[-1,0], it follows from (84) that (87)Hz,z¯=-ΦθΨθFUt∗,0,f0·f0=-p1θ+p2θ2,p1θ-p2θ2iΨ10,Ψ20T×FUt∗,0,f0·f0=-12p1θΨ10-iΨ20+p2θΨ10+iΨ20×FUt∗,0,f0·f0=-14g20p1θ+g¯02p2θz2·f0-12g11p1θ+g¯11p2θzz¯·f0+⋯.Thus, for θ∈[-1,0], we have (88)H20θ=-12g20p1θ+g¯02p2θ·f0,H11θ=-12g11p1θ+g¯11p2θ·f0,Hz,z¯0=FUt∗,0-ΦΨ,FUt∗,0,f0·f0,H200=τ∗2-r1K-αζaγZ∗γ+T∗3ξ2-r2Kξ2-aγγ+T∗2ξζβζ-2dζ2-2bγξζe-iω∗τ∗γ+T∗2+2bγZ∗ζ2e-2iω∗τ∗γ+T∗3-g20p10+g¯02p202·f0,and (89)H110=τ∗4-2r1K-αζ+ζ¯2aγZ∗γ+T∗3ξξ¯-2r2Kξξ¯-aγγ+T∗2ξζ¯+ξ¯ζβζ+ζ¯-2dζζ¯-bγγ+T∗2ξζ¯e-iω∗τ∗+ξ¯ζeiω∗τ∗+4bγZ∗γ+T∗3ζζ¯-g11p10+g¯11p202·f0.Combining with the definition of Aτ∗ and (86), we have (90)W˙20θ=2iω∗τ∗W20θ+g20p10+g¯02p202·f0,-1≤θ≤0.Since p1(θ)=p1(0)eiω∗τ∗, we have (91)W20θ=i2g20ω∗τ∗p1θ+g¯023ω∗τ∗p2θ·f0+Ee2iω∗τ∗θ,where (92)E=W200-i2g20ω∗τ∗p10+g¯023ω∗τ∗p20·f0.From the definition of Aτ∗ and (86) and (92), we have (93)2iω∗τ∗ig202ω∗τ∗p10·f0+ig¯206ω∗τ∗p20·f0+E-τ∗DΔig202ω∗τ∗p10·f0+ig¯206ω∗τ∗p20·f0+E-Lτ∗ig202ω∗τ∗p10·f0+ig¯206ω∗τ∗p20·f0+Ee2iω∗τ∗θ=τ∗2-r1K-αζaγZ∗γ+T∗3ξ2-r2Kξ2-aγγ+T∗2ξζβζ-2dζ2-2bγξζe-iω∗τ∗γ+T∗2+2bγZ∗ζ2e-2iω∗τ∗γ+T∗3-g20p10+g¯02p202·f0.Noticing that (94)τ∗DΔp10·f0+Lτ∗p1θ·f0=iω∗τ∗p10·f0,τ∗DΔp20·f0+Lτ∗p2θ·f0=-iω∗τ∗p20·f0,we have (95)2iω∗τ∗E-τ∗DΔE-Lτ∗Ee2iω∗τ∗θ=τ∗2-r1K-αζaγZ∗γ+T∗3ξ2-r2Kξ2-aγγ+T∗2ξζβζ-2dζ2-2bγξζe-iω∗τ∗γ+T∗2+2bγZ∗ζ2e-2iω∗τ∗γ+T∗3.From previous formulas, we can obtain (96)E=122iω∗τ∗+r1KP∗0αP∗02iω∗τ∗-aT∗Z∗γ+T∗2+r2KP2∗aT∗γ+T∗-βZ∗bZ∗γ+T∗e-2iω∗τ∗-bT∗Z∗γ+T∗22iω∗τ∗+dZ∗-1×-r1K-αζaγZ∗γ+T∗3ξ2-r2Kξ2-aγγ+T∗2ξζβζ-2dζ2-2bγξζe-iω∗τ∗γ+T∗2+2bγZ∗ζ2e-2iω∗τ∗γ+T∗3.With the same method, we also have(97)W˙11θ=g11p10+g¯11p202·f0,-1≤θ≤0,and (98)W11θ=i2ω∗τ∗-g11p1θ+g¯11p2θ·f0+F,where (99)F=14r1KP∗0αP∗0r2KT∗-aT∗Z∗γ+T∗2aT∗γ+T∗-βZ∗bγZ∗γ+T∗2dZ∗-1×-2r1K-αζ+ζ¯2aγZ∗γ+T∗3ξξ¯-2r2Kξξ¯-aγγ+T∗2ξζ¯+ξ¯ζβζ+ζ¯-2dζζ¯-bγγ+T∗2ξζ¯e-iω∗τ∗+ξ¯ζeiω∗τ∗+4bγZ∗γ+T∗3ζζ¯.At this point, we are able to completely establish the value of g21 and thereby can determine the properties of Hopf bifurcation.System (3) has the following Poincaré normal form: (100)ϖ˙=iω∗τ∗ϖ+C10ϖϖ2+oϖ5,where (101)C10=i2ω∗τ∗g20g11-2g112-g0223+g213.Hence, (102)σ2=-ReC10Reλ′τ∗,β2=2ReC10,T2=-ImC10+σ2Imλ′τ∗ω∗τ∗.By the Hopf bifurcation theory [25], we know that σ2 determines the bifurcation direction: if σ2>0 (σ2<0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ>τ0 (τ<τ0); β2 determines the stability of bifurcating periodic solutions: if β2>0 (β2<0), then the periodic solutions are stable (unstable); T2 determines the monotonicity of the period of periodic solutions: if T2>0 (T2<0), then the period increases (decreases).


## 6. Numerical Simulation

In this section, we conduct the numerical simulations with the help of MATLAB. We first choose the following parameter value:(103)d1=d2=0.5,d3=2,Ω=0,π,r1=5,r2=0.2,K=10,α=0.068,β=0.08,a=0.02,b=0.005,γ=0.5,d=0.7.By direct calculation, we have that system (3) has the unique positive equilibrium solution E∗=(9.84786,9.78241,1.11867) and the corresponding Hopf bifurcation value is τ0≈22.3146.Figures1–3 show that the constant equilibrium solution of (3) is asymptotically stable when the time delay is zero or appropriately small. On the other hand, once the time delay is larger than the critical value τ0, the equilibrium solution would no longer be stable and spatially homogeneous periodic solution will bifurcate at the equilibrium solution (see Figure 4).Figure 1
The equilibrium solutionE∗ of (3) is asymptotically stable when τ=0.
(a) (b) (c)Figure 2
The equilibrium solutionE∗ of (3) is asymptotically stable when τ=0.65<τ0.
(a) (b) (c)Figure 3
The equilibrium solutionE∗ of (3) is asymptotically stable when τ=15<τ0.
(a) (b) (c)Figure 4
The spatially homogeneous periodic solution exists whenτ=25>τ0.
(a) (b) (c)We reselectd1=d2=0.0005, d3=0.005, d=0.1, and keep other coefficients the same. From Figures 5-6, it is shown that the equilibrium solution E∗ is asymptotically stable when τ=2.5 and spatially inhomogeneous periodic solution exists when τ=5.42. It confirms that both spatial diffusion and toxin delay have significant effects on the spatiotemporal dynamics of system (3). Besides, if we set Ω=(0,2π), then the spatially inhomogeneous periodic solution vanishes and the equilibrium solution becomes asymptotically stable (see Figure 7). It can be concluded that spatially inhomogeneous pattern is more prone to occurring in small space.Figure 5
The equilibrium solutionE∗ of (3) is asymptotically stable when τ=2.5 with new coefficient values.
(a) (b) (c)Figure 6
The spatially inhomogeneous periodic solution exists whenτ=5.42 with new coefficient values.
(a) (b) (c)Figure 7
The equilibrium solutionE∗ of (3) is asymptotically stable when τ=5.42 and Ω=(0,2π).
(a) (b) (c)


## 7. Conclusions

In this paper, we have proposed a delayed reaction-diffusion model incorporating three plankton populations. Mathematical analysis indicates that the three populations can coexist when the natural growth rate of nontoxic phytoplankton is large and the death rate of zooplankton by toxin is small. In this case, the balance is finally achieved by interdependence and mutual restraint. Otherwise, the zooplankton will become extinct if the natural growth rate of nontoxic phytoplankton is small and the death rate of zooplankton by toxin is large. In the latter case, the biomasses of nontoxic phytoplankton and toxin-producing phytoplankton reach saturation and algae bloom occurs.And when considering the toxin delay, the spatiotemporal dynamics of system is almost unaffected when the time delay is sufficiently small. However, when the time delay passes through some critical value, spatially homogeneous or inhomogeneous periodic solution may arise. This means that algal bloom erupts periodically under certain conditions. Therefore, some measures can be adopted to control or defer the occurrence of Hopf bifurcation, such as reducing the toxin delay or adding feedback control. Numerical simulations complementally illustrate that when there is no or only tiny time delay, the population distribution pattern is eventually spatially homogeneous even if the initial distribution is inhomogeneous. When the time delay is sufficiently large, the distribution pattern is time periodic. Moreover, nontoxic phytoplankton has the largest oscillation amplitude, toxin-producing phytoplankton has the smallest oscillation amplitude, and the zooplankton has modest amplitude. These phenomena show that toxin delay has the greatest effects on nontoxic phytoplankton.



---

*Source: 4590915-2019-07-11.xml*