# Stability Switches and Hopf Bifurcation Induced by Nutrient Recycling Delay in a Reaction-Diffusion Nutrient-Phytoplankton Model

Authors: Kejun Zhuang; Ying Li; Bin Gong
Journal: Complexity
 (2021)

---

## Abstract

In this paper, we study a reaction-diffusion nutrient-phytoplankton model with nutrient recycling delay and toxin-producing phytoplankton. For the delay-free case, we prove that the release rate of toxic chemicals by the toxin-producing phytoplankton can destabilize the system and cause periodic oscillation. For the time-delayed case, we observe that nutrient recycling delay may bring about stability switches and Hopf bifurcation. Moreover, we derive the formula for determining the direction and stability of the bifurcating periodic solutions. Finally, we give some numerical simulations to support the theoretical analyses.


---

## Body

## 1. Introduction

Aquatic ecosystem has a significant effect on human’s living environment. It not only plays an important role in maintaining the global biological matter and water cycles but also undertakes the functions of water source, power source, transportation pollution purification place, and so on. Obviously, it is of great importance to investigate the dynamic behaviors of planktonic ecosystem. In the last few decades, lots of mathematical models have been constructed in terms of differential equations to describe the dynamical interactions between the dissolved nutrient and plankton population, see [1–5] and references cited therein. These models can deepen our knowledge of aquatic ecosystem cycle. Some interesting dynamical behaviors have been investigated, such as uniform persistence, asymptotic stability, Hopf bifurcation, Turing instability, and travelling waves.In nature, the dead phytoplankton can be partly transformed into the nutrient. By taking into account the recycling factor, Chakraborty et al. [6] proposed the following nutrient-phytoplankton model:(1)dNdt=a−bNP−dN+βP,dPdt=cNP−αP−θP2H2+P2,where N and P denote the concentrations of nutrient and phytoplankton population, respectively. All the coefficients are positive constants, a is the constant nutrient input flow, b is the nutrient uptake rate of phytoplankton, c is the nutrient conversion rate of phytoplankton, d is the per capita-loss rate of nutrient, α is the per capita-mortality rate of phytoplankton, β is the portion of the phytoplankton recycled back to the nutrient concentration, θ is the release rate of toxic chemicals by the toxin-producing phytoplankton population, and H is the half-saturation constant. In view of the conversion efficiency in ecological system, we always assume that α>β and b>c throughout the paper.For system (1), the limit cycle behavior is obtained by numerical simulations in [6]. It is found that inclusion of nutrient recycling enhances the chances of recurring bloom dynamics. After that, some expanded models were investigated by incorporating zooplankton population in [7–11].As indicated in [12], no life is possible if spatial and temporal processes do not interact. Actually, the individuals can go anywhere in the spatial domain since a well-known fact about Fickian diffusion leading to a Laplacian term. In the lakes or oceans, dissolved nutrient and phytoplankton population can move with the tide. The spatial diffusion has been considered in modelling the biological and ecological systems, such as [13–15]. Reaction-diffusion equations are widely used for modelling chemical reactions, biological systems, population dynamics, and nuclear reactor physics. This kind of models can better describe the spatial and temporal interactions. So, it is more realistic to use the reaction-diffusion equations to model the aquatic system. Then, we have the following model with spatial diffusion:(2)∂Nx,t∂t=d1ΔNx,t+a−bNx,tPx,t−dNx,t+βPx,t,x∈Ω,t>0,∂Px,t∂t=d2ΔPx,t+cNx,tPx,t−αPx,t−θP2x,tH2+P2x,t,x∈Ω,t>0,∂Nx,t∂ν=∂Px,t∂ν=0,x∈∂Ω,t>0,Nx,0=N0x≥0,Px,0=P0x≥0,x∈Ω¯,where d1 and d2 are the diffusion coefficients of nutrient and phytoplankton, respectively. The initial functions N0x and P0x are nonnegative and nonidentical to zero. The homogeneous Neumann boundary conditions mean that the nutrient and phytoplankton do not leave the living region. For simplicity, we specify the spatial location as an interval Ω=0,π, where π denotes the depth of the water column and can be changed into any depth through a suitable transformation. Moreover, some other nutrient-plankton models incorporating spatial effects are investigated in [16–21].In models (1) and (2), it is assumed that the decomposition process of dead phytoplankton is instantaneous. However, it is not accurate enough. Usually, the interaction between nutrient and phytoplankton is not instantaneous and it does take time to transform dead phytoplankton into nutrient. Time delay can also bring rich dynamical behaviors, especially the stability switches, and hence it cannot be ignored [22–24]. Therefore, we further build the following model with nutrient recycling delay motivated by [25, 26]:(3)∂Nx,t∂t=d1ΔNx,t+a−bNx,tPx,t−dNx,t+βPx,t−τ,x∈Ω,t>0,∂Px,t∂t=d2ΔPx,t+cNx,tPx,t−αPx,t−θP2x,tH2+P2x,t,x∈Ω,t>0,∂Nx,t∂ν=∂Px,t∂ν=0,x∈∂Ω,t>0,Nx,t=N1x,t≥0,Px,t=P1x,t≥0,x∈Ω¯,t∈−τ,0,where the nutrient recycling delay τ>0 is the decomposition time of dead phytoplankton population, Nx,t,Px,t∈CΩ¯×−τ,∞,ℝ+2. The initial functions N1x,t and P1x,t are nonnegative and nonidentical to zero. The aim of this paper is to investigate the joint effects of spatial diffusion and time delay on the spatiotemporal dynamics of the nutrient-phytoplankton system.The highlights of this paper mainly contain two aspects. First, we synthetically consider the nutrient recycling delay and spatial diffusion in modelling the interactions between nutrient concentration and phytoplankton population, and the model is more closer to the real world. Second, by regarding nutrient recycling delay as the bifurcation parameter, we find that time delay may induce stability switches and sustained oscillations, which can help to understand the cyclical burst of algae bloom.The rest of this paper is organized as follows. In Section2, we establish the existence of nonnegative constant equilibria. In Section 3, we analyze the stability of positive equilibrium in the absence of time delay, the stability switches, and Hopf bifurcation in the presence of time delay, respectively. In Section 4, the algorithm for determining the Hopf bifurcating periodic solutions is derived. In Section 5, some numerical simulations are presented to verify the theoretical results. Finally, a brief conclusion is given in Section 6.


## 2. Existence of Equilibria

In this section, we give some preliminary results. For convenience, we make the following assumptions:(H1) bθ>ac−αd>0.(H2) T0<0,D0>0.(H3) T0=0 and P∗−H≠0.(H4) d2bP∗+d−d1θP∗P∗2−H2/H2+P∗22>0.(H5) d1θP∗P∗2−H2/H2+P∗22−d2bP∗+d>2d1d2D0.(H6) Tn4−4AnTn2+4B2>0 and Tn2−2An<0 for some n∈ℕ0.Spatial diffusion and time delay do not affect the number and location of differential system, so systems (1)–(3) have the same nonnegative constant equilibria and the equilibria are the roots of the following algebraic equation:(4)a−bNP−dN+βP=0,cNP−αP−θP2H2+P2=0.Obviously, there is always the boundary equilibriumE0a/d,0. Next, we establish the existence of positive constant equilibrium. From the second equation of (4), we have(5)N=αc+θPcH2+P2.Substituting it into the first equation of (4) leads to(6)a+βP=d+bPαc+θPcH2+P2,that is,(7)a3P3+a2P2+a1P+a0=0,where a3=αb−βc>0, a2=αd+bθ−ac, a1=dθ+αbH2−βcH2>0, a0=H2αd−ac. Based on the Descartes’ rule of signs, equation (7) has one or three positive roots and assumption H1 can guarantee the existence of unique positive real root P∗ of equation (7).Thus, we can establish the existence of positive equilibrium.Lemma 1.
If (H1) holds, then systems (1)–(3) have the unique positive equilibrium E∗N∗,P∗.Remark 1.
From Lemma1, if the nutrient uptake rate of phytoplankton and release rate of toxic chemicals by toxin-producing phytoplankton population are appropriately large, and the net growth rates of nutrient and phytoplankton are appropriately small, then the unique positive equilibrium of system may exist.


## 3. Stability Analysis and Occurrence of Hopf Bifurcation

In the following, we analyze the stability ofE∗. Linearizing system (1) at E∗, we have(8)∂Nx,t∂t∂Px,t∂t=DΔNx,tΔPx,t+L1Nx,tPx,t+L2Nx,t−τPx,t−τ,where(9)D=d100d2,L1=−bP∗−d−bN∗cP∗θP∗P∗2−H2H2+P∗22,L2=0β00.Then, the characteristic equation of (8) is(10)detλI−Mn−L1−L2e−λτ=0,where I=diag1,1, Mn=−n2diag1,1, n∈ℕ0. Equation (10) can be turned into(11)λ2−Tnλ+An+Be−λτ=0,where(12)Tn=−d1+d2n2+θP∗P∗2−H2H2+P∗22−bP∗−d,An=d1d2n4+d2bP∗+d−d1θP∗P∗2−H2H2+P∗22n2+bcN∗P∗−θP∗P∗2−H2bP∗+dH2+P∗22,B=−βcN∗.


### 3.1. The Case ofτ=0

Ifτ=0 and n=0, then characteristic equation (11) can be simplified to(13)λ2−T0λ+D0=0,which corresponds to ordinary differential system (1) and D0=A0+B. By Routh–Hurwitz criterion, the positive equilibrium E∗ of (1) is asymptotically stable if the condition (H2) is satisfied.Next, we explore the influence of parameterθ on the stability of E∗. Let T0=0, then θ∗=bP∗+dH2+P∗22/P∗P∗2−H2 is the potential Hopf bifurcation value. In addition, we have(14)dReλθdθ|θ=θ∗=P∗P∗2−H22H2+P∗22.From above, we can conclude the stability results of system (1).Theorem 1.
(i)
If (H1) and (H2) are satisfied, then the positive equilibrium E∗ of (1) is asymptotically stable.(ii)
If (H1) and (H3) are satisfied, then the positive equilibrium E∗ of (1) is asymptotically stable when θ<θ∗ and unstable when θ>θ∗ for P∗>H. Conversely, E∗ of (1) is asymptotically stable when θ>θ∗ and unstable when θ<θ∗ for P∗<H. Moreover, θ∗ is the Hopf bifurcation value.Remark 2.
From Theorem1, if the equilibrium state of phytoplankton P∗ is smaller than the half-saturation constant, then condition (H2) can easily be satisfied and the nutrient and phytoplankton can achieve a balanced coexistence in the absence of diffusion. Otherwise, the release rate of toxic chemicals may induce the periodic variations of nutrient and phytoplankton.Here, we mainly investigate the Turing instability induced by spatial diffusion. The linearization of system (2) at the positive equilibrium E∗ is(15)∂N∂t=d1ΔN−bP∗+dN+β−bN∗P,∂P∂t=d2ΔP+cP∗N+θP∗P∗2−H2H2+P∗22P,and the corresponding characteristic equation is(16)λ2−Tnλ+Dn=0,where(17)Tn=−d1+d2n2+θP∗P∗2−H2H2+P∗22−bP∗−d,Dn=d1d2n4+d2bP∗+d−d1θP∗P∗2−H2H2+P∗22n2+cP∗bN∗−β−θP∗P∗2−H2bP∗+dH2+P∗22.It is obvious thatTn<0 when T0<0. Then, the positive equilibrium E∗ of diffusive system (2) is stable if Dn>0 for any n and unstable if Dn<0 for some n>0. It is not difficult to determine the sign of Dn by using the properties of parabola.Following Turing’s idea in [27], different diffusion rates of two chemical reactants could sometimes destabilize an otherwise stable equilibrium of the reaction-diffusion system and lead to nonuniform spatial patterns. This kind of diffusion-driven instability is usually called Turing instability.Theorem 2.
Suppose (H1) and (H2) hold, we have the following statements:(i)
If (H4) is satisfied, then the positive equilibrium E∗ of (2) is locally asymptotically stable.(ii)
If (H5) is satisfied, then the positive equilibrium E∗ of (2) is unstable. It implies that Turing instability occurs.Remark 3.
From Theorem2, the diffusion coefficients have significant effect on the stability of positive equilibrium in the presence of spatial diffusion. If the diffusion rate of phytoplankton is larger, then the aquatic system is stable. On the contrary, the system may be unstable.


### 3.2. The Case ofτ>0

Forτ>0, let λ=iωω>0 be a root of (11). We have(18)−ω2−iωTn+An+Bcosωτ−isinωτ=0.Separating the real and imaginary parts can lead to(19)Bcosωτ=ω2−An,Bsinωτ=−ωTn.Squaring both sides of the above equations, we get(20)ω4+Tn2−2Anω2+An2−B2=0.FromDn=An+B and B<0, it can be found that Dn>0 implies An−B>0. Moreover, equation (20) has two positive real roots under the assumption (H6). Thus, characteristic equation (11) has two pairs of purely imaginary roots λ=±iωn±, where(21)ωn±=2An−Tn2±Tn4−4AnTn2+4B22,n∈ℕ0,τnj±=1ωn±2π−arccosωn±2−AnB+2jπ,j=0,1,2,….By direct computation, we haveωn−<ωn+ and ωn−2−An/B>ωn+2−An/B>0. Thus, τnj−>τnj+.In what follows, we shall verify the transversality condition.Lemma 2.
Redλ/dττ=τnj+−1>0 and Redλ/dττ=τnj−−1<0.Proof.
Differentiating both sides of equation (11) with respect to τ yields(22)dλdτ−1=2λ−TnBλe−λτ−τλ=2λ−TneλτBλ−τλ.
Then,(23)Redλdττ=τnj±−1=Re2iωn±−Tncosωn±τnj±+isinωn±τnj±iBωn±=2ωn±cosωn±τnj±−Tnsinωn±τnj±Bωn±=2ωn±ωn±2−An+Tn2ωn±B2ωn±=2ωn±2−An+Tn2B2=±Tn4−4AnTn2+4B2B2.Through the above analysis and motivated by Corollary 2.4 of Ruan and Wei in [28], we can establish the distribution of roots of characteristic equation (11) in the case of n=0.Lemma 3.
If (H1), (H2), (H4), and (H6) are satisfied, then we have the following:(i)
Ifτ01+<τ00−, then all roots of (11) have negative real parts for τ∈0,τ00+; equation (11) has a pair of purely imaginary roots and other roots have negative real parts for τ=τ00+; equation (11) has at least a root with positive real part for τ>τ00+.(ii)
Ifτ01+>τ00−, then there exists a positive integer k such that all roots of (11) have negative real parts for τ∈0,τ00+∪τ00−,τ01+∪⋯∪τ0,k−1−,τ0k+, and equation (11) has at least one root with positive real part for τ∈τ00+,τ00−∪τ01+,τ01−∪⋯∪τ0,k−1+,τ0,k−1−∪τ0k+,+∞.As time delay changes, the stability of the equilibrium solution may also change. The stability may transit from stable to unstable at some critical value and may also transit from unstable to stable at some other critical value. Such phenomena are often referred to as stability switches, i.e., the state can be switched from stable to unstable and back to stable and so on just by progressive increase of the time delay. Then, we can obtain the stability switches and existence of Hopf bifurcation forE∗ in (3).Theorem 3.
If (H1), (H2), (H4), and (H6) are satisfied, then the following conclusions are true:(i)
Ifτ01+<τ00−, then positive equilibrium E∗ of (3) is asymptotically stable for τ∈0,τ00+ and is unstable for τ>τ00+. System (3) undergoes Hopf bifurcations near E∗ when τ=τ00+.(ii)
Ifτ01+>τ00−, then there exists a positive integer k such that E∗ switches k times from stability to instability to stability and so on such that E∗ is asymptotically stable for τ∈0,τ00+∪τ00−,τ01+∪⋯∪τ0,k−1−,τ0k+ and is unstable for τ∈τ00+,τ00−∪τ01+,τ01−∪⋯∪τ0,k−1+,τ0,k−1−∪τ0k+,+∞. System (3) undergoes Hopf bifurcation around E∗ for every τ=τ0j±.


## 3.1. The Case ofτ=0

Ifτ=0 and n=0, then characteristic equation (11) can be simplified to(13)λ2−T0λ+D0=0,which corresponds to ordinary differential system (1) and D0=A0+B. By Routh–Hurwitz criterion, the positive equilibrium E∗ of (1) is asymptotically stable if the condition (H2) is satisfied.Next, we explore the influence of parameterθ on the stability of E∗. Let T0=0, then θ∗=bP∗+dH2+P∗22/P∗P∗2−H2 is the potential Hopf bifurcation value. In addition, we have(14)dReλθdθ|θ=θ∗=P∗P∗2−H22H2+P∗22.From above, we can conclude the stability results of system (1).Theorem 1.
(i)
If (H1) and (H2) are satisfied, then the positive equilibrium E∗ of (1) is asymptotically stable.(ii)
If (H1) and (H3) are satisfied, then the positive equilibrium E∗ of (1) is asymptotically stable when θ<θ∗ and unstable when θ>θ∗ for P∗>H. Conversely, E∗ of (1) is asymptotically stable when θ>θ∗ and unstable when θ<θ∗ for P∗<H. Moreover, θ∗ is the Hopf bifurcation value.Remark 2.
From Theorem1, if the equilibrium state of phytoplankton P∗ is smaller than the half-saturation constant, then condition (H2) can easily be satisfied and the nutrient and phytoplankton can achieve a balanced coexistence in the absence of diffusion. Otherwise, the release rate of toxic chemicals may induce the periodic variations of nutrient and phytoplankton.Here, we mainly investigate the Turing instability induced by spatial diffusion. The linearization of system (2) at the positive equilibrium E∗ is(15)∂N∂t=d1ΔN−bP∗+dN+β−bN∗P,∂P∂t=d2ΔP+cP∗N+θP∗P∗2−H2H2+P∗22P,and the corresponding characteristic equation is(16)λ2−Tnλ+Dn=0,where(17)Tn=−d1+d2n2+θP∗P∗2−H2H2+P∗22−bP∗−d,Dn=d1d2n4+d2bP∗+d−d1θP∗P∗2−H2H2+P∗22n2+cP∗bN∗−β−θP∗P∗2−H2bP∗+dH2+P∗22.It is obvious thatTn<0 when T0<0. Then, the positive equilibrium E∗ of diffusive system (2) is stable if Dn>0 for any n and unstable if Dn<0 for some n>0. It is not difficult to determine the sign of Dn by using the properties of parabola.Following Turing’s idea in [27], different diffusion rates of two chemical reactants could sometimes destabilize an otherwise stable equilibrium of the reaction-diffusion system and lead to nonuniform spatial patterns. This kind of diffusion-driven instability is usually called Turing instability.Theorem 2.
Suppose (H1) and (H2) hold, we have the following statements:(i)
If (H4) is satisfied, then the positive equilibrium E∗ of (2) is locally asymptotically stable.(ii)
If (H5) is satisfied, then the positive equilibrium E∗ of (2) is unstable. It implies that Turing instability occurs.Remark 3.
From Theorem2, the diffusion coefficients have significant effect on the stability of positive equilibrium in the presence of spatial diffusion. If the diffusion rate of phytoplankton is larger, then the aquatic system is stable. On the contrary, the system may be unstable.


## 3.2. The Case ofτ>0

Forτ>0, let λ=iωω>0 be a root of (11). We have(18)−ω2−iωTn+An+Bcosωτ−isinωτ=0.Separating the real and imaginary parts can lead to(19)Bcosωτ=ω2−An,Bsinωτ=−ωTn.Squaring both sides of the above equations, we get(20)ω4+Tn2−2Anω2+An2−B2=0.FromDn=An+B and B<0, it can be found that Dn>0 implies An−B>0. Moreover, equation (20) has two positive real roots under the assumption (H6). Thus, characteristic equation (11) has two pairs of purely imaginary roots λ=±iωn±, where(21)ωn±=2An−Tn2±Tn4−4AnTn2+4B22,n∈ℕ0,τnj±=1ωn±2π−arccosωn±2−AnB+2jπ,j=0,1,2,….By direct computation, we haveωn−<ωn+ and ωn−2−An/B>ωn+2−An/B>0. Thus, τnj−>τnj+.In what follows, we shall verify the transversality condition.Lemma 2.
Redλ/dττ=τnj+−1>0 and Redλ/dττ=τnj−−1<0.Proof.
Differentiating both sides of equation (11) with respect to τ yields(22)dλdτ−1=2λ−TnBλe−λτ−τλ=2λ−TneλτBλ−τλ.
Then,(23)Redλdττ=τnj±−1=Re2iωn±−Tncosωn±τnj±+isinωn±τnj±iBωn±=2ωn±cosωn±τnj±−Tnsinωn±τnj±Bωn±=2ωn±ωn±2−An+Tn2ωn±B2ωn±=2ωn±2−An+Tn2B2=±Tn4−4AnTn2+4B2B2.Through the above analysis and motivated by Corollary 2.4 of Ruan and Wei in [28], we can establish the distribution of roots of characteristic equation (11) in the case of n=0.Lemma 3.
If (H1), (H2), (H4), and (H6) are satisfied, then we have the following:(i)
Ifτ01+<τ00−, then all roots of (11) have negative real parts for τ∈0,τ00+; equation (11) has a pair of purely imaginary roots and other roots have negative real parts for τ=τ00+; equation (11) has at least a root with positive real part for τ>τ00+.(ii)
Ifτ01+>τ00−, then there exists a positive integer k such that all roots of (11) have negative real parts for τ∈0,τ00+∪τ00−,τ01+∪⋯∪τ0,k−1−,τ0k+, and equation (11) has at least one root with positive real part for τ∈τ00+,τ00−∪τ01+,τ01−∪⋯∪τ0,k−1+,τ0,k−1−∪τ0k+,+∞.As time delay changes, the stability of the equilibrium solution may also change. The stability may transit from stable to unstable at some critical value and may also transit from unstable to stable at some other critical value. Such phenomena are often referred to as stability switches, i.e., the state can be switched from stable to unstable and back to stable and so on just by progressive increase of the time delay. Then, we can obtain the stability switches and existence of Hopf bifurcation forE∗ in (3).Theorem 3.
If (H1), (H2), (H4), and (H6) are satisfied, then the following conclusions are true:(i)
Ifτ01+<τ00−, then positive equilibrium E∗ of (3) is asymptotically stable for τ∈0,τ00+ and is unstable for τ>τ00+. System (3) undergoes Hopf bifurcations near E∗ when τ=τ00+.(ii)
Ifτ01+>τ00−, then there exists a positive integer k such that E∗ switches k times from stability to instability to stability and so on such that E∗ is asymptotically stable for τ∈0,τ00+∪τ00−,τ01+∪⋯∪τ0,k−1−,τ0k+ and is unstable for τ∈τ00+,τ00−∪τ01+,τ01−∪⋯∪τ0,k−1+,τ0,k−1−∪τ0k+,+∞. System (3) undergoes Hopf bifurcation around E∗ for every τ=τ0j±.


## 4. Properties of Hopf Bifurcation

Here, we derive the Hopf bifurcation properties by using the center manifold theorem and normal form theory for partial functional differential equations in [29, 30]. For simplicity, we denote the Hopf bifurcation value by τ∗ and the corresponding simply purely imaginary characteristic roots by ±iω∗, respectively.LetX=N,P∈W2,2Ω:∂Nx,t/∂x=∂Px,t/∂x=0,x=0,π, N˜⋅,t=N⋅,τt, P˜⋅,t=P⋅,τt, U˜t=N˜⋅,t,P˜⋅,t, and τ=τ∗+μ with μ∈ℝ, then μ=0 is the Hopf bifurcation value of system (3). For simplicity, we drop the tilde and rewrite system (3) in the form(24)dUtdt=τDΔUt+LμUt+fUt,μ,where D=diagd1,d2, φ=φ1,φ2T∈C=C−1,0,X, and Lμ⋅:C⟶X, f:C×ℝ⟶X are given by(25)Lμφ=τ∗+μ−bP∗−dφ10−bN∗φ20+βφ2−1cP∗φ10+θP∗P∗2−H2H2+P∗22φ20,fφ,μ=τ∗+μ∑i+j+l≥21i!j!l!Fijlφ1i0φ2j0φ2l−1∑i+j≥21i!j!Gijφ1i0φ2j0,and(26)Fu,v,w=a−buv−du+βw,Gu,v=cuv−αv−θv2H2+v2,Fijl=∂i+j+lF∂ui∂vj∂wlN∗,P∗,P∗,Gij=∂i+jG∂ui∂vjN∗,P∗.From the previous section, we know that±iω∗τ∗ is a pair of simple purely imaginary eigenvalues of the following linear differential equation:(27)U˙t=τDΔUt+LμUt.Due to the Riesz representation theorem, there exists a2×2 matrix function ηϑ,μ−1≤ϑ≤0, whose elements are of bounded variation functions such that(28)−τ∗Dn2φ0+Lμφ=∫−10dηϑ,τ∗φθforφ∈C−1,0,ℝ2,where(29)ηϑ,τ∗=τ∗−bP∗−d−d1n2−bN∗cP∗θP∗P∗2−H2H2+P∗22−d2n2,ϑ=0,0,ϑ∈−1,0,−τ∗0β00,ϑ=−1.ForΦ∈C1−1,0,ℝ2, Ψ∈C10,1,ℝ2, we define(30)A1Φϑ=dΦϑdϑ,ϑ∈−1,0,∫−10dηϑ,τ∗Φϑ,ϑ=0,A1∗Ψs=−dΨsds,s∈0,1,∫−10dηTϑ,τ∗Ψ−ϑ,s=0.Then,A1∗ and A1 are adjoint operators under the bilinear form(31)Ψ,Φ0=Ψ¯0Φ0−∫−10∫0ϑΨ¯ζ−ϑdηϑ,0Φζdζ=Ψ¯0Φ0+τ∗∫−10Ψ¯ζ+100r−rΦζdζ.It is calculated thatqϑ=1,q1Teiω∗τ∗ϑϑ∈−1,0 and q∗s=Mq2,1Teiω∗τ∗ss∈0,1 are eigenvectors of A1 and A1∗ corresponding to iω∗τ∗ and −iω∗τ∗, respectively, where(32)q1=iω∗+d1n2+bP∗+dbN∗−βe−iω∗τ∗,q2=cP∗−iω∗+d1n2+bP∗+d,M=η+ξ¯+βηξ¯τ∗e−iω∗τ∗−1.Then,Q=spanqϑ,qϑ¯ is the center subspace of system (3), and Q∗=spanq∗s,q∗s¯ is the adjoint subspace.Leth⋅fn=h1βn1+h2βn2, fn=βn1,βn2, and βn1=cosnx,0T, βn2=0,cosnxT. The complex-valued L2 inner product on Hilbert space XC is(33)U1,U2=1π∫0πu1v1¯+u2v2¯dx,for U1=u1,u2,U2=v1,v2∈XC, and β0i,β0i=1, βni,βni=1/2, i=1,2, n=1,2,…,(34)Φ,fn=Φ,βn1,Φ,βn2,where Φ∈C−1,0,X. Then, the center subspace of system (27) at μ=0 is(35)PCNℒ=qϑz+qϑ¯z¯⋅fn,z∈ℂ.Ifμ=0, then we can obtain(36)Wz,z¯,θ=W20z22+W11zz¯+W02z¯22+⋯.The flow of system (24) on the center manifold can be written as(37)Ut=qθzt+qθ¯z¯t⋅fn+Wzt,z¯t,θ.Moreover, forUt∈C0 of (24) at τ=τ∗, we have z˙=iω∗τ∗z+gz,z¯, where(38)gz,z¯=q∗0¯fUt,0,fn=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.Following the calculation procedures in [29, 30], we can get(39)g20=τ∗M¯−2bξη¯+2cξ+G22ξ2,g11=τ∗M¯−2bReξη¯+2cReξ+G22ξ2,g02=g20¯,g21=−M¯τ∗π∫0πbη¯cos2nx⋅W2020+ξ¯W2010+2ξW1110+2W1120dx+M¯τ∗π∫0πccos2nx⋅W2020+ξ¯W2010+2ξW1110+2W1120+G22cos2nx⋅ξ¯W2020+2ξW1120dx,where(40)W20ϑ=ig20ω∗τ∗qϑ+ig02¯3ω∗τ∗qϑ¯⋅fn+E1e2iω∗τ∗ϑ,W11ϑ=−ig11ω∗τ∗qϑ+ig11¯qϑ¯ω∗τ∗⋅fn+E2,E1=E1′×−2bξη¯2cξ+G22ξ2cos2nx,E1′=2iω∗+d1n2+bP∗+dbN∗−βe−2iω∗τ∗−cP∗2iω∗+d2n2−θP∗P∗2−H2H2+P∗22−1,E2=E2′×−2bReξ2cReξ+G22ξ2cos2nx,E2′=d1n2+bP∗+dbN∗−β−cP∗d2n2−θP∗P∗2−H2H2+P∗22−1.Thus, we can compute(41)c10=i2ω∗τ∗g20g11−2g112−13g022+g212,ℓ2=−Rec10Reλ′τ∗,ι2=2Rec10,〉2=−1ω∗τ∗Imc10+ℓ2Imλ′τ∗.Then, by the Hopf bifurcation theory in [29, 30], we have the following conclusions about the properties of the Hopf bifurcation.Theorem 4.
For system (3),(i)
ℓ2 determines the direction of Hopf bifurcation: if ℓ2>0, then the supercritical bifurcation bifurcates stable limit cycles; if ℓ2<0, then the subcritical bifurcation bifurcates unstable limit cycles.(ii)
ι2 determines the stability of bifurcating periodic solutions: the periodic solutions are orbitally asymptotically stable if ι2<0 and are unstable if ι2>0.(iii)
ϱ2 determines the period of the bifurcating periodic solutions: the period is monotonically increasing at the time delay τ when ϱ2>0 and is monotonically decreasing at the time delay τ when ϱ2<0.


## 5. Numerical Simulations

In this section, we mainly verify the theoretical results obtained in the previous sections by some numerical examples. For system (3), we set Ω=0,π, d1=d2=1, a=0.03, b=1.5, c=1, d=0.003, k=0.05, H=0.035, α=0.1, and β=0.05.Example 1.
Ifθ=0.04 and τ=0 with the initial conditions N1x,t=0.5+0.2sinx+t and P1x,t=0.05+0.05sinx+t, then the positive constant equilibrium E∗0.6636,0.0296 of system (2) is asymptotically stable, as shown in Figure 1.Figure 1
The positive equilibriumE∗ of (2) is asymptotically stable when θ=0.02.Example 2.
Ifθ=0.02 and τ=0 with the initial conditions N1x,t=0.5+0.2sinx+t and P1x,t=0.05+0.02sinx+t, then the positive constant equilibrium E∗0.6636,0.0296 of system (2) is unstable and Hopf bifurcation occurs, see Figure 2.
Forτ>0, we fix θ=0.029. Then, system (3) has the unique positive equilibrium solution E∗0.5110,0.0397. Moreover, equation (20) has two positive roots ω0+≈0.2294 and ω0−=0.0408, τ00+≈14<τ00−≈41<τ01+≈69<τ01−≈96<⋯. According to Theorem 3, if time delay τ varies, then the stability of positive equilibrium E∗0.5110,0.0397 of system (3) would change.Figure 2
The positive equilibriumE∗ of (2) is unstable when θ=0.02.Example 3.
Forτ00+≈14, we can compute c10≈−2.3543−1.2902i, then ℓ2>0 and ι2<0. From Theorem 4, we can know that the Hopf bifurcation is supercritical and the spatially homogeneous periodic solutions are stable. The numerical solutions with initial conditions N1x,t=0.51 and P1x,t=0.035 are presented in Figures 3 and 4.Figure 3
The positive equilibriumE∗ of (3) is stable when τ=10.Figure 4
The positive equilibriumE∗ of (3) is unstable when τ=30.Example 4.
Forτ01+≈69, we can compute c10≈−3.8745−1.6482i, then ℓ2>0 and ι2<0. The Hopf bifurcation is supercritical and the spatially homogeneous periodic solutions are stable. The numerical solutions with initial conditions N1x,t=0.51 and P1x,t=0.035 are presented in Figures 5 and 6.
It can be found that nutrient recycling delay can bring about the stability switches in the nutrient-phytoplankton system.Figure 5
The positive equilibriumE∗ of (3) is stable when τ=50.Figure 6
The positive equilibriumE∗ of (3) is unstable when τ=75.


## 6. Conclusion

In this paper, we have expanded the model proposed in [6] and investigated a novel reaction-diffusion nutrient-phytoplankton ecological model incorporating nutrient recycling delay. We proved that the release rate of toxic chemicals by the toxin-producing phytoplankton population θ can destabilize the positive constant equilibrium and may generate time-periodic solutions. We also revealed the effect of nutrient recycling delay on our model. By taking time delay as the bifurcation parameter, multiple stability switches of the positive constant equilibrium and Hopf bifurcation phenomena can be observed under some conditions, which contrast sharply with the case in the absence of time delay. The results can help to understand the periodic outbreak of algae.It is observed that the model is more comprehensive by taking account of the zooplankton population, and rich dynamical behaviors have been presented in [3, 8, 10, 26], such as uniform persistence, Turing bifurcation, and global Hopf bifurcation. However, these models do not investigate the joint effect of spatial diffusion and time delay. In addition, in a reaction-diffusion system, when the Hopf bifurcation curve and Turing bifurcation curve intersect, the Turing-Hopf bifurcation will occur and generate spatiotemporal patterns, see [31–33]. We will do further research on spatiotemporal patters via Turing-Hopf bifurcation in the delayed reaction-diffusion model incorporating three components in the future following the aforementioned references.



---

*Source: 7943788-2021-10-20.xml*