Stability switches and chaos induced by delay in a reaction-diffusion nutrient-plankton model

Abstract

In this paper, we investigate a reaction-diffusion model incorporating dynamic variables for nutrient, phytoplankton, and zooplankton. Moreover, we account for the impact of time delay in the growth of phytoplankton following nutrient uptake. Our theoretical analysis reveals that the time delay can trigger the emergence of persistent oscillations in the model via a Hopf bifurcation. We also analytically track the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions. Our simulation results demonstrate stability switches occurring for the positive equilibrium with an increasing time lag. Furthermore, the model exhibits homogeneous periodic-2 and 3 solutions, as well as chaotic behaviour. These findings highlight that the presence of time delay in the phytoplankton growth can bring forth dynamical complexity to the nutrient-plankton system of aquatic habitats.

Keywords:

• Nutrient-plankton model
• time delay
• Hopf bifurcation
• stability switches
• chaos

MATHEMATICS SUBJECT CLASSIFICATIONS:

• 65E05
• 37M05
• 37M10
• 37N25
• 39A05

1. Introduction

Algal blooms occur in aquatic systems due to the massive growth of phytoplankton, leading to major water-quality problems and significant negative impacts on human health [1–3]. Algal blooms are the consequences of the interplay of various hydrodynamics, chemical, and biological processes [4]. Nutrient levels and zooplankton predation have been identified as key factors influencing the occurrence of algal blooms [5,6]. Vanni [7] reported through experimental observations that even small changes in zooplankton size can have a tremendous impact on the phytoplankton community, while nutrient levels play a more significant role in the growth of phytoplankton population. Despite considerable research on phytoplankton growth dynamics, the specific process of phytoplankton bloom formation remains unclear. Previous studies have emphasized the need for in-depth investigations into nutrient-plankton interaction.

In recent years, extensive experimental studies and field observations have been conducted to study plankton dynamics. However, due to the complexity and nonlinearity of aquatic ecosystems, understanding the mechanisms behind phytoplankton growth solely through experimental and/or field observations remains challenging [8]. Mathematical models serve as a powerful tool, providing quantitative insights into the complex dynamics of plankton communities, which have been alternatively employed by researchers [9–11]. Notably, all the aforementioned researches have been conducted with an assumption that the growth process of phytoplankton is instantaneous. However, experimental evidence indicates that under high nutrient availability, a time lag occurs in the growth response of Isochrysis galbana [12]. That is, the interaction between nutrient and phytoplankton will not be essentially instantaneous, suggesting that there is a delay in the growth of phytoplankton. Therefore, it is reasonable to consider the effect of time delay in the phytoplankton growth. Both experimental and field observations have demonstrated that changes in plankton population density usually exhibit oscillatory behaviour due to the existence of internal factors [13,14], supporting that the population density of phytoplankton in aquatic reservoirs is not constant but changes over time in realistic scenarios [15]. In recent years, the effect of time delay has been widely studied in the dynamics of phytoplankton growth [16–20]. Growing evidences suggest that plankton models with time delay exhibit rich and complex dynamical behaviours. Previous studies have demonstrated that the delay can destabilize a system via Hopf bifurcation and induce periodic solutions [21–23], create stability switches [24,25], and even lead to chaotic dynamics [26,27].

Moreover, it is well established that nutrients and plankton can disperse in aquatic reservoirs due to currents and turbulent diffusion. Actually, both nutrient and plankton exhibit constant movement, leading to population densities that fluctuate spatially as well as temporally. Growing evidences suggest that plankton growth and distribution can be characterized by its spatial variation. Building on [1,28], a variety of ecological phenomena, including fundamental mechanisms of complex spatio-temporal plankton dynamics, can be modelled using appropriate partial differential equations. Segel and Jackson [29] were the first to introduce spatial structure in population dynamics. Subsequently, numerous plankton models have been scrutinized by considering the effect of diffusion [30–34]. Chakraborty et al. [35] have investigated the spatial dynamics of a nutrient-phytoplankton model with the toxicity of phytoplankton. Their findings showed that in the presence of toxicity, the distributions of nutrient and phytoplankton in the aquatic ecosystem become inhomogeneous in space, and produce different patterns, such as stripes, spots, and mixtures, depending on the level of toxicity. They have also observed that the distributions of nutrient and phytoplankton show spatio-temporal oscillations for certain toxicity levels.

In recent years, numerous studies have been conducted on delay-induced dynamical behaviours in diffusive phytoplankton growth models [36–38]. The findings of Dai et al. [38] showed that the presence of delay not only triggers instability in a positive equilibrium but also promotes the emergence of patchiness through Hopf bifurcation. Moreover, the recent developments in dynamical model theory have considered delay-induced chaotic fluctuations of a dynamical model as highly desirable because fluctuations allow for easier control of such a system. Additionally, while various mechanisms, such as competition [39], predation [40], and resource and consumer digested delay [41], contribute to chaotic dynamics, the phenomenon of delay-induced chaos in spatio-temporal plankton growth models has drawn considerable scientific attention. Gökce et al. [42] demonstrated that delay is capable of switching the model from stable to unstable through Hopf bifurcation, resulting in spatio-temporal chaos. Despite the extensive discussion about the characteristics of delay-induced chaos and its presence in phytoplankton growth models, evidence for its existence in spatio-temporal phytoplankton growth models remains scarce.

It is worth noting that despite extensive investigations into the spatio-temporal dynamics of nutrient-plankton models, the influence of time delay involved in the growth of phytoplankton after the nutrient uptake on the spatial distributions of nutrient, phytoplankton, and zooplankton in aquatic ecosystems remains unclear. Therefore, in the present study, we aim to explore the role of such time delay in a diffusive model for the nutrient-plankton interactions in an aquatic reservoir, considering the release of toxic chemicals by certain phytoplankton species to mitigate predation pressure from zooplankton.

The remaining portion of this paper is organized as follows. In Section 2, we propose our model for the combined actions of delay and diffusion in a nutrient-plankton system of an aquatic ecosystem. In Section 3, we conduct mathematical analysis of the model's dynamics. Firstly, we obtain the equilibrium points of the temporal model in the absence of time delay. Then, we discuss the stability behaviour of the model in the presence of delay and diffusion. We also demonstrate the existence of Hopf bifurcation analytically and analyse the direction and stability of the bifurcating periodic solutions. In Section 4, we perform numerical simulations to illustrate the analytical findings and gain further insights into the dynamics of the delayed-diffusive nutrient-plankton model. Finally, in Section 5, we present the conclusions of this study.

2. The model formulation

One approach to the study of the dynamics of ecological community is via its food web and the coupling of interacting species with each other [43]. We develop a mathematical model to explore the impact of time lag concealed in the growth of phytoplankton after the uptake of nutrient in a diffusive nutrient-plankton system of an aquatic ecosystem. The model encompasses three dynamic variables: nutrient concentration, phytoplankton biomass, and zooplankton biomass. The growth dynamics of the phytoplankton and zooplankton populations mainly depend on the availability of nutrient and phytoplankton, respectively. The aquatic system's nutrient concentration is determined by input and washout rates, as well as uptake by the phytoplankton population. For the phytoplankton population, its biomass is determined by nutrient uptake, grazing by zooplankton, and the mortality, either natural or due to nutrient scarcity. The predation of phytoplankton by zooplankton is described by a more appropriate functional form [44]. Finally, the zooplankton biomass in the aquatic reservoir is influenced by the grazing of phytoplankton, toxin liberation by phytoplankton, and mortality.

Let $N(t,x)$, $P(t,x)$, and $Z(t,x)$ denote the concentration of nutrient, biomass of phytoplankton, and the biomass of zooplankton, respectively, at any given time t>0 and position x. To build our model, we make the following ecological assumptions.

1. There is a constant input of nutrient into the aquatic system; the nutrients depleting naturally at a rate b and being taken up by phytoplankton through bilinear interaction.

2. The predation rate of zooplankton on phytoplankton follows the Holling type II functional form [44], $\frac{cPZ}{h+P}$; with c is denoting the maximum uptake rate of phytoplankton by zooplankton and h is representing the half-saturation constant.

3. The parameters m and k express the natural mortality rates of phytoplankton and zooplankton species, respectively. The phytoplankton species also experience mortality due to competition for available nutrient in the aquatic system; the parameter r signifies the strength of intraspecific competition among phytoplankton species owing to limited resources.

4. The biomass of phytoplankton is converted into the biomass of zooplankton at a rate d. Toxin-producing phytoplankton species release toxic chemicals at a rate ρ, leading to a decrement in the biomass of zooplankton at a rate $\frac{\rho PZ}{h+P}$.

5. The growth response of phytoplankton following nutrient uptake is mediated by some time lag.

6. In the aquatic system, nutrient and plankton dispersal occur due to currents and turbulent diffusion.

The above assumptions yield the following delayed reaction-diffusion nutrient-plankton model: (1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) In the above model, the parameter τ represents the time which is required by phytoplankton to absorb nutrient and reproduce. The constants $d_1,d_2$, and $d_3$ denote the self-diffusion coefficients of nutrient, phytoplankton and zooplankton, respectively. Let $x\in \mathrm{\Omega }=[0,l\pi ]$, and for $t\in [-\tau ,0]$, $N(x,t)={\varphi }_1(x,t)\ge 0$, $P(x,t)={\varphi }_2(x,t)\ge 0$ and $Z(x,t)={\varphi }_3(x,t)\ge 0$. All the parameters in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) are assumed to be positive, and their biological meanings are listed in Table 1.

Table 1. Biological meanings of parameters involved in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ).

Parameters	Descriptions	Values
α	Nutrient input rate to the aquatic system	$2^{\left(\mathrm{a}\right)}$
b	Natural washout rate of nutrient	${0.4}^{\left(\mathrm{a}\right)}$
e	Uptake rate of nutrient by phytoplankton	${0.1}^{\left(\mathrm{b}\right)}$
β	Conversion rate of nutrient concentration into the biomass of phytoplankton	${0.9677}^{\left(\mathrm{c}\right)}$
c	Uptake rate of phytoplankton by zooplankton	${0.1}^{\left(\mathrm{b}\right)}$
d	Conversion rate of phytoplankton biomass into the biomass of zooplankton	${0.9661}^{\left(\mathrm{c}\right)}$
h	Half-saturation constant	${3.5}^{\left(\mathrm{b}\right)}$
m	Natural mortality rate of phytoplankton	$0.01-{0.65}^{\left(\mathrm{d}\right)}$
k	Natural mortality rate of zooplankton	${0.2}^{\left(\mathrm{a}\right)}$
r	Strength of intraspecific competition among the phytoplankton population	Estimated
ρ	Rate of toxin liberation by toxin-producing phytoplankton	$0.22/0.3/{0.92}^{\left(\mathrm{e}\right)}$

Note: (a) [45]; (b) [25]; (c) [46]; (d) [47]; (e) [48].

3. Analytical results of model (1)

At first, we discuss about the existence of ecologically meaningful equilibria of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) by ignoring the delay and diffusion factors. Then, we perform the linear stability and bifurcation analysis by taking the delay factor as a bifurcation parameter.

3.1. The existence of equilibria

Model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) does not have any trivial solution, but it has a non-trivial solution $E^0=(N^0,P^0,Z^0)$, where $P^0=\frac{\text{β}{N}^0-m}{r}$, $Z^0=0$ and $N^0$ is the positive solution of the following quadratic equation: $e\text{β}{N}^2+(br-em)N-\alpha r=0.$The only one positive root of above equation is given by $N^0=\frac{(em-br)+\sqrt{(br-em{)}^2+4e\text{β}\alpha r}}{2e\text{β}}.$Additionally, if the following conditions are satisfied: $\alpha >\frac{\left\{m\right(d-k-\rho )+hkr\}\left\{b\right(d-k-\rho )+ehk\}}{\text{β}(d-k-\rho {)}^2},d>k+\rho ,$then model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) admits a unique positive equilibrium $E^{\ast }=(N^{\ast },P^{\ast },Z^{\ast })$ whose components are given by $N^{\ast }=\frac{\alpha }{b+e{P}^{\ast }},P^{\ast }=\frac{hk}{d-k-\rho },Z^{\ast }=\frac{(\text{β}{N}^{\ast }-m-r{P}^{\ast })(h+P^{\ast })}{c}.$From an ecological point of view, the positive equilibrium $E^{\ast }$ is very important as all the considered dynamical variables are present here. So, we will focus only on the dynamics of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) around this equilibrium point in the forthcoming sections.

3.2. Stability and bifurcation

Here, we will discuss about the stability behaviour of the positive equilibrium $E^{\ast }$ and the existence of Hopf bifurcation in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ). As we have considered zero-flux boundary conditions, let $0={\mu }_0<{\mu }_1<{\mu }_2<\cdots ({\mu }_j=j^2/l^2$, $j\in {\mathbb{N}}_0=\{0,1,2,\dots \})$ be the eigenvalues of the operator $-\mathrm{\Delta }$ on Ω with the homogeneous Neumann boundary condition. Following [38], we represent model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) as an abstract functional differential equation in the abstract space $C\left(\right[-\tau ,0],W)$. The linear form of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) around the equilibrium $E^{\ast }$ is given by (2) $\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}=D\mathrm{\Delta }U\left(t\right)+\mathrm{\Phi }\left(U_t\right),$(2) where $D=\text{diag}(d_1,d_2,d_3)$, and $\begin{array}{rl}\text{dom}\left(D\mathrm{\Delta }\right)& =\left\{\right(N,P,Z{)}^T:N,P,Z\in C^3\left(\right[0,l\pi ],R),\frac{\mathrm{\partial }N}{\mathrm{\partial }x}=\frac{\mathrm{\partial }P}{\mathrm{\partial }x}=\frac{\mathrm{\partial }Z}{\mathrm{\partial }x}=0,x=0,l\pi \},\\ \mathrm{\Phi }\left(\phi \right)& =\left(\begin{array}{c}-\left[\right(b+e{P}^{\ast }\left){\phi }_1\right(0)+(e{N}^{\ast }\left){\phi }_2\right(0\left)\right]\\ \text{β}{P}^{\ast }{\phi }_1(-\tau )+(\frac{c{P}^{\ast }{Z}^{\ast }}{(h+P^{\ast }{)}^2}-r{P}^{\ast }){\phi }_2\left(0\right)-\frac{c{P}^{\ast }}{h+P^{\ast }}{\phi }_3\left(0\right)\\ \frac{(d-\rho )h{Z}^{\ast }}{(h+P^{\ast }{)}^2}{\phi }_2\left(0\right)\end{array}\right).\end{array}$Following [49], we write the characteristic equation corresponding to Equation (2(2) $\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}=D\mathrm{\Delta }U\left(t\right)+\mathrm{\Phi }\left(U_t\right),$(2) ) as follows: (3) $\lambda y-D\mathrm{\Delta }y-\mathrm{\Phi }\left({\mathrm{e}}^{\lambda }y\right)=0,y\in \text{dom}\left(D\mathrm{\Delta }\right),y\ne 0.$(3) Now, we define (4) $y=\sum _{j=0}^{+\mathrm{\infty }}{\mathrm{\Psi }}_j\left(x\right)\left(\begin{array}{c}{y}_{1,j}\\ y_{2,j}\\ y_{3,j}\end{array}\right),$(4) where ${\mathrm{\Psi }}_j\left(x\right)$ comprises of the eigenfunctions corresponding to the eigenvalues ${\mu }_j$. On substituting Equations (4(4) $y=\sum _{j=0}^{+\mathrm{\infty }}{\mathrm{\Psi }}_j\left(x\right)\left(\begin{array}{c}{y}_{1,j}\\ y_{2,j}\\ y_{3,j}\end{array}\right),$(4) ) into (3(3) $\lambda y-D\mathrm{\Delta }y-\mathrm{\Phi }\left({\mathrm{e}}^{\lambda }y\right)=0,y\in \text{dom}\left(D\mathrm{\Delta }\right),y\ne 0.$(3) ), we get $\begin{array}{r}\left(\begin{array}{ccc}-(b+e{P}^{\ast }+d_1{\mu }_j)& -\left(e{N}^{\ast }\right)& 0\\ \text{β}{P}^{\ast }{\mathrm{e}}^{-\lambda \tau }& \frac{c{P}^{\ast }{Z}^{\ast }}{(h+P^{\ast }{)}^2}-r{P}^{\ast }-d_2{\mu }_j& -\frac{c{P}^{\ast }}{h+P^{\ast }}\\ 0& \frac{(d-\rho )h{Z}^{\ast }}{(h+P^{\ast }{)}^2}& -d_3{\mu }_j\end{array}\right)\left(\begin{array}{c}{y}_{1,j}\\ y_{2,j}\\ y_{3,j}\end{array}\right)=\lambda \left(\begin{array}{c}{y}_{1,j}\\ y_{2,j}\\ y_{3,j},\end{array}\right)\end{array}$for $j=0,1,2,\dots .$

Thus, an equivalent equation of Equation (3(3) $\lambda y-D\mathrm{\Delta }y-\mathrm{\Phi }\left({\mathrm{e}}^{\lambda }y\right)=0,y\in \text{dom}\left(D\mathrm{\Delta }\right),y\ne 0.$(3) ) is obtained as follows: (5) ${\lambda }^3+B_j{\lambda }^2+C_j\lambda +D_j{\mathrm{e}}^{-\lambda \tau }+E\lambda {\mathrm{e}}^{-\lambda \tau }+F_j=0,j=0,1,2,\dots ,$(5) where $\begin{array}{rl}{B}_j& =-(a_{11}+a_{22})+(d_1+d_2+d_3)\frac{j^2}{l^2},\\ C_j& =a_{11}{a}_{22}-a_{23}{a}_{32}-(a_{11}{d}_2+a_{22}{d}_1+a_{11}{d}_3+a_{22}{d}_3)\frac{j^2}{l^2}+(d_1{d}_3+d_2{d}_3+d_1{d}_2)\frac{j^4}{l^4},\\ D_j& =-a_{21}{a}_{12}{d}_3\frac{j^2}{l^2},E=-a_{12}{a}_{21},\\ F_j& =a_{11}{a}_{23}{a}_{32}+(a_{11}{a}_{22}{d}_3-a_{23}{a}_{32}{d}_1)\frac{j^2}{l^2}-(a_{11}{d}_2{d}_3+a_{22}{d}_1{d}_3)\frac{j^4}{l^4}+d_1{d}_2{d}_3\frac{j^6}{l^6}\end{array}$with $\begin{array}{rl}{a}_{11}& =-b-e{P}^{\ast },a_{12}=-e{N}^{\ast },a_{32}=\frac{(d-\rho )h{Z}^{\ast }}{(h+P^{\ast }{)}^2},a_{21}=\text{β}{P}^{\ast },a_{22}=\frac{c{P}^{\ast }{Z}^{\ast }}{(h+P^{\ast }{)}^2},\\ a_{23}& =-\frac{c{P}^{\ast }}{h+P^{\ast }}.\end{array}$(H1) Let $\frac{a_{11}}{a_{22}}{d}_2<d_1<\left(\frac{a_{11}{a}_{22}+a_{12}{a}_{21}}{a_{23}{a}_{32}}\right)d_3$. Then, we have the following result.

Theorem 3.1

If $b>max\{d-p,\frac{\alpha \text{β}}{m}\}$ and (H1) holds, then the positive equilibrium $E^{\ast }$ of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) without delay is locally asymptotically stable for any $j\in {\mathbb{N}}_0$.

Proof.

For $\tau =0$, Equation (5(5) ${\lambda }^3+B_j{\lambda }^2+C_j\lambda +D_j{\mathrm{e}}^{-\lambda \tau }+E\lambda {\mathrm{e}}^{-\lambda \tau }+F_j=0,j=0,1,2,\dots ,$(5) ) reduces to the following form: (6) ${\lambda }^3+B_j{\lambda }^2+(C_j+E)\lambda +D_j+F_j=0.$(6) Since $a_{21}{a}_{12}<0$, we have $D_j+F_j>0$ in view of (H1). Obviously, $\lambda =0$ is not a solution of Equation (5(5) ${\lambda }^3+B_j{\lambda }^2+C_j\lambda +D_j{\mathrm{e}}^{-\lambda \tau }+E\lambda {\mathrm{e}}^{-\lambda \tau }+F_j=0,j=0,1,2,\dots ,$(5) ). If $b>max\{d-p,\frac{\alpha \text{β}}{m}\}$, then $B_j>0$ and $C_j+E>0$ for any $j\in {\mathbb{N}}_0$. Thus, Equation (6(6) ${\lambda }^3+B_j{\lambda }^2+(C_j+E)\lambda +D_j+F_j=0.$(6) ) has no positive roots, and hence Theorem 3.1.

Now, we consider the nutrient-plankton model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) with $\tau >0$. Following [50], we assume that $i\omega (\omega >0)$ is a root of Equation (5(5) ${\lambda }^3+B_j{\lambda }^2+C_j\lambda +D_j{\mathrm{e}}^{-\lambda \tau }+E\lambda {\mathrm{e}}^{-\lambda \tau }+F_j=0,j=0,1,2,\dots ,$(5) ). Putting $\lambda =i\omega$ into Equation (5(5) ${\lambda }^3+B_j{\lambda }^2+C_j\lambda +D_j{\mathrm{e}}^{-\lambda \tau }+E\lambda {\mathrm{e}}^{-\lambda \tau }+F_j=0,j=0,1,2,\dots ,$(5) ), and separating real and imaginary parts, we get the following two equations: (7) $\begin{array}{rl}-{\omega }^3+C_j\omega & =D_j\sin \left(\omega \tau \right)-E\omega \cos \left(\omega \tau \right),\\ -B_j{\omega }^2+F_j& =-D_j\cos \left(\omega \tau \right)-E\omega \sin \left(\omega \tau \right).\end{array}$(7) Squaring and adding both sides of the above equation, and simplifying, we get (8) ${\omega }^6+M_{1,j}{\omega }^4+M_{2,j}{\omega }^2+Z_j=0,j=0,1,2,\dots$(8) with $M_{1,j}=B_j^2-2C_j$, $M_{2,j}=C_j^2-2B_j{F}_j-E^2$, $Z_j=F_j^2-D_j^2$. Now, we present the following lemma.

Lemma 3.2

Let $b>max\{d-p,\frac{\alpha \text{β}}{m}\}$ and (H1) holds. We have the following two cases.

(i)	If ${\mathrm{\Delta }}_1=M_{1,j}^2-3M_{2,j}\le 0$ for any $j\in {\mathbb{N}}_0$, then Equation (8(8) ${\omega }^6+M_{1,j}{\omega }^4+M_{2,j}{\omega }^2+Z_j=0,j=0,1,2,\dots$(8) ) does not have any positive root.
(ii)	Let a $j_0\in {\mathbb{N}}_0$ be such that ${\mathrm{\Delta }}_1=M_{1,j_0}^2-3M_{2,j_0}>0$, then Equation (8(8) ${\omega }^6+M_{1,j}{\omega }^4+M_{2,j}{\omega }^2+Z_j=0,j=0,1,2,\dots$(8) ) possesses two positive roots if and only if $x_{j_0,1}>0$ and $f\left(x_{j_0,1}\right)<0$.

Proof.

Let ${\omega }^2=v$, then from Equation (8(8) ${\omega }^6+M_{1,j}{\omega }^4+M_{2,j}{\omega }^2+Z_j=0,j=0,1,2,\dots$(8) ), we get (9) $f\left(v\right)=v^3+M_{1,j}{v}^2+M_{2,j}v+Z_j.$(9) Clearly, $Z_j>0$ in view of (H1). We denote by ${△}_1=M_{1,j}^2-3M_{2,j}$. Evidently, for ${△}_1<0$, Equation (9(9) $f\left(v\right)=v^3+M_{1,j}{v}^2+M_{2,j}v+Z_j.$(9) ) has no positive root. However, if there exists a $j_0\in {\mathbb{N}}_0$ such that ${△}_1>0$, then equation $f^{\mathrm{\prime }}\left(v\right)=0$ has the following two real roots: $x_{j_0,1}=\frac{-M_{1,j_0}+\sqrt{{\mathrm{\Delta }}_1}}{3},x_{j_0,2}=\frac{-M_{1,j_0}-\sqrt{{\mathrm{\Delta }}_1}}{3}.$For $Z_j$, ${\mathrm{\Delta }}_1>0$, Equation (9(9) $f\left(v\right)=v^3+M_{1,j}{v}^2+M_{2,j}v+Z_j.$(9) ) has two positive real roots if and only if $x_{j_0,1}>0$ and $f\left(x_{j_0,1}\right)<0$.

Now, we present the following lemma about the transversality condition that is required for the presence of Hopf bifurcation in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ).

Lemma 3.3

If the second condition of Lemma 3.2 holds, then $\frac{\mathrm{d}Re\lambda \left(\tau \right)}{\mathrm{d}\tau }{\mid }_{\tau ={\tau }_{j_0,1}^s}>0,\frac{\mathrm{d}Re\lambda \left(\tau \right)}{\mathrm{d}\tau }{\mid }_{\tau ={\tau }_{j_0,2}^s}<0,s=0,1,2,\dots .$

Proof.

Denoting the two positive roots of Equation (9(9) $f\left(v\right)=v^3+M_{1,j}{v}^2+M_{2,j}v+Z_j.$(9) ) by ${\chi }_{j_0,1}$ and ${\chi }_{j_0,2}$, then the two positive roots of Equation (8(8) ${\omega }^6+M_{1,j}{\omega }^4+M_{2,j}{\omega }^2+Z_j=0,j=0,1,2,\dots$(8) ) will be ${\omega }_{j_0,1}=\sqrt{{\chi }_{j_0,1}}$ and ${\omega }_{j_0,2}=\sqrt{{\chi }_{j_0,2}}$. From Equation (7(7) $\begin{array}{rl}-{\omega }^3+C_j\omega & =D_j\sin \left(\omega \tau \right)-E\omega \cos \left(\omega \tau \right),\\ -B_j{\omega }^2+F_j& =-D_j\cos \left(\omega \tau \right)-E\omega \sin \left(\omega \tau \right).\end{array}$(7) ), we obtain the critical value of time delay as, (10) $\begin{array}{rl}{\tau }_{j_0,p}^s& =\frac{1}{{\omega }_{j_0,p}}({\cos }^{-1}\frac{E{\omega }_{j_0,p}^4+(B_{j_0}{D}_{j_0}-C_{j_0}E){\omega }_{j_0,p}^2-D_{j_0}{F}_{j_0}}{D_{j_0}^2+E^2{\omega }_{j_0,p}^2}+2s\pi ),\\ & (s=0,1,2,\dots ,p=1,2).\end{array}$(10) Obviously, the right side of Equation (10(10) $\begin{array}{rl}{\tau }_{j_0,p}^s& =\frac{1}{{\omega }_{j_0,p}}({\cos }^{-1}\frac{E{\omega }_{j_0,p}^4+(B_{j_0}{D}_{j_0}-C_{j_0}E){\omega }_{j_0,p}^2-D_{j_0}{F}_{j_0}}{D_{j_0}^2+E^2{\omega }_{j_0,p}^2}+2s\pi ),\\ & (s=0,1,2,\dots ,p=1,2).\end{array}$(10) ) is an increasing function of s. Thus, we have ${\tau }_c=\underset{p=1,2}{min}{\tau }_{j_0,p}^0$. Let $\lambda \left(\tau \right)=\mu \left(\tau \right)+i\omega \left(\tau \right)$ be a root of Equation (5(5) ${\lambda }^3+B_j{\lambda }^2+C_j\lambda +D_j{\mathrm{e}}^{-\lambda \tau }+E\lambda {\mathrm{e}}^{-\lambda \tau }+F_j=0,j=0,1,2,\dots ,$(5) ) such that (11) $\mu \left({\tau }_{j_0,p}^s\right)=0,\omega \left({\tau }_{j_0,p}^s\right)={\omega }_k,(s=0,1,2,\dots ,p=1,2).$(11) On differentiating both sides of Equation (5(5) ${\lambda }^3+B_j{\lambda }^2+C_j\lambda +D_j{\mathrm{e}}^{-\lambda \tau }+E\lambda {\mathrm{e}}^{-\lambda \tau }+F_j=0,j=0,1,2,\dots ,$(5) ) with respect to τ, we obtain (12) ${\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right)}^{-1}=\frac{3{\lambda }^2+2B_j\lambda +C_j+E{e}^{-\lambda \tau }}{\lambda (D_j{\mathrm{e}}^{-\lambda \tau }+E\lambda ){\mathrm{e}}^{-\lambda \tau }}-\frac{\tau }{\lambda }.$(12) This yields, (13) $\frac{\mathrm{d}Re\left(\lambda \right({\tau }_{j_0,p}^s\left)\right)}{\mathrm{d}\tau }=\frac{{\omega }_k^2}{{\mathrm{\Delta }}_2}{f}^{\mathrm{\prime }}\left({\omega }_k^2\right),$(13) where ${\mathrm{\Delta }}_2=\left(D_{j_0}{\omega }_{k}sin\right({\omega }_k{\tau }_{j_0,p}^s)-E{\omega }_k^{2}cos({\omega }_k{\tau }_{j_0,p}^s){)}^2+(D_{j_0}{\omega }_{k}cos\left({\omega }_k{\tau }_{j_0,p}^s\right)+E{\omega }_k^{2}sin\left({\omega }_k{\tau }_{j_0,p}^s\right){)}^2.$For convenience, we assume that ${\tau }_{j_0,1}^s<{\tau }_{j_0,2}^s$. Notably, the second point of Lemma 3.2 yields that $f^{\mathrm{\prime }}\left({\chi }_{j_0,1}\right)>0$ and $f^{\mathrm{\prime }}\left({\chi }_{j_0,2}\right)<0$. Hence, $\frac{\mathrm{d}Re\lambda \left(\tau \right)}{\mathrm{d}\tau }{\mid }_{\tau ={\tau }_{j_0,1}^s}>0$ and $\frac{\mathrm{d}Re\lambda \left(\tau \right)}{\mathrm{d}\tau }{\mid }_{\tau ={\tau }_{j_0,2}^s}<0$ for $s=0,1,2,\dots$.

By Lemma 3.3, Hopf bifurcation occurs in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) at the critical value $\tau ={\tau }_c$, i.e. the positive equilibrium $E^{\ast }$ becomes unstable for $\tau >{\tau }_c$. We summarize these results in the following theorem.

Theorem 3.4

For model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ), we find that

(i)

If Lemma 3.2 (i) holds, then the positive equilibrium $E^{\ast }$ of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) is locally asymptotically stable for all $\tau \ge 0$.

(ii)

If Lemma 3.2 (ii) holds, then there exists a non-negative integer n such that the positive equilibrium $E^{\ast }$ of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) is locally asymptotically stable when $\tau \in [0,{\tau }_{j_0,1}^0)\cup ({\tau }_{j_0,2}^0,{\tau }_{j_0,1}^1)\cup \cdots \cup ({\tau }_{j_0,2}^{n-1},{\tau }_{j_0,1}^n)$, and is unstable when $\tau \in ({\tau }_{j_0,1}^0,{\tau }_{j_0,2}^0)\cup ({\tau }_{j_0,1}^1,{\tau }_{j_0,2}^1)\cup \cdots \cup ({\tau }_{j_0,1}^{n-1},{\tau }_{j_0,2}^{n-1})\cup ({\tau }_{j_0,1}^n,+\mathrm{\infty })$. Further, model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) undergoes a Hopf bifurcation around the equilibrium point $E^{\ast }$ at every $\tau ={\tau }_{j_0,p}^s$ for p = 1, 2 and $s=0,1,2,\dots$.

Now, we state the following results which are based on the centre manifold theorem and the normal form theory.

Theorem 3.5

The Hopf bifurcation in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) is characterized by the signs of ${\mu }_2,{\text{β}}_2$ and $T_2$ as follows.

(i)	If ${\mu }_2>0(<0)$, Hopf bifurcation is supercritical (subcritical) at $\tau ={\tau }^{\ast }$.
(ii)	If ${\text{β}}_2<0(>0)$, then the bifurcating periodic solutions are stable (unstable).
(iii)	If $T_2>0(<0)$, then period of the bifurcating periodic solutions increases (decreases).

The proof of this theorem is placed in Appendix.

4. Numerical simulations

In this section, we numerically investigate the influence of time delay on the dynamics of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ). For the numerical investigations of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ), we adopt the following set of parameter values that is in accordance with the values mentioned in Table 1: $\begin{array}{rl}b& =0.4,e=0.09,\text{β}=0.95,c=0.09,h=3,m=0.01,r=0.2,\\ d& =0.95,k=0.2,d_1=d_2=0.01,d_3=\mathrm{0.02.}\end{array}$Unless specified otherwise, these parameter values are used consistently throughout the simulation section. Here, we choose $\alpha ,\rho$, and τ as control parameters in the nutrient-plankton model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ). Further, we choose the initial values for nutrient, phytoplankton, and zooplankton as 0.9, 0.9, and 22, respectively.

At first, we examine the dynamics of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) in the absence of time delay ($\tau =0$). Figure 1 illustrates the temporal and spatial variations in the biomass of phytoplankton and zooplankton populations. As depicted in the figure, the corresponding equilibrium point is asymptotically stable, implying that the plankton populations remain constant over time. Recall that the theoretical analysis indicates that Hopf bifurcation may occur when there is a time lag in the growth of phytoplankton after nutrient uptake. Moreover, stability switches may occur with gradual increments in the value of time delay. Numerically, we found that the natural mortality of phytoplankton species does not affect the model's dynamics much. Thus, in order to highlight the dynamical behaviours of stability switches, we choose m = 0.19, and present the number of stability switches experienced by the model by varying the parameters ρ and α (see Figure 2(a)). In the figure, region I represents the unstable domain of the positive equilibrium, while region X indicates stability if the positive equilibrium exists. In region II, the positive equilibrium loses its stability when τ exceeds the critical value ${\tau }_0$. Further, region III exhibits the emergence of stability switches. To provide a closer examination of the stability switches, an enlarged view of Figure 2(a) is presented in Figure 2(b).

Figure 1. Biomass distributions of (a) phytoplankton and (b) zooplankton populations in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) over time and space for $\tau =0$.

Figure 2. (a) The number of stability switches of the positive equilibrium of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) in the $\rho -\alpha$ plane. Stability switches occur once in region III, twice in region IV, three times in region V, four times in region VI, five times in region VII, six times in region VIII and seven times in region IX. (b) Enlarged diagram of regions V to X in (a).

Following Ref. [38], we perform a sensitive analysis that is based on one-dimensional space solutions, Figure 3. The figures depict that model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) exhibits sensitivity to variations in the parameters α and ρ. Figure 3(a) illustrates that when the nutrient input rate is low (region I), model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) does not exhibit any positive equilibrium point. However, as the value of α increases, model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) experiences stability switches. Further, Figure 3(b) displays that when the toxin liberation rate is low, the biomass of phytoplankton oscillates. On further increasing the value of ρ, the dynamics of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) change from unstable to stable via Hopf bifurcation. However, once ρ enters region I, the positive equilibrium of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) does not exist. The results from Figure 3 suggest that the nutrient input rate has the potential to generate as well as suppress the unstable coexistence of plankton populations, while a higher toxin liberation rate contributes to the stability of the system. However, very high toxin liberation rate can pose a threat to the entire ecosystem.

Figure 3. Sensitive analysis of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) based on one-dimensional space solutions. The solid blue line represents the stable equilibrium, and the solid red line represents the maximal and minimal amplitudes of phytoplankton biomass. Figures show the effects of (a) nutrient input rate (α) on the model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) with $\rho =0.05$ and $\tau =15$, and (b) rate of toxin liberation (ρ) on the model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) with $\alpha =0.35$ and $\tau =15$.

Further, we draw bifurcation diagrams of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) in the ρ-τ plane for $\alpha =0.45$ (see Figure 4(a)) and $\alpha =0.35$ (see Figure 4(b)). These figures illustrate the existence of stability switches in the system before ρ enters the green zone (region I). Within the green zone, the positive equilibrium remains stable regardless of the value of time delay. However, if ρ falls within region II, the positive equilibrium does not exist for any value of τ. As our aim is to examine the effect of time lag on the dynamics of the nutrient-plankton system, we present the bifurcation diagram of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) with respect to the delay parameter τ, focusing solely on the biomass of the phytoplankton population (see Figure 5). In Figure 5(a), we find that for $\alpha =0.45$, the stability switch occurs once and periodic-2 solutions emerge as the value of time delay increases. The biomass distributions of phytoplankton over time and space at $\tau =2$, 10, and 18 are depicted in Figure 6, where the positive equilibrium is stable at $\tau =2$ and 18 but unstable at $\tau =10$. Figure 5(b) shows that for $\alpha =0.35$ stability switches occur twice as the value of time delay gradually increases, indicating the emergence of transitions from stable to unstable via Hopf bifurcation. The biomass distributions at $\tau =$3, 18, and 42 (stable) and $\tau =$10, 36, and 50 (unstable) are plotted in Figure 7. Overall, our findings indicate that the nonlinear analysis is in agreement with the results predicted by the linear analysis (see Figures 4 and 5). The results indicate that delay in the growth of phytoplankton is the only reason for the rise in the instability around the positive equilibrium. Moreover, Figures 6 and 7 suggest that delay can generate as well as suppress oscillations around the coexistence equilibrium point.

Figure 4. Bifurcation diagrams of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) with respect to ρ and τ for (a) $\alpha =0.45$ and (b) $\alpha =0.35$. In the figures, the solid, dashed, dash-dot, and dotted curves represent the critical values of τ for j = 0, 1, 2 and 3, respectively.

Figure 5. Bifurcation diagrams of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) with respect to τ for (a) $\alpha =0.45$ and $\rho =0.1$, and (b) $\alpha =0.35$ and $\rho =0.01$. In the figure, the solid blue line and the dotted red line represent the stable positive equilibrium and the maximum biomass of phytoplankton, respectively.

Figure 6. Biomass distribution of phytoplankton population in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) over time and space for $\alpha =0.45$, $\rho =0.1$, and different values of τ: (a) $\tau =2$, (b) $\tau =10$ and (c) $\tau =18$.

Figure 7. Biomass distribution of phytoplankton population in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) over time and space for $\alpha =0.35$, $\rho =0.05$, and different values of τ: (a) $\tau =3$, (b) $\tau =10$, (c) $\tau =18$, (d) $\tau =36$, (e) $\tau =42$ and (f) $\tau =50$.

Finally, we present a bifurcation diagram for a wide range of the delay parameter τ, focusing on the biomass of the phytoplankton population (see Figure 8(a)). The figure clearly shows that model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) exhibits periodic-1, 2, and 3 solutions as the value of time delay gradually increases. It can be noted in the figure that for a very large value of time delay, the model showcases chaotic dynamics. The periodic-1, 2, and 3 solutions at $\tau =10$, 30, and 56 are marked with green, blue, and yellow solid diamonds, respectively. A vertical black dashed line is drawn at $\tau =74$, indicating the entrance into the chaotic regime. The biomass distributions of the phytoplankton population over time and space at these critical values of time delay are presented in Figure 8(b–e).

Figure 8. (a) Bifurcation diagram of model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) with respect to τ for $\alpha =0.6$ and $\rho =0.1$. In the figure, the green, blue and yellow solid diamonds respectively represent the periodic−1, 2 and 3 solutions at $\tau =10$, 30 and 56; the black dashed line denotes chaotic solution at $\tau =74$. Biomass distributions of phytoplankton in model (1(1) $\begin{array}{rl}\frac{\mathrm{\partial }N}{dt}& =\alpha -bN-eNP+d_1\mathrm{\Delta }N,\\ \frac{\mathrm{\partial }P}{dt}& =\text{β}N(t-\tau )P-\frac{cPZ}{h+P}-mP-r{P}^2+d_2\mathrm{\Delta }P,\\ \frac{\mathrm{\partial }Z}{dt}& =\frac{dPZ}{h+P}-kZ-\frac{\rho PZ}{h+P}+d_3\mathrm{\Delta }Z.\end{array}$(1) ) over time and space at (b) $\tau =10$, (c) $\tau =30$, (d) $\tau =56$ and (e) $\tau =74$.

5. Conclusions

The phenomenon of planktonic blooms has gathered significant attention from both ecologists and mathematicians in the recent decades [2]. Mathematical models have been widely used to investigate this phenomenon. Huppert et al. [51] investigated a simple mathematical model and found that high nutrient concentrations can trigger planktonic blooms in the aquatic reservoirs. Further, Steele and Henderson [52] explored the role of predation in the plankton models. Additionally, the effect of time delay is inevitable in the aquatic environments. Researchers have also recognized the significance of nutrient and plankton diffusion on the spatial distribution of phytoplankton in the aquatic ecosystems [53]. However, the dynamics induced by delay, such as oscillation, chaotic behaviour, and stability switches, in the phytoplankton growth models with diffusion require further exploration.

Here, we have investigated a nutrient-plankton model to examine the effects of time delay on the growth of phytoplankton and the interplay between nutrient and plankton populations. We also considered the effect of diffusion in the model to incorporate the spatial movements of nutrient and plankton populations in the aquatic ecosystem. Our theoretical results for the model without delay and diffusion showed that the positive equilibrium is locally asymptotically stable under some specific conditions on the parameters. In such a situation, the nutrient and plankton populations persist in the ecosystem. However, when the time delay involved in the growth of phytoplankton exceeds critical values, the positive equilibrium loses its stability, leading to the emergence of periodic solutions and oscillatory dynamics in the nutrient-plankton system. This indicates that the biomass of plankton populations cannot maintain a steady level. That is to say, time delay plays a significant role in inducing oscillations in the nutrient-plankton model. We also discussed the direction of Hopf bifurcation around the positive equilibrium and the stability of the bifurcating periodic solutions by employing the centre manifold theorem and the normal form theory. While many studies have demonstrated unstable coexistence of species in the ecosystems when the time delay exceeds a critical value [46,54,55], our simulation results reveal that time delay can generate as well as suppress oscillations in the system. This phenomenon is known as stability switches in the plankton system [41]. Obviously, these findings contribute to the control of plankton biomass in the aquatic reservoirs when the model exhibits stable coexistence of nutrient and plankton populations.

Furthermore, we explored the effects of time delay on the model's dynamics through bifurcation diagrams. Figures 4 and 5 illustrate the number of stability switches that occurred around the positive equilibrium. Unlike other nutrient-plankton models that only demonstrate a destabilizing role of time delay, our model exhibits both destabilizing as well as stabilizing effects of delay on the system's dynamics. Specifically, our model shows the occurrence of multiple stability switches around the positive equilibrium as the delay involved in the growth of phytoplankton increases. Figures 6 and 7 display the evolution of phytoplankton population concerning different time delay values. The results suggest that the time delay possesses the capability to generate as well as suppress oscillations in the plankton biomass, indicating the crucial role of time delay plays in the stability of nutrient-plankton system.

We observed that the system generates periodic-2 and 3 solutions, and ultimately becomes chaotic if there is a large time delay in the growth of phytoplankton following nutrient uptake. Thus, delay can be considered as a primary factor contributing to the chaotic behaviour of the plankton system. It is worth noting that the emergence of chaos may imply the unpredictability of biomass distributions of plankton populations over time and space. Therefore, the investigated nutrient-plankton model cannot accurately predict bloom events, but it contributes to enhancing our understanding of the influences of time delay and diffusion on the interplay between nutrient and plankton populations in the aquatic reservoirs. Our research provides important empirical evidence and specific insights into delay-induced oscillations within the context of the investigated nutrient-plankton model. Additionally, our findings offer valuable insights into the effect of delay on the emergence of chaos in the nutrient-plankton system.

Acknowledgments

All the authors thank the associate editor and anonymous reviewers for their valuable comments, which contributed to the improvement in the presentation of the paper. Qing Guo: Conceptualization, Methodology, Software, Writing – Original Draft. Lijun Wang: Writing – Review and Editing. He Liu: Writing – Review and Editing. Yi Wang: Writing – Review and Editing. Jianbing Li: Supervision, Resources, Writing – Review and Editing. Pankaj Kumar Tiwari: Writing – Review and Editing. Min Zhao: Writing – Review and Editing, Supervision, Resources, Funding acquisition. Chuanjun Dai: Software, Resources, Writing – Review and Editing.

Disclosure statement

The author declares that there are no conflicts of interest.

Additional information

Funding

This work was supported by the National Key Research and Development Program of China [grant number 2018YFE0103700], the National Natural Science Foundation of China [grant numbers 61901303, 61871293].

Appendix. The proof of Theorem 3.5

We will track the direction of Hopf bifurcation around the positive equilibrium and the stability of the bifurcating periodic solutions by using the normal formal theory and the centre manifold theorem of differential equations [56]. Let $U=(u,v,w),U_t=(u_t,v_t,w_t)$ and $\tau ={\tau }^{\ast }+\mu ,\mu \in R$. In space $\mathcal{C}$, Equation (2(2) $\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}=D\mathrm{\Delta }U\left(t\right)+\mathrm{\Phi }\left(U_t\right),$(2) ) can be rewritten as (A1) $\frac{\mathrm{d}U}{\mathrm{d}t}=\tilde{D}\mathrm{\Delta }U+L_{\mu }\left(U_t\right)+f(\mu ,U_t),$(A1) where $\tilde{D}=({\tau }^{\ast }+\mu )D$, $L_{\mu }\left(\phi \right)=({\tau }^{\ast }+\mu )\left(B_1\phi \right(0)+B_2\phi (-1\left)\right)$ and $\begin{array}{rl}f(\mu ,\phi )& =({\tau }^{\ast }+\mu )\\ & \times \left(\begin{array}{c}-e{\phi }_1\left(0\right){\phi }_2\left(0\right)+\cdots \\ \text{β}{\phi }_1(-1){\phi }_2\left(0\right)+(\frac{ch{Z}^{\ast }}{(h+P^{\ast }{)}^3}-r){{\phi }_2}^2\left(0\right)-ch(h+P^{\ast }{)}^{-2}{\phi }_2(0\left){\phi }_3\right(0)+\cdots \\ (\rho -d)h{Z}^{\ast }(h+P^{\ast }{)}^{-3}{{\phi }_2}^2(0)+(d-\rho \left)h\right(h+P^{\ast }{)}^{-2}{\phi }_2\left(0\right){\phi }_3\left(0\right)+\cdots \end{array}\right)\end{array}$with $B_1=\left(\begin{array}{ccc}{a}_{11}& a_{12}& 0\\ 0& a_{22}& a_{23}\\ 0& a_{32}& 0\end{array}\right),B_2=\left(\begin{array}{ccc}0& 0& 0\\ a_{21}& 0& 0\\ 0& 0& 0\end{array}\right)$for $\phi =({\phi }_1,{\phi }_2,{\phi }_3{)}^T\in \mathcal{C}$. Hence, $\mu =0$ is the bifurcation point of Equation (A1(A1) $\frac{\mathrm{d}U}{\mathrm{d}t}=\tilde{D}\mathrm{\Delta }U+L_{\mu }\left(U_t\right)+f(\mu ,U_t),$(A1) ).

On linearizing Equation (A1(A1) $\frac{\mathrm{d}U}{\mathrm{d}t}=\tilde{D}\mathrm{\Delta }U+L_{\mu }\left(U_t\right)+f(\mu ,U_t),$(A1) ) about $(0,0)$, we get (A2) $\frac{\mathrm{d}U}{\mathrm{d}t}=\tilde{D}\mathrm{\Delta }U+L_{\mu }\left(U_t\right).$(A2) Note that Equation (5(5) ${\lambda }^3+B_j{\lambda }^2+C_j\lambda +D_j{\mathrm{e}}^{-\lambda \tau }+E\lambda {\mathrm{e}}^{-\lambda \tau }+F_j=0,j=0,1,2,\dots ,$(5) ) has a pair of purely imaginary roots ${\mathrm{\Lambda }}_j=\{i{\omega }_j{\tau }^{\ast },-i{\omega }_j{\tau }^{\ast }\}$ when $\mu =0$ and other roots have negative real parts. The solution operator of Equation (A2(A2) $\frac{\mathrm{d}U}{\mathrm{d}t}=\tilde{D}\mathrm{\Delta }U+L_{\mu }\left(U_t\right).$(A2) ) is a semigroup whose infinitesimal generator is given by $A_{\mu }\phi =\{\begin{array}{ll}\dot{\phi }\left(\theta \right),& \theta \in [-r,0),\\ \tilde{D}\mathrm{\Delta }\phi \left(0\right)+L_{\mu }\left(\theta \right),& \theta \in 0.\end{array}$Hence, Equation (A2(A2) $\frac{\mathrm{d}U}{\mathrm{d}t}=\tilde{D}\mathrm{\Delta }U+L_{\mu }\left(U_t\right).$(A2) ) can be rewritten in the abstract form as, (A3) ${\dot{U}}_t=A_{\mu }{U}_t+R(\mu ,U_t),$(A3) where $R(\mu ,U_t)=\{\begin{array}{ll}0,& \theta \in [-1,0),\\ f(\mu ,U_t),& \theta =0.\end{array}$Define $f_j=({\text{β}}_j^1,{\text{β}}_j^2,{\text{β}}_j^3)$, where $\begin{array}{rl}{\text{β}}_j^1& =(b_j,0,0{)}^T,{\text{β}}_j^2=(0,b_j,0{)}^T,{\text{β}}_j^3=(0,0,b_j{)}^T,\\ b_j& =\frac{\cos \left(jx/l\right)}{\parallel \cos \left(jx/l\right)\parallel },\parallel \cos \left(jx/l\right)\parallel ={({\int }_0^{l\pi }{\cos }^2(jx/l\left)\mathrm{d}x\right)}^{\frac{1}{2}}.\end{array}$Define $\begin{array}{rl}c\cdot f_j& =c_1{\text{β}}_j^1+c_2{\text{β}}_j^2+c_3{\text{β}}_j^3,c=(c_1,c_2,c_3{)}^T\in \mathcal{C}([-\tau ,0],X),\\ <u,v>& ={\int }_0^{l\pi }[u_1{\overline{v}}_1+u_2{\overline{v}}_2+u_3{\overline{v}}_3]\mathrm{d}x,\end{array}$where $u=(u_1,u_2,u_3),v=(v_1,v_2,v_3),u,v\in X$.

We also define $<\phi ,f_j>=(<{\phi }_1,{\text{β}}_j^1>,<{\phi }_2,{\text{β}}_j^2>,<{\phi }_3,{\text{β}}_j^3>{)}^T.$Further, we define $A_{\mu ,j}{\phi }_j{b}_j=\{\begin{array}{ll}{\dot{\phi }}_j\left(\theta \right)b_j,& \theta \in [-1,0),\\ {\int }_{-r}^0\mathrm{d}{\eta }_j(\mu ,\theta ){\phi }_j\left(\theta \right)b_j,& \theta =0.\end{array}$Thus, we have (A4) $-{\mu }_j\tilde{D}{\phi }_j\left(0\right)+L_{\mu ,j}\left({\phi }_j\right)={\int }_{-1}^0\mathrm{d}{\eta }_j(\mu ,\theta ){\phi }_j\left(\theta \right),$(A4) where ${\eta }_j(\mu ,\theta )=\{\begin{array}{ll}-({\tau }^{\ast }+\mu )B_2,& \theta =-1\\ 0,& \theta \in (-1,0)\\ ({\tau }^{\ast }+\mu )(B_1-{\mu }_{j}D),& \theta =0.\end{array}$For ${\mathcal{C}}^{\ast }=\mathcal{C}\left(\right[0,\tau ],X)$, we define a bilinear form $(\cdot ,\cdot )$ on ${\mathcal{C}}^{\ast }\times \mathcal{C}$ as follows: $(\psi ,\phi )=\sum _{k,j=0}^{\mathrm{\infty }}({\psi }_k,{\phi }_j{)}_c{\int }_0^{l\pi }{b}_k{b}_j\mathrm{d}x,$where $\psi =\sum _{j=0}^{\mathrm{\infty }}{\psi }_j{b}_j\in {\mathcal{C}}^{\ast }$, $\phi =\sum _{j=0}^{\mathrm{\infty }}{\phi }_j{b}_j\in \mathcal{C}$, ${\phi }_j\in \mathcal{C}$, ${\psi }_j\in {\mathcal{C}}^{\ast }$. Note that ${\int }_0^{l\pi }{b}_k{b}_j=0$ for $k\ne j$. Thus, we have $(\psi ,\phi )=\sum _{j=0}^{\mathrm{\infty }}({\psi }_j,{\phi }_j{)}_c\mid b_j{\mid }^2,$where $(\cdot ,\cdot {)}_c$ is the bilinear form defined on ${\mathcal{C}}^{\ast }\times \mathcal{C}$, and satisfies the following equation: $({\psi }_j,{\phi }_j{)}_c={\overline{\psi }}_j(0\left){\phi }_j\right(0)-{\int }_{-1}^0{\int }_0^{\theta }{\overline{\psi }}_j(\xi -\theta \left)\mathrm{d}{\eta }_j\right(0,\theta \left){\phi }_j\right(\xi )\mathrm{d}\xi .$Now, we define an adjoint operator $A^{\ast }$ as follows: $A^{\ast }\psi \left(s\right)=\{\begin{array}{ll}-\dot{\psi }\left(s\right),& s\in (0,1],\\ \sum _{j=0}^{\mathrm{\infty }}{\int }_{-1}^0\mathrm{d}{\eta }_j(0,t){\phi }_j(-t)b_j,& s=0.\end{array}$Let $q\left(\theta \right)b_j=(1,q_1,q_2{)}^T{\mathrm{e}}^{i{\omega }_j{\tau }^{\ast }\theta }{b}_j$ and $q^{\ast }\left(s\right)b_j=M(q_3,q_4,1){\mathrm{e}}^{i{\omega }_j{\tau }^{\ast }s}{b}_j$ be the eigenfunctions of the operators A and $A^{\ast }$ corresponding to the eigenvalues $i{\omega }_j$ and $-i{\omega }_j$, respectively. Thus, one can obtain $\begin{array}{rl}\left(\begin{array}{ccc}i{\omega }_j-a_{11}+\frac{j^2}{l^2}{d}_1& -a_{12}& 0\\ -a_{21}{\mathrm{e}}^{-i{\omega }_j{\tau }^{\ast }}& i{\omega }_j-a_{22}+\frac{j^2}{l^2}{d}_2& -a_{23}\\ 0& -a_{32}& i{\omega }_j+\frac{j^2}{l^2}{d}_3\end{array}\right)\left(\begin{array}{c}1\\ q_1\\ q_2\end{array}\right)& =0,\\ \left(\begin{array}{ccc}i{\omega }_j+a_{11}-\frac{j^2}{l^2}{d}_1& a_{21}{\mathrm{e}}^{-i{\omega }_j{\tau }^{\ast }}& 0\\ a_{12}& i{\omega }_j+a_{22}-\frac{j^2}{l^2}{d}_2& a_{32}\\ 0& a_{23}& i{\omega }_j-\frac{j^2}{l^2}{d}_3\end{array}\right)\left(\begin{array}{c}{q}_3\\ q_4\\ 1\end{array}\right)& =0.\end{array}$By a direct calculation, we get $\begin{array}{rl}{q}_1& =\frac{i{\omega }_j{l}^2-a_{11}{l}^2+j^2{d}_1}{a_{12}{l}^2},q_2=\frac{a_{32}{l}^2{q}_1}{i{\omega }_j{l}^2+j^2{d}_3};\\ q_3& =\frac{l^2{a}_{21}{\mathrm{e}}^{-i{\omega }_j{\tau }^{\ast }}{q}_4}{j^2{d}_1-i{\omega }_j{l}^2-a_{11}{l}^2},q_4=\frac{j^2{d}_3-i{\omega }_j{l}^2}{a_{23}{l}^2}.\end{array}$We choose an M such that $(q^{\ast },q{)}_c=1,(q^{\ast },\overline{q}{)}_c=0$. Thus, we have $M=(q_3+{\overline{q}}_1{q}_4+{\overline{q}}_2+{\tau }^{\ast }{\mathrm{e}}^{i{\omega }_j{\tau }^{\ast }}{a}_{21}{q}_4{)}^{-1}.$This yields the decomposition $\mathcal{C}=P⨁Q$ by ${\mathrm{\Lambda }}_j$ with $P=\{zq{b}_j+\overline{z}\overline{q}{b}_j\mid z\in \mathcal{C}\},Q=\{\phi \in \mathcal{C}\mid ({\overline{q}}^{\ast }{b}_j,\phi )=0,(q^{\ast }{b}_j,\phi )=0\}.$One can rewrite the solution of Equation (A1(A1) $\frac{\mathrm{d}U}{\mathrm{d}t}=\tilde{D}\mathrm{\Delta }U+L_{\mu }\left(U_t\right)+f(\mu ,U_t),$(A1) ) as follows: (A5) $U_t=(qz+\overline{q}\overline{z})\cdot f_j+W(z,\overline{z}),$(A5) where $W=(W^{\left(1\right)},W^{\left(2\right)},W^{\left(3\right)}{)}^T\in Q$. On the centre manifold ${\mathcal{C}}_0$, we have (A6) $W(t,\theta )=W_{20}\left(\theta \right)\frac{z^2}{2}+W_{11}\left(\theta \right)z\overline{z}+W_{02}\left(\theta \right)\frac{{\overline{z}}^2}{2}+\cdots .$(A6) At $\tau ={\tau }^{\ast },z$ satisfies the following relation: (A7) $\dot{z}\left(t\right)=i{\omega }_{j}z\left(t\right)+{\overline{q}}^{{\ast }^T}\left(0\right)<f(0,U_t),f_j>.$(A7) Next, we define $f(0,U_t){\mid }_{{\mathcal{C}}_0}=F(0,z,\overline{z})=F_{zz}^{{}^{″}}\frac{z^2}{2}+F_{z\overline{z}}^{{}^{″}}z\overline{z}+F_{\overline{z}\overline{z}}^{{}^{″}}\frac{{\overline{z}}^2}{2}+\cdots .$For our convenience, we rewrite Equation (A7(A7) $\dot{z}\left(t\right)=i{\omega }_{j}z\left(t\right)+{\overline{q}}^{{\ast }^T}\left(0\right)<f(0,U_t),f_j>.$(A7) ) in the following form: (A8) $\dot{z}=i{\omega }_{j}z+g(z,\overline{z}),$(A8) where (A9) $g(z,\overline{z})=g_{20}\frac{z^2}{2}+g_{11}z\overline{z}+g_{02}\frac{{\overline{z}}^2}{2}+\cdots .$(A9) From the above, we obtain $\begin{array}{rl}{g}_{20}& =2\overline{M}{\tau }^{\ast }[-e{q}_1{\overline{q}}_3+\text{β}{q}_1{\overline{q}}_4{\mathrm{e}}^{-i{\omega }_j{\tau }^{\ast }}+(ch{Z}^{\ast }(h+P^{\ast }{)}^{-3}-r){\overline{q}}_4{q}_1^2\\ & -ch(h+P^{\ast }{)}^{-2}{q}_1{q}_2{\overline{q}}_4+(\rho -d\left)h{Z}^{\ast }\right(h+P^{\ast }{)}^{-3}{q}_1^2\\ & +(d-\rho )h(h+P^{\ast }{)}^{-2}{q}_1{q}_2]{\int }_0^{l\pi }{b}_j^3\mathrm{d}x,\\ g_{11}& =\overline{M}{\tau }^{\ast }\{-e(q_1+{\overline{q}}_1){\overline{q}}_3+{\overline{q}}_4[\text{β}({\overline{q}}_1{\mathrm{e}}^{-i{\omega }_j{\tau }^{\ast }}+q_1{\mathrm{e}}^{i{\omega }_j{\tau }^{\ast }})\\ & +\left(ch{Z}^{\ast }\right(h+P^{\ast }{)}^{-3}-r)2q_1{\overline{q}}_1-ch(h+P^{\ast }{)}^{-2}(q_1{\overline{q}}_2+{\overline{q}}_1{q}_2)]\\ & +\left[\right(\rho -d\left)h{Z}^{\ast }\right(h+P^{\ast }{)}^{-3}2q_1{\overline{q}}_1\\ & +(d-\rho )h(h+P^{\ast }{)}^{-2}(q_1{\overline{q}}_2+{\overline{q}}_1{q}_2\left)\right]\}{\int }_0^{l\pi }{b}_j^3\mathrm{d}x,\\ g_{02}& =2\overline{M}{\tau }^{\ast }\{-e{\overline{q}}_1{\overline{q}}_3+{\overline{q}}_4[\text{β}{\overline{q}}_1{\mathrm{e}}^{i{\omega }_j{\tau }^{\ast }}+\left(ch{Z}^{\ast }\right(h+P^{\ast }{)}^{-3}-r){\overline{q}}_1^2\\ & -ch(h+P^{\ast }{)}^{-2}{\overline{q}}_1{\overline{q}}_2]+(\rho -d)h{Z}^{\ast }(h+P^{\ast }{)}^{-3}{\overline{q}}_1^2\\ & +(d-\rho )h(h+P^{\ast }{)}^{-2}{\overline{q}}_1{\overline{q}}_2\}{\int }_0^{l\pi }{b}_j^3\mathrm{d}x,\\ g_{21}& =2\overline{M}{\tau }^{\ast }\left\{{\overline{q}}_3\right(-e)\left[W_{11}^{\left(2\right)}\right(0)+\frac{1}{2}{W}_{20}^{\left(2\right)}(0)+\frac{1}{2}{W}_{20}^{\left(1\right)}(0){\overline{q}}_1+W_{11}^{\left(1\right)}(0\left)q_1\right]{\int }_0^{l\pi }{b}_j^2\mathrm{d}x\\ & +{\overline{q}}_4[\text{β}\left(W_{11}^{\left(1\right)}\right(-1)q_1+\frac{1}{2}{W}_{20}^{\left(1\right)}(-1){\overline{q}}_1+\frac{1}{2}{W}_{20}^{\left(2\right)}(0){\mathrm{e}}^{i{\omega }_j{\tau }^{\ast }}+W_{11}^2(0\left){\mathrm{e}}^{-i{\omega }_j{\tau }^{\ast }}\right)\\ & +\left(ch{Z}^{\ast }\right(h+P^{\ast }{)}^{-3}-r)\left(2q_1{W}_{11}^{\left(2\right)}\right(0)+W_{20}^{\left(2\right)}(0\left){\overline{q}}_1\right)+(-ch(h+P^{\ast }{)}^{-2})\\ & \times \left(q_1{W}_{11}^3\right(0)+\frac{{\overline{q}}_1{W}_{20}^{\left(3\right)}\left(0\right)}{2}+\frac{W_{20}^{\left(2\right)}\left(0\right)}{2}{\overline{q}}_2+W_{11}^{\left(2\right)}(0\left)q_2\right)]{\int }_0^{l\pi }{b}_j^2\mathrm{d}x\\ & +\left[\right(\rho -d\left)h{Z}^{\ast }\right(h+P^{\ast }{)}^{-3}\left(2q_1{W}_{11}^{\left(2\right)}\right(0)+{\overline{q}}_1{W}_{20}^{\left(2\right)}(0\left)\right)+(d-\rho )h(h+P^{\ast }{)}^{-2}\\ & \times \left(q_1{W}_{11}^{\left(3\right)}\right(0)+\frac{1}{2}{W}_{20}^{\left(3\right)}(0){\overline{q}}_1+{\overline{q}}_2\frac{W_{20}^{\left(2\right)}\left(0\right)}{2}+q_2{W}_{11}^{\left(2\right)}(0\left)\right)]{\int }_0^{l\pi }{b}_j^2\mathrm{d}x\}.\end{array}$Now, we will compute $W_{20}\left(\theta \right)$ and $W_{11}\left(\theta \right)$. From Equation (A7(A7) $\dot{z}\left(t\right)=i{\omega }_{j}z\left(t\right)+{\overline{q}}^{{\ast }^T}\left(0\right)<f(0,U_t),f_j>.$(A7) ), we have (A10) $\begin{array}{rl}\dot{W}& ={\dot{U}}_t-\dot{z}q{b}_j-\dot{\overline{z}}\overline{q}{b}_j\\ & =\{\begin{array}{l}AW-2Re\left\{g\right(z,\overline{z}\left)q\right(\theta \left)\right\}{b}_j,\theta \in [-r,0)\\ AW-2Re\left\{g\right(z,\overline{z}\left)q\right(\theta \left)\right\}{b}_j+F,\theta =0\end{array}\\ & =AW+H(z,\overline{z},\theta ),\end{array}$(A10) where (A11) $H(z,\overline{z},\theta )=H_{20}\left(\theta \right)\frac{z^2}{2}+H_{11}\left(\theta \right)z\overline{z}+H_{02}\left(\theta \right)\frac{{\overline{z}}^2}{2}+\cdots .$(A11) Further, we obtain $\begin{array}{rl}{H}_{20}\left(\theta \right)& =\{\begin{array}{ll}-g_{20}q\left(\theta \right)b_j+{\overline{g}}_{02}\overline{q}\left(\theta \right)b_j,& \theta \in [-r,0),\\ -g_{20}q\left(0\right)b_j+{\overline{g}}_{02}\overline{q}\left(0\right)b_j+F_{zz}^{{}^{″}},& \theta =0,\end{array}\\ H_{11}\left(\theta \right)& =\{\begin{array}{ll}-g_{20}q\left(\theta \right)b_j+{\overline{g}}_{11}\overline{q}\left(\theta \right)b_j,& \theta \in [-r,0),\\ -g_{11}q\left(0\right)b_j+{\overline{g}}_{11}\overline{q}\left(0\right)b_j+F_{z\overline{z}}^{{}^{″}},& \theta =0.\end{array}\end{array}$After differentiating both sides of Equation (A6(A6) $W(t,\theta )=W_{20}\left(\theta \right)\frac{z^2}{2}+W_{11}\left(\theta \right)z\overline{z}+W_{02}\left(\theta \right)\frac{{\overline{z}}^2}{2}+\cdots .$(A6) ) and comparing the coefficients with Equation (A10(A10) $\begin{array}{rl}\dot{W}& ={\dot{U}}_t-\dot{z}q{b}_j-\dot{\overline{z}}\overline{q}{b}_j\\ & =\{\begin{array}{l}AW-2Re\left\{g\right(z,\overline{z}\left)q\right(\theta \left)\right\}{b}_j,\theta \in [-r,0)\\ AW-2Re\left\{g\right(z,\overline{z}\left)q\right(\theta \left)\right\}{b}_j+F,\theta =0\end{array}\\ & =AW+H(z,\overline{z},\theta ),\end{array}$(A10) ), we get (A12) $(A_0-2i{\omega }_{0}I)W_{20}\left(\theta \right)=-H_{20}\left(\theta \right),A_0{W}_{11}\left(\theta \right)=-H_{11}\left(\theta \right),$(A12) where $\theta \in [-1,0)$. By the definition of $A_{\mu }$ and Equation (A12(A12) $(A_0-2i{\omega }_{0}I)W_{20}\left(\theta \right)=-H_{20}\left(\theta \right),A_0{W}_{11}\left(\theta \right)=-H_{11}\left(\theta \right),$(A12) ), we have (A13) $\begin{array}{rl}{W}_{20}\left(\theta \right)& =-\frac{g_{20}}{i{\omega }_j{\tau }^{\ast }}q\left(\theta \right)b_j-\frac{{\overline{g}}_{02}}{3i{\omega }_j{\tau }^{\ast }}\overline{q}\left(\theta \right)b_j+E_1{\mathrm{e}}^{2i{\omega }_j{\tau }^{\ast }\theta },\\ W_{11}\left(\theta \right)& =\frac{g_{11}}{i{\omega }_j{\tau }^{\ast }}q\left(\theta \right)b_j-\frac{{\overline{g}}_{11}}{i{\omega }_j{\tau }^{\ast }}\overline{q}\left(\theta \right)b_j+E_2.\end{array}$(A13) Denote by $E_1=\sum _{j=1}^{\mathrm{\infty }}{E}_1^j{b}_j$ and $E_2=\sum _{j=1}^{\mathrm{\infty }}{E}_2^j{b}_j$. Obviously, $E_1^j$ and $E_2^j$ satisfies the following relations: (A14) $\begin{array}{rl}& (A_0-2i{\omega }_j{\tau }^{\ast }I)E_1^j{b}_j{\mathrm{e}}^{2i{\omega }_0{\tau }^{\ast }}{\mid }_{\theta =0}=-<F_{zz}^{{}^{″}},f_j>b_j,\\ & A_0{E}_2^j{b}_j{\mid }_{\theta =0}=-<F_{z\overline{z}}^{{}^{″}},f_j>b_j,j=1,2,\dots ,\end{array}$(A14) where $\begin{array}{rl}{F}_{zz}^{{}^{″}}& =\sum _{j=1}^{\mathrm{\infty }}<F_{zz}^{\mathrm{\prime }},f_j>b_j,F_{z\overline{z}}^{{}^{″}}=\sum _{j=1}^{\mathrm{\infty }}<F_{z\overline{z}}^{\mathrm{\prime }},f_j>b_j,\\ E_1^j& ={(2i{\omega }_j{\tau }^{\ast }I-{\int }_{-1}^0{\mathrm{e}}^{2i{\omega }_j{\tau }^{\ast }\theta }\mathrm{d}{\eta }_j(0,\theta \left)\right)}^{-1}<F_{zz}^{{}^{″}},f_j>,\\ E_2^j& =-{\left({\int }_{-1}^0\mathrm{d}{\eta }_j\right(0,\theta \left)\right)}^{-1}<F_{z\overline{z}}^{{}^{″}},f_j>.\end{array}$Thus, we have $\begin{array}{rl}{E}_1^j& ={\left(\begin{array}{ccc}2i{\omega }_j{\tau }^{\ast }+\frac{j^2{d}_1}{l^2}-a_{11}& -a_{12}& 0\\ -a_{21}{\mathrm{e}}^{-2i{\omega }_j{\tau }^{\ast }}& 2i{\omega }_j{\tau }^{\ast }+\frac{j^2{d}_2}{l^2}-a_{22}& -a_{23}\\ 0& -a_{32}& 2i{\omega }_j{\tau }^{\ast }+\frac{j^2{d}_3}{l^2}\end{array}\right)}^{-1}<F_{zz}^{{}^{″}},f_j>,\\ E_2^j& ={\left(\begin{array}{ccc}\frac{j^2{d}_1}{l^2}-a_{11}& -a_{12}& 0\\ -a_{21}& \frac{j^2{d}_2}{l^2}-a_{22}& -a_{23}\\ 0& -a_{32}& \frac{j^2{d}_3}{l^2}\end{array}\right)}^{-1}<F_{z\overline{z}}^{{}^{″}},f_j>.\end{array}$From these, one can compute the value of $g_{21}$. Finally, we will compute the values of following quantities: $\begin{array}{rl}{C}_1\left(0\right)& =\frac{i}{2{\omega }_j}(g_{20}{g}_{11}-2\mid g_{11}{\mid }^2-\frac{1}{3}\mid g_{02}{\mid }^2)+\frac{g_{21}}{2},\\ {\mu }_2& =-\frac{Re\left\{C_1\right(0\left)\right\}}{Re{\lambda }^{\mathrm{\prime }}\left(0\right)},\\ {\text{β}}_2& =2Re\left\{C_1\right(0\left)\right\},\\ T_2& =-\frac{Im\left\{C_1\right(0\left)\right\}+{\mu }_{2}Im{\lambda }^{\mathrm{\prime }}\left(0\right)}{{\omega }_j}.\end{array}$The signs of ${\mu }_2$, ${\text{β}}_2$ and $T_2$ will determine the properties of Hopf bifurcation. Particularly, ${\mu }_2$ signifies the direction of the Hopf bifurcation; for ${\mu }_2>0$ $(<0)$, the Hopf bifurcation is supercritical (subcritical). Further, ${\text{β}}_2$ determines the stability of the bifurcation periodic solutions; the solutions are stable (unstable) when ${\text{β}}_2<0$ $(>0)$. Finally, $T_2$ characterizes the period of the bifurcating periodic solutions; the period increases (decreases) when $T_2>0$ $(<0)$.