Role of Delay on Planktonic Ecosystem in the Presence of a Toxic Producing Phytoplankton

Abstract

A mathematical model is proposed to study the role of distributed delay on plankton ecosystem in the presence of a toxic producing phytoplankton. The model includes three state variables, namely, nutrient concentration, phytoplankton biomass, and zooplankton biomass. The release of toxic substance by phytoplankton species reduces the growth of zooplankton and this plays an important role in plankton dynamics. In this paper, we introduce a delay (time-lag) in the digestion of nutrient by phytoplankton. The stability analysis of all the feasible equilibria are studied and the existence of Hopf-bifurcation for the interior equilibrium of the system is explored. From the above analysis, we observe that the supply rate of nutrient and delay parameter play important role in changing the dynamical behaviour of the underlying system. Further, we have derived the explicit algorithm which determines the direction and the stability of Hopf-bifurcation solution. Finally, numerical simulation is carried out to support the theoretical result.

1. Introduction

The study of plankton system is an important area of research in marine ecology. Phytoplankton perform great service for earth. They provide food for marine life, oxygen for human being and also absorb half of the carbon dioxide which may be contributing to the global warming [1]. The dynamics of rapid (massive) increase or almost equal decrease of phytoplankton population is a common feature in marine plankton ecology and known as bloom. Blooms of red tide produce chemical toxin, a type of paralytic poison which can be harmful to zooplankton, fin fish, shellfish, fish, birds, marine mammals, and humans also. Some species, such as the dinoflagellate Alexandrium tamarense and the diatom Pseudo-nitzschia australis [2] produce potent toxins which are liberated into the water before they are eaten, and they may well affect zooplankton when they are in water. It is now well established that a significant number of phytoplankton species produce toxin, such as Pseudo-nitzschia, Gambierdiscus toxicus, Prorocentrum, Ostrepsis, Coolia monotis, Thecadinium, Amphidinium carterae, Dinophysis, Gymnodinium breve, Alexandrium, Gymodinium catenatum, Pyrodinium bahamense, Pfiesteria piscicida, Chrysochromulina polylepis, Prymnesium patelliferum, and P. parvum [3, 4]. Reduction of grazing pressure due to toxin is an important phenomena in plankton ecology [5, 6]. In aquatic system, toxin-producing phytoplankton may act as controlling factor in the phytoplankton-zooplankton interaction dynamics. Efforts have been made to study the role of toxin producing phytoplankton on the phytoplankton-zooplankton dynamics [7–9]. Toxicity may be the strong mediator of zooplankton feeding rate as shown by field studies [4, 10] and laboratory study [11].

Many researchers have been showing keen interest to investigate the direction and stability of Hopf-bifurcation of the system as refer their in [12, 13].

Our current study is motivated by [7, 8, 14], who have considered nutrient interaction with phytoplankton, phytoplankton interaction with zooplankton. The study of interacting species system with nutrient cycling which contributes to the growth of nutrient is carried out in [9, 14, 15]. A model can be more realistic if a delay effect (or time-lag) is being considered in the conversion from one species to another species [15–17]. The prey-predator systems with time delay are deeply considered by [18, 19] and many researchers have used distributed delay in their models [15, 16].

2. Mathematical Model Formulation

Let 𝑁(𝑡) denote the concentration of nutrient at time 𝑡, 𝑥(𝑡) denotes the biomass of toxic producing phytoplankton in the habitat that are partially harmful to zooplankton biomass, 𝑦(𝑡). We assume a constant supply rate of nutrient (i.e., 𝑁₀) in the system. The loss of nutrient due to leaching is assumed to be given by the term 𝑎𝑁. We take 𝑎₁ as the growth rate of phytoplankton biomass, 𝑤 as the rate of predation of phytoplankton by zooplankton, and 𝑤₁ denotes the corresponding conversion rate of zooplankton. The phytoplankton and zooplankton interaction is assumed to follow the Holling type-II functional response [20, 21] with 𝐷 as half saturation constant. Again, the specific rate of nutrient uptake by per unit biomass of phytoplankton in unit time is considered to be 𝑏𝑁 and depletion of zooplankton biomass due to toxin-producing phytoplankton is given by the term 𝑐₁𝑥𝑦. We further assume that phytoplankton and zooplankton biomass deplete due to natural mortality at the rate of 𝑏₁ and 𝑎₂, respectively, and 𝑘 is the fraction of dead phytoplankton biomass that is being recycled back to the nutrient pool. With these assumption and notations, the resultant dynamics of the system under consideration is given by the following set of differential equations [14]: 𝑑𝑁𝑑𝑡=𝑁⁰𝑑𝑥−𝑎𝑁−𝑏𝑁𝑥+𝑘𝑏₁𝑥,𝑑𝑡=𝑎¹𝑁𝑥−𝑏₁𝑥−𝑤𝑥𝑦𝑑𝑦𝐷+𝑥,𝑑𝑡=𝑤¹𝑥𝑦𝐷+𝑥−𝑎²𝑦−𝑐₁𝑥𝑦,(2.1) with nonnegative initial conditions 𝑁(0),𝑥(0),𝑦(0)>0. The above system of equations can be nondimensionalised using the relations 𝑥₁=𝑁/𝐷, 𝑥₂=𝑥/𝐷, 𝑥₃=𝑦/𝐷, 𝜏=𝑎𝑡, and introducing the new parameters 𝛼=𝑁₀/𝑎𝐷, 𝛽=𝑏𝐷/𝑎, 𝛾=𝑏₁/𝑎, 𝛽₁=𝑎₁𝐷/𝑎, 𝜇=𝑤/𝑎, 𝜇₁=𝑤₁/𝑎, 𝜉=𝑎₂/𝑎, and 𝜉₁=𝑐₁𝐷/𝑎. The non-dimensionalised equations of the above system (2.1) are as follows in which we have replaced 𝜏 by 𝑡 and get 𝑑𝑥₁𝑑𝑡=𝛼−𝑥¹𝑑𝑥₂−𝛽𝑥₁𝑥₂+𝑘𝛾𝑥₂,𝑑𝑡=𝛽¹𝑥₁𝑥₂−𝛾𝑥₂−𝜇𝑥²𝑥₃1+𝑥₂𝑑𝑥₃,𝑑𝑡=𝜇¹𝑥₂𝑥₃1+𝑥₂−𝜉𝑥₃−𝜉₁𝑥₂𝑥₃,(2.2) with initial conditions: 𝑥₁(0)>0,𝑥₂(0)>0,𝑥₃(0)>0, where, 𝛼, 𝛽, 𝛾, 𝛽₁, 𝜇, 𝜇₁, 𝜉, 𝜉₁, and 0<𝑘<1 are positive constants.

In this paper, our mathematical model is an extension of the system (2.2) which is studied earlier [14]. Now, we have introduced distributed delay in the digestion of nutrient by phytoplankton. The system (2.2) can be written as 𝑑𝑥₁𝑑𝑡=𝛼−𝑥¹𝑑𝑥₂−𝛽𝑥₁𝑥₂+𝑘𝛾𝑥₂,𝑑𝑡=𝛽¹𝑥₂^𝑡−∞𝛼₁exp−𝛼₁(𝑡−𝑠)𝑓(𝑠)𝑑𝑠−𝛾𝑥₂−𝜇𝑥²𝑥₃1+𝑥₂,(2.3) where 𝑓(𝑠)=𝑓(𝑥₁)=𝑥₁(𝑠), 𝛼₁ is delay parameter, 𝑑𝑥₃𝑑𝑡=𝜇¹𝑥₂𝑥₃1+𝑥₂−𝜉𝑥₃−𝜉₁𝑥₂𝑥₃,(2.4) putting 𝑅(𝑡)=^𝑡−∞𝛼₁exp−𝛼₁(𝑡−𝑠)𝑓(𝑠)𝑑𝑠.(2.5) The above system of delay differential equation can be written as 𝑑𝑥₁𝑑𝑡=𝛼−𝑥¹𝑑𝑥₂−𝛽𝑥₁𝑥₂+𝑘𝛾𝑥₂,𝑑𝑡=𝛽¹𝑥₂𝑅−𝛾𝑥₂−𝜇𝑥²𝑥₃1+𝑥₂𝑑𝑥₃,𝑑𝑡=𝜇¹𝑥₂𝑥₃1+𝑥₂𝑑𝑅−𝜉𝑥₃−𝜉₁𝑥₂𝑥₃,𝑑𝑡=𝛼¹𝑥₁−𝑅,(2.6) with the additional initial condition: 𝑅(0)>0.

3. Boundedness and Equilibria of the System

In this section, we will establish that the system (2.6) is bounded. We begin with the following lemma.

Lemma 3.1. The system (2.6) is uniformly bounded in Ω₁, where Ω₁=𝑥₁,𝑥₂,𝑥₃,𝑅∶0≤𝑥₁(𝑡)+𝑥₂(𝑡)+𝑥₃(𝑡)+𝑅(𝑡)≤𝛼𝜂₁,(3.1)𝜂₁=min{(1−𝛼),𝛾(1−𝑘𝛽₁/𝛽),𝜉,𝛼₁}.
Proof. Let us consider a time dependent function: 𝑊₁(𝑡)=𝑥₁(𝑡)+𝛽𝛽₁𝑥₂(𝑡)+𝛽𝜇𝛽₁𝜇₁𝑥₃(𝑡)+𝑅(𝑡).(3.2) Clearly, 𝑑𝑊₁𝑑𝑡=𝑑𝑑𝑥𝑡¹+𝛽𝛽₁𝑑𝑥₂𝑑𝑡+𝛽𝛽₁𝜇𝜇₁𝑑𝑥₃𝑑𝑡+𝑑𝑑𝑅𝑡.(3.3) Using (2.6) in the above expression we obtain 𝑑𝑊₁𝑑𝑡=𝛼−𝑥¹−𝛽𝑥₁𝑥₂+𝑘𝛾𝑥₂+𝛽𝛽₁𝛽¹𝑥₂𝑅−𝛾𝑥₂−𝜇𝑥²𝑥₃+𝜇𝛽1+𝑥₂𝜇₁𝛽₁𝜇¹𝑥₂𝑥₃1+𝑥₂−𝜉𝑥₃−𝜉₁𝑥₂𝑥₃+𝛼¹𝑥₁=𝛼−1−𝛼₁−𝑅,𝑥₁−𝛾1−𝑘𝛽𝛽¹𝛽𝛽₁𝑥₂−𝜉𝜇𝜇₁𝛽𝛽₁𝑥₃−𝜉₁𝜇𝛽𝜇₁𝛽₁𝑥₂𝑥₃−𝛼₁𝑅+𝑅−𝑥₁𝛽𝑥₂≤𝛼−𝜂₁𝑊₁(𝑡)+𝑅−𝑥₁𝛽𝑥₂,(3.4) where 𝜂₁ is chosen as the minimum of {(1−𝛼₁),𝛾(1−𝑘𝛽₁/𝛽),𝜉,𝛼₁}. Thus, 𝑑𝑊₁𝑑𝑡+𝜂¹𝑊₁≤𝛼+𝑅−𝑥₁𝛽𝑥₂.(3.5) Now applying the theorem of differential inequalities [22], we obtain 0<𝑊₁(𝑡)≤𝑊₁(0)𝑒^−𝜂1𝑡+𝛼𝜂₁+𝑅−𝑥₁𝛽𝑥₂𝑒^−𝜂1𝑡𝑑𝑡,(3.6) as 𝑡→∞, 𝑅→𝑥₁, which gives 0≤𝑊₁≤𝛼𝜂₁.(3.7) Hence the solution of the system (2.6) is bounded in Ω₁.

We now consider the existence of possible equilibria of the system (2.6). The system of (2.6) has three feasible equilibria, namely,(i)boundary equilibrium: 𝐸₁≡(𝛼,0,0,𝛼), (ii)boundary equilibrium: 𝐸₂≡(𝛾/𝛽₁,(𝛼𝛽₁−𝛾)/(𝛽−𝑘𝛽₁)𝛾,0,𝛾/𝛽₁). The boundary equilibrium 𝐸₂ exists if either 𝛾/𝛼<𝛽₁<𝛽/𝑘 or 𝛽/𝑘<𝛽₁<𝛾/𝛼 is satisfied, that is, growth rate of phytoplankton biomass lies between the fraction of natural death rate of phytoplankton biomass to constant supply rate of nutrient and the fraction of nutrient uptake by phytoplankton biomass to fraction of dead phytoplankton biomass,(iii)positive interior equilibrium: 𝐸₃≡(𝑥^∗1,𝑥^∗2,𝑥^∗3,𝑅^∗), where 𝑥^∗1=𝛼+𝑘𝛾𝑥^∗21+𝛽𝑥^∗2,𝑥^∗3=1𝜇𝛽¹𝑥^∗1−𝛾1+𝑥^∗2,𝑅^∗=𝛼+𝑘𝛾𝑥^∗21+𝛽𝑥^∗2𝑥_∗2,=−𝑠²±𝑠²²−4𝑠₁𝑠₃2𝑠₁,𝑠₁=𝜉₁,𝑠₂=−𝜇₁−𝜉+𝜉₁,𝑠₃=𝜉.(3.8) The positive interior equilibrium 𝐸₃ exists if 𝜇₁>𝜉+𝜉₁, that is, the growth rate of zooplankton biomass is greater than the sum of natural death rate and death due to harmful phytoplankton, and also if 𝑥^∗1>𝛾/𝛽₁ means that concentration of nutrient at equilibrium is greater than the fraction of natural death rate of phytoplankton biomass to the depletion rate of nutrient uptake by phytoplankton biomass.

4. Dynamic Behaviour and Hopf-Bifurcation

In the previous section we observed that the system of (2.6) have three equilibria, namely, 𝐸₁(𝛼,0,0,𝛼), 𝐸₂(𝛾/𝛽,(𝛼𝛽₁−𝛾)/(𝛽−𝑘𝛽₁)𝛾,0,𝛾/𝛽), and 𝐸₃(𝑥^∗1,𝑥^∗2,𝑥^∗3,𝑅^∗). We will now examine the dynamical behaviour of the system about all the feasible equilibria.

The variational matrix for the system of (2.6) evaluated at 𝐸₁ is𝑉₁=⎡⎢⎢⎢⎢⎢⎢⎣−01−𝛼(𝛼𝛽𝛽¹−𝑘𝛾)00𝛼00−𝛾−0𝜉00₁00−𝛼₁⎤⎥⎥⎥⎥⎥⎥⎦.(4.1)

The eigenvalues of the characteristic equation of 𝑉₁ are 𝜆₁=−1, 𝜆₂=(𝛼𝛽₁−𝛾), 𝜆₃=−𝜉, and 𝜆₄=−𝛼₁. It is seen from these eigenvalues that the equilibrium 𝐸₁ is locally asymptotically stable if 𝛽₁<𝛾/𝛼, which means that the growth rate of phytoplankton due to the availability of nutrient is less than the fraction of natural death rate of phytoplankton biomass to the constant input rate of nutrient.

The variational matrix for the system of (2.6) evaluated about 𝐸₂ is𝑉₂=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣−𝛽¹(𝛼𝛽−𝑘𝛾)𝛾𝛽−𝑘𝛽₁−𝛽−𝑘𝛽¹𝛾𝛽₁00−𝜇𝛼𝛽0₁0−𝛾𝛾𝛽−𝑘𝛽₁+𝛼𝛽₁−𝛾𝛽¹𝛼𝛽₁−𝛾00𝑗₃₃𝛾𝛽−𝑘𝛽₁𝛼₁000−𝛼₁⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦,(4.2) where𝑗₃₃=−√𝜉𝛾𝛽−𝑘𝛽₁−√𝜉₁𝛼𝛽₁−𝛾²+√𝜉+√𝜉₁²−𝜇₁𝛾𝛼𝛽¹−𝛾𝛽−𝑘𝛽₁𝛾𝛽−𝑘𝛽₁+𝛼𝛽₁−𝛾𝛾𝛽−𝑘𝛽₁.(4.3)

The eigenvalues 𝜆₁, 𝜆₂, and 𝜆₃ of the above matrix are the roots of the following cubic equation: 𝜆³+𝛽¹(𝛼𝛽−𝑘𝛾)𝛾𝛽−𝑘𝛽₁+𝛼¹𝜆²+𝛼¹𝛽₁(𝛼𝛽−𝑘𝛾)𝛾𝛽−𝑘𝛽₁𝜆+𝛼¹𝛼𝛽₁−𝛾=0,(4.4) and the fourth eigenvalue is given by 𝜆₄=𝑗₃₃.

Clearly, 𝜆₁, 𝜆₂, and 𝜆₃ have negative real parts if (𝛼𝛽−𝑘𝛾)(𝛽−𝑘𝛽₁)>0, 𝛽₁>𝛾/𝛼, that is, growth rate of phytoplankton biomass is greater then the fraction of natural death rate of phytoplankton biomass to constant supply rate of nutrient and 𝛽₁(𝛼𝛽−𝑘𝛾)(𝛽₁(𝛼𝛽−𝑘𝛾)+𝛼₁𝛾(𝛽−𝑘𝛽₁))>(𝛽−𝑘𝛽₁)²(𝛼𝛽₁−𝛾)𝛾² are satisfied, the eigenvalue 𝜆₄ is negative when 𝜉+𝜉₁+2√𝜉√𝜉₁>𝜇₁. Thus, we can state the following theorem.

Theorem 4.1. The second boundary equilibrium 𝐸₂ is linearly asymptotically stable if ((iii))𝛽(𝛼₁𝛽−𝑘𝛾)𝛽−𝑘𝛽₁>0,>𝛾(iii)𝜉+𝜉𝛼₁,(iv)𝛽₁+2√𝜉√𝜉₁>𝜇₁,(𝛼𝛽−𝑘𝛾)𝛽₁(𝛼𝛽−𝑘𝛾)+𝛼₁𝛾𝛽−𝑘𝛽₁>𝛽−𝑘𝛽₁²𝛼𝛽₁−𝛾𝛾²,(4.5) are satisfied.

For the sake of convenience, the equilibrium points (𝑥^∗1,𝑥^∗2,𝑥^∗3,𝑅^∗) of the system is shifted to new points (𝑛₁,𝑛₂,𝑛₃,𝑛₄) through transformations 𝑛₁=𝑥₁−𝑥^∗1, 𝑛₂=𝑥₂−𝑥^∗2, 𝑛₃=𝑥₃−𝑥^∗3, 𝑛₄=𝑅−𝑅^∗.

In terms of the new variables, the dynamical equations (2.6) can be written as in matrix form as 𝑋̇=𝐴𝑋+𝐵,(4.6) where dot(⋅) cover 𝑋 denotes the derivative with respect to time. Here 𝐴𝑋 is the linear part of the system and 𝐵 represents the nonlinear part. Moreover,𝑋=⎛⎜⎜⎜⎜⎜⎜⎝𝑛¹𝑛₂𝑛₃𝑛₄⎞⎟⎟⎟⎟⎟⎟⎠,𝐴=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝−1+𝛽𝑥^∗20−𝛽𝜇𝑥𝑥^∗1_∗2−𝑘𝛾00𝑥^∗31+𝑥^∗2²−𝜇𝑥^∗21+𝑥^∗2𝛽₁𝑥^∗20𝜇¹1+𝑥^∗2²−𝜉₁𝑥^∗3𝛼00₁00−𝛼₁𝐵=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝−𝛽𝑛₁⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠,𝑛₂𝛽₁𝑛₂𝑛₄+𝑝₂𝑝₃𝑝²⁴𝑛³²−𝑝₂𝑝₃𝑛²²+𝑝₂𝑝₄𝑛₂𝑛₃−𝑝₂𝑝²⁴𝑛²²𝑛₃+𝑝₂𝑝³⁴𝑛³²𝑛₃𝑝₃+𝑝₃𝑝₄𝜇𝑛²²−𝜇𝑝²⁴𝑝₃𝑛³²−𝜇𝑝₄𝑛₂𝑛₃+𝜇𝑝²⁴𝑛²²𝑛₃𝑝₅𝑥^∗2𝑛²²−𝑝₃𝑝₄𝑝₅𝑥^∗2𝑛³²−𝑝₅𝑥^∗2𝑛₂𝑛₃+𝑝₄𝑝₅𝑥^∗2𝑛²²𝑛₃−𝑝²⁴𝑝₅𝑥^∗2𝑛³²𝑛₃−𝑝₅𝑥^∗3𝑛²²+𝑝₅0𝑝₃𝑛³²−𝜉₁𝑛₂𝑛₃⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠.(4.7) The eigenvalues of the matrix 𝐴 help to understand the stability of the system. The characteristic equation for the variational matrix 𝐴 is given by 𝜆⁴+𝐴₁𝜆³+𝐴₂𝜆²+𝐴₃𝜆+𝐴₄=0,(4.8) where𝐴₃𝐴₁=𝑝₁−𝑝₂𝑝₃,𝐴₂=𝑔₂𝛼₁+𝑝₂𝑝₅𝑥^∗3−𝑝₂𝑝₁𝑝₃+𝜉₁𝑥^∗3,𝑝₁=𝛼₁𝛽₁𝑝₆𝑥^∗2+𝑝₁𝑝₂𝑝₅𝑥^∗3−𝑝₂𝑝₁𝜉₁𝑥^∗3+𝑔₂𝑝₃𝛼₁,𝐴₄=𝑝₂𝑝₅−𝜉₁𝛼₁𝑔₂𝑥^∗3,=1+𝛽𝑥^∗2+𝛼₁,𝑝₂=𝜇𝑥^∗21+𝑥^∗2,𝑝₃=𝑥^∗31+𝑥^∗2,𝑝₄=11+𝑥^∗2𝑝₅𝑔₂=1+𝛽𝑥^∗2,=𝜇¹1+𝑥^∗2²,𝑝₆=𝛼𝛽−𝑘𝛾1+𝛽𝑥^∗2.(4.9)

Using the Routh-Hurwitz criteria [21, 23], we derive that the equilibrium 𝐸₃ is locally asymptotically stable, if 𝐴₁>0, 𝐴₃>0, 𝐴₁𝐴₂>𝐴₃, and 𝐴₁𝐴₂𝐴₃−𝐴²³−𝐴²¹𝐴₄>0. Here the conditions 𝐴₁>0, 𝐴₃>0, 𝐴₁𝐴₂>𝐴₃, and 𝐴₁𝐴₂𝐴₃−𝐴²³−𝐴²¹𝐴₄>0 requires 𝜇₁𝑝₁>𝑟₁𝜉₁𝜇𝑥^∗2𝑥^∗3,𝑟₁+𝑟₃>𝑟₂+𝑟₄,𝑝₁𝑝₄𝛼₁+𝐿₁>𝑝²²𝑝₃𝑝₅𝑥^∗3+𝐿₂,𝐿₁+𝑟₂𝐿₃+𝑟₃𝐿₄+𝑟₄𝐿₂>𝑟₁𝐿₃+𝑟₂𝐿₁+𝑟₃𝐿₃+𝑟₄𝐿₄,(4.10) where 𝑟₁=𝑝₁𝑝₂𝑝₅𝑥^∗3, 𝑟₂=𝑝₁𝑝₂𝜉₁𝑥^∗3, 𝑟₃=𝛼₁𝛽₁𝑝₆𝑥^∗2, 𝑟₄=𝑝₂𝑝₃𝑔₂𝛼₁, 𝐿₁=𝑝₁𝑝²²𝑝₃+𝑝²²𝑝₃𝜉₁𝑥^∗3, 𝐿₂=𝑝²¹𝑝₂𝑝₃+𝛼₁𝛽₁𝑝₆𝑥^∗2, 𝐿₃=𝐿₂+𝑝²²𝑝₃𝑝₅𝑥^∗3, 𝐿₄=𝐿₁+𝑝₁𝑔₂𝛼₁.

Thus, we can state the following theorem.

Theorem 4.2. The interior equilibrium 𝐸₃ if it exists is linearly asymptotically stable when ((iii))𝜇𝑟₁₁𝑝₁>𝜉₁𝜇𝑥^∗2𝑥^∗3,(iii)𝑝₁+𝑟₃>𝑟₂+𝑟₄,(iv)𝑟₁𝑔₂𝛼₁+𝐿₁>𝑝²²𝑝₃𝑝₅𝑥^∗3+𝐿₂,𝐿₁+𝑟₂𝐿₃+𝑟₃𝐿₄+𝑟₄𝐿₂>𝑟₁𝐿₃+𝑟₂𝐿₁+𝑟₃𝐿₃+𝑟₄𝐿₄,(4.11) are satisfied.

Now, we will study the Hopf-bifurcation [24, 25] of the system given by (2.6), taking 𝛼 (i.e., constant input rate of nutrient) as the bifurcation parameter. The necessary and sufficient conditions for the existence of the Hopf-bifurcation for 𝛼=𝛼^∗, if it exists, are (i) 𝐴_𝑖(𝛼^∗)>0, 𝑖=1,2,3,4 (ii) 𝐴₁(𝛼^∗)𝐴₂(𝛼^∗)>𝐴₃(𝛼^∗) (iii) 𝐴₁(𝛼^∗)𝐴₂(𝛼^∗)𝐴₃(𝛼^∗)−𝐴²³(𝛼^∗)−𝐴²¹(𝛼^∗)𝐴₄(𝛼^∗)=0, and (iv) the eigenvalues of the characteristic equation (4.8) should be of the form 𝜆_𝑖=𝑢_𝑖+𝑖𝑣_𝑖, where 𝑑𝑢_𝑖/𝑑𝛼≠0, 𝑖=1,2,3,4. The condition 𝐴₁𝐴₂𝐴₃−𝐴²³−𝐴²¹𝐴₄=0 becomes −𝐺₁𝛼³+𝐺₂𝛼²+𝐺₃𝛼+𝐺₄=0,(4.12) where 𝐺₁=ℎ₁+𝑓₁, 𝐺₂=ℎ₂−𝑞²⁵−𝑓₂, 𝐺₃=ℎ₃−2𝑞₅𝑞₆−𝑓₃, 𝐺₄=ℎ₄−𝑞²⁶−𝑓₄, 𝑓₁=𝑞²²𝑞₇, 𝑓₂=𝑞²²𝑞₈−2𝑞₁𝑞₂𝑞₇, 𝑓₃=𝑞²¹𝑞₇−2𝑞₁𝑞₂𝑞₈, 𝑓₄=𝑞²¹𝑞₈, ℎ₁=𝑙₂𝑞₅, ℎ₂=𝑙₁𝑞₅−𝑙₂𝑞₆, ℎ₃=𝑙₁𝑞₆+𝑙₃𝑞₅, ℎ₄=𝑙₃𝑞₆, 𝑞₁=𝑝₁−𝑘₂𝑔₆, 𝑞₂=𝑘₂𝑔₅, 𝑞₃=(𝑝₂𝑘₆−𝑝₁𝑘₂)𝑔₅, 𝑞₄=(𝑝₂𝑘₆−𝑝₁𝑘₂)𝑔₆+𝑔₂𝛼₁, 𝑞₅=(𝑝₂𝑘₆𝑝₁−𝑔₂𝛼₁𝑘₂)𝑔₅+𝛽𝛽₁𝛼₁𝑔₃𝑥^∗2, 𝑞₆=(𝑝₂𝑘₆𝑝₁−𝑔₂𝛼₁𝑘₂)𝑔₆+(𝛽𝑔₄−𝑘𝛾)𝛽₁𝛼₁𝑥^∗2, 𝑞₇=𝑝₂𝑔₂𝛼₁𝑘₆𝑔₅, 𝑞₈=𝑝₂𝑘₆𝑔₂𝛼₁𝑔₆, 𝑔₁=𝑘𝛾𝑥^∗2, 𝑔₃=1/𝑔₂, 𝑔₄=𝑔₁/𝑔₂, 𝑔₅=(1+𝑥^∗2)𝑔₃𝛽₁/𝜇₁, 𝑔₆=(1+𝑥^∗2)(𝑔₄𝛽₁−𝛾)/𝜇, 𝑘₂=𝜇𝑥^∗2/(1+𝑥^∗2)², 𝑘₆=𝜇₁/(1+𝑥^∗2)²−𝜉₁. Therefore, one pair of eigenvalues of the characteristic equation (4.8) at 𝛼=𝛼^∗ are of the form 𝜆_1,2=±𝑖𝑣, where 𝑣 is positive real number.

Now, we will verify the Hopf-bifurcation condition (iv), putting 𝜆=𝑢+𝑖𝑣 in (4.8), we get (𝑢+𝑖𝑣)⁴+𝐴₁(𝑢+𝑖𝑣)³+𝐴₂(𝑢+𝑖𝑣)²+𝐴₃(𝑢+𝑖𝑣)+𝐴₄=0.(4.13) On separating the real and imaginary parts and eliminating 𝑣 between real and imaginary parts, we have 𝑢⁴+𝑣⁴−−4𝑢6𝑣𝑢²₃𝑣²+𝐴₁𝑢³−3𝐴₁𝑢𝑣²+𝐴₂𝑢²−𝑣²+𝐴₃𝑢+𝐴₄=0,(4.14)+4𝑢³𝑣−𝐴₁𝑣³+3𝐴₁𝑢²𝑣+2𝐴₂𝑢𝑣+𝐴₃𝑣=0.(4.15) Substituting the value of 𝑣² from (4.15) in (4.14), we get𝑢⁴4×𝑢+−𝐴6₁𝑢²₂+4𝑢³+3𝐴₁𝑢²+2𝐴₂𝑢+𝐴₃²+4𝑢³+3𝐴₁𝑢²+2𝐴₂𝑢+𝐴₃−3𝐴₁𝑢−𝐴₂4𝑢+𝐴₁+𝐴₁𝑢³+𝐴₂𝑢²+𝐴₃𝑢+𝐴₄4𝑢+𝐴₁²=0,(4.16) differentiating with respect to 𝛼 and putting 𝛼=𝛼^∗,we have 𝑑𝑑𝛼𝑢_𝛼=𝛼^∗=(−𝑑/𝑑𝛼)𝐴¹𝐴₂𝐴₃−𝐴²³−𝐴²¹𝐴₄4𝐴₃𝐴₂−𝐴₃−2𝐴₁𝐴²²−𝐴₁𝐴₃+8𝐴₁𝐴₂≠0.(4.17) Hence we can state the following theorem.

Theorem 4.3. The system (2.6) has a Hopf-bifurcation at 𝛼^∗>0 such that 𝐴₁𝛼^∗𝐴₂𝛼^∗𝐴₃𝛼^∗−𝐴²³𝛼^∗²−𝐴²¹𝛼^∗𝐴₄𝛼^∗=0,𝑑𝑢𝑑𝛼_𝛼=𝛼^∗≠0.(4.18)

At the Hopf-bifurcation point, the equilibrium state loses its stability and bifurcates to a periodic orbit. We obtain the value of 𝛼 at the Hopf-bifurcation point denoted as 𝛼^∗ and solve the equation 𝐴₁𝐴₂𝐴₃−𝐴²³−𝐴²¹𝐴₄=0.(4.19) At the Hopf-bifurcation point, where the real parts of complex conjugate eigenvalues are zero, the roots of (4.8) are 𝜆_1,2=±𝑣𝑖,𝜆₃=−𝑚₁,𝜆₄=−𝑚₂,(4.20) where 𝑣=√𝐴₃/√𝐴₁, 𝑚₁=(−𝐴²¹−√𝐴₁𝐴³¹−4𝐴₁𝐴₂+4𝐴₃)/2𝐴₁, and 𝑚₂=(−𝐴²¹+√𝐴₁𝐴³¹−4𝐴₁𝐴₂+4𝐴₃)/2𝐴₁.

Next, we seek a transformation matrix 𝑃 which reduces the matrix 𝐴 to the form 𝑃⁻¹𝐴𝑃=⎛⎜⎜⎜⎜⎜⎜⎝𝑣00−00𝑣−00𝑚₁00000−0𝑚₂⎞⎟⎟⎟⎟⎟⎟⎠,(4.21) where the nonsingular matrix 𝑃 is given as 𝑃=⎛⎜⎜⎜⎜⎜⎜⎝𝑎¹¹𝑎₁₂𝑎₁₃𝑎₁₄0𝑎₂₂𝑎₂₃𝑎₂₄𝑎₃₁0𝑎₃₁𝑎₃₁𝑎₄₁𝑎₄₂𝑎₄₃𝑎₄₄⎞⎟⎟⎟⎟⎟⎟⎠,(4.22) where, 𝑎₁₁=𝑣²𝑐−𝑒𝑑𝛼₁−𝛼₁𝑣, 𝑎₁₂=𝛼₁𝑣𝑐+𝑣³+𝑒𝑑𝜔, 𝑎₁₃=(𝛼₁−𝑚₁)(𝑐𝑚₁+𝑚²¹−𝑒𝑑), 𝑎₁₄=(𝛼₁−𝑚₂)(𝑐𝑚₂+𝑚²²−𝑒𝑑), 𝑎₂₂=𝛼₁𝑣𝛽₁𝑥^∗2, 𝑎₂₃=𝛼₁𝑚₁𝛽₁𝑥^∗2, 𝑎₂₄=𝛼₁𝑚₂𝛽₁𝑥^∗2, 𝑎₃₁=𝛼₁𝑒𝛽₁𝑥^∗2, 𝑎₄₁=(𝜔²+𝑒𝑑)𝛼₁, 𝑎₄₂=𝑣𝑐𝛼₁, 𝑎₄₃=(𝑒𝑑−𝑐𝑚₁−𝑚²¹)𝛼₁, 𝑎₄₄=(𝑒𝑑−𝑐𝑚₂−𝑚²²)𝛼₁, 𝑐=𝜇𝑥^∗2𝑥^∗3/(1+𝑥^∗2)², 𝑑=−𝜇𝑥^∗2/(1+𝑥^∗2), 𝑒=(𝜇₁/(1+𝑥^∗2)²−𝜉₁)𝑥^∗3.

To achieve the normal form of (4.6), we make another change of variable, that is, 𝑋=𝑃𝑌, where 𝑌=^𝑦1𝑦2_𝑦3𝑦4.

Through some algebraic manipulations, (4.6) takes the form𝑌̇=Ω𝑌+𝐹,(4.23) where Ω=𝑃⁻¹𝐴𝑃 and 𝐹=𝑃⁻¹𝑓=⎛⎜⎜⎜⎜⎜⎜⎝𝐹¹𝐹₂𝑦₁,𝑦₂,𝑦₃,𝑦₄𝐹₃𝑦₁,𝑦₂,𝑦₃,𝑦₄𝐹₄𝑦₁,𝑦₂,𝑦₃,𝑦₄𝑦₁,𝑦₂,𝑦₃,𝑦₄⎞⎟⎟⎟⎟⎟⎟⎠,𝑓isgivenby𝑓=⎛⎜⎜⎜⎜⎜⎜⎝𝑓¹𝑓₂𝑦₁,𝑦₂,𝑦₃,𝑦₄𝑓₃𝑦₁,𝑦₂,𝑦₃,𝑦₄𝑓₄𝑦₁,𝑦₂,𝑦₃,𝑦₄𝑦₁,𝑦₂,𝑦₃,𝑦₄⎞⎟⎟⎟⎟⎟⎟⎠,(4.24) where, 𝑓¹(𝑦₁,𝑦₂,𝑦₃,𝑦₄)=−𝛽(𝑎₁₁𝑦₁+𝑎₁₂𝑦₂+𝑎₁₃𝑦₃+𝑎₁₄𝑦₄)(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄), 𝑓²(𝑦₁,𝑦₂,𝑦₃,𝑦₄) = 𝑝−𝜇𝑝₄(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)(𝑎₃₁𝑦₁+𝑎₃₁𝑦₃+𝑎₃₁𝑦₄)+𝜇𝑝²⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)²(𝑎₃₁y₁+𝑎₃₁𝑦₃+𝑎₃₁𝑦₄)−₂𝑎₂₃𝑝²⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)²(𝑎₃₁𝑦₁+𝑎₃₁𝑦₃+𝑎₃₁𝑦₄)+𝑝₂𝑝³⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)³(𝑎₃₁𝑦₁+𝑎₃₁𝑦₃+𝑎₃₁𝑦₄)+(𝑎₂₂𝑦₂+𝑝₂𝑦₃+𝑎₂₄𝑦₄)(𝑎₄₁𝑦₁+𝑎₄₂𝑦₂+𝑎₄₃𝑦₃+𝑎₄₄𝑦₄)𝛽₁+𝑝₂𝑝₄(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)𝑥^∗3+𝜇𝑝²⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)²𝑥^∗3−𝑝²⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)²𝑥^∗3−𝜇𝑝³⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)³𝑥^∗3+𝑝₂𝑝³⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)³𝑥^∗3, 𝑎𝑓³₃₁(𝑦₁,𝑦₂,𝑦₃,𝑦₄)=−𝑝₂𝑝₄(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)(𝑎₃₁𝑦₁+𝑎₃₁𝑦₃+𝑎₃₁𝑦₄)+𝑝₂𝑝²⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)²(𝑎₃₁𝑦₁+𝑎₃₁𝑦₃+𝑎₃₁𝑦₄)−(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)(𝑎₃₁𝑦₁+𝑎₃₁𝑦₃+𝑎₃₁𝑦₄)𝜉₁−𝑝₂𝑝³⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)²(𝑎₃₁𝑦₁+𝑎₃₁𝑦₃+𝑝₂𝑦₄)𝑥^∗2−𝑝₂𝑝₄(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)²𝑥^∗3+𝑝₂𝑝²⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)²𝑥^∗3+𝑝₂𝑝²⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)³𝑥^∗3−𝑝³⁴(𝑎₂₂𝑦₂+𝑎₂₃𝑦₃+𝑎₂₄𝑦₄)³𝑥^∗3, 𝑓⁴(𝑦₁,𝑦₂,𝑦₃,𝑦₄)=0.

Equation (4.23) is the normal form of (4.6) from which the stability coefficient can be computed. In (4.6), on the right hand side the first term is linear and the second one is nonlinear in 𝑦’s. For evaluating the direction of bifurcating solution, we can evaluate the following quantities at 𝛼=𝛼^∗ and origin 𝑔₁₁=14𝜕²𝐹¹𝜕𝑦²¹+𝜕²𝐹¹𝜕𝑦²²+𝑖𝜕²𝐹²𝜕𝑦²¹+𝜕²𝐹²𝜕𝑦²²,(4.25)𝑔₀₂=14𝜕²𝐹¹𝜕𝑦²¹−𝜕²𝐹¹𝜕𝑦²²−2𝜕²𝐹²𝜕𝑦₁𝜕𝑦₂+𝑖𝜕²𝐹²𝜕𝑦²¹−𝜕²𝐹²𝜕𝑦²²+2𝜕²𝐹¹𝜕𝑦₁𝜕𝑦₂,(4.26)𝑔₂₀=14𝜕²𝐹¹𝜕𝑦²¹−𝜕²𝐹¹𝜕𝑦²²+2𝜕²𝐹²𝜕𝑦₁𝜕𝑦₂+𝑖𝜕²𝐹²𝜕𝑦²¹−𝜕²𝐹²𝜕𝑦²²−2𝜕²𝐹¹𝜕𝑦₁𝜕𝑦₂,(4.27)𝐺₂₁=18𝜕³𝐹¹𝜕𝑦³¹+𝜕³𝐹²𝜕𝑦³²+𝜕³𝐹²𝜕𝑦²¹𝜕𝑦₂+𝜕³𝐹¹𝜕𝑦₁𝜕𝑦²²+𝑖𝜕³𝐹²𝜕𝑦³¹+𝜕³𝐹¹𝜕𝑦³²+𝜕³𝐹¹𝜕𝑦²¹𝜕𝑦₂+𝜕³𝐹²ℎ𝜕𝑦₁𝜕𝑦²²₁₁₁,=14𝜕²𝐹³𝜕𝑦²¹+𝜕²𝐹³𝜕𝑦²²,ℎ²¹¹=14𝜕²𝐹⁴𝜕𝑦²¹+𝜕²𝐹⁴ℎ₁₂₀𝜕𝑦²²,=14𝜕²𝐹³𝜕𝑦²¹−𝜕²𝐹³𝜕𝑦²²−2𝑖𝜕²𝐹³𝜕𝑦₁𝜕𝑦₂,ℎ²²⁰=14𝜕²𝐹⁴𝜕𝑦²¹−𝜕²𝐹⁴𝜕𝑦²²−2𝑖𝜕²𝐹⁴𝑤₁₁₁𝜕𝑦₁𝜕𝑦₂,=ℎ¹¹¹𝑚₁,𝑤²¹¹=ℎ²¹¹𝑚₂,𝑤¹²⁰=ℎ¹²⁰𝑚₁+2𝑖𝑤,𝑤²²⁰=ℎ²²⁰𝐺𝑚₂+2𝑖𝑤₁₁₁₀,=12𝜕²𝐹¹𝜕𝑦₁𝜕𝑦₃+𝜕²𝐹²𝜕𝑦₂𝜕𝑦₃+𝑖𝜕²𝐹²𝜕𝑦₁𝜕𝑦₃−𝜕²𝐹¹𝐺𝜕𝑦₂𝜕𝑦₃₂₁₁₀,=12𝜕²𝐹¹𝜕𝑦₁𝜕𝑦₄+𝜕²𝐹²𝜕𝑦₂𝜕𝑦₄+𝑖𝜕²𝐹²𝜕𝑦₁𝜕𝑦₄−𝜕²𝐹¹𝐺𝜕𝑦₂𝜕𝑦₄₁₁₀₁,=12𝜕²𝐹¹𝜕𝑦₁𝜕𝑦₃−𝜕²𝐹²𝜕𝑦₂𝜕𝑦₃+𝑖𝜕²𝐹²𝜕𝑦₁𝜕𝑦₃+𝜕²𝐹¹𝐺𝜕𝑦₂𝜕𝑦₃₂₁₀₁,=12𝜕²𝐹¹𝜕𝑦₁𝜕𝑦₄−𝜕²𝐹²𝜕𝑦₂𝜕𝑦₄+𝑖𝜕²𝐹²𝜕𝑦₁𝜕𝑦₄+𝜕²𝐹¹𝜕𝑦₂𝜕𝑦₄,(4.28)𝑔₂₁=𝐺₂₁+2𝐺¹¹¹⁰𝑤¹¹¹+𝐺¹¹⁰¹𝑤¹²⁰+2𝐺²¹¹⁰𝑤²¹¹+𝐺²¹⁰¹𝑤²²⁰.(4.29) Thus, we can determine 𝑔₁₁, 𝑔₂₀, 𝑔₀₂, 𝑔₂₁ from (4.25), (4.26), (4.27), and (4.29), respectively. Thus, we can compute the following values: 𝐶₁(0)=𝑖2𝑤𝑔²⁰𝑔₁₁−||𝑔₁₁||²−13||𝑔⁰²||²𝜇₂+𝑔2²¹,=−Re𝐶¹(0)Re𝜆^(𝛼^∗),𝜏²=−Im𝐶¹(0)+𝜇₂Im𝜆^𝛼^∗𝑣𝛼^∗,𝛽₂=2Re𝐶₁(0),(4.30) which determine the qualities of bifurcation periodic solution in the center manifold at the critical value 𝛼^∗.

Theorem 4.4. The parameter 𝜇₂ determines the direction of the Hopf-bifurcation if 𝜇₂>0(𝜇₂<0), then the Hopf-bifurcation is supercritical (subcritical) and the bifurcation periodic solutions exist for 𝛼>𝛼^∗(𝛼>𝛼^∗); 𝛽₂ determines the stability of bifurcating periodic solution; the bifurcation periodic solutions are orbitally asymptotically stable (unstable) if 𝛽₂<0(𝛽₂>0); 𝜏₂ determines the periodic of the bifurcating periodic solution; the period increases (decreases) if 𝜏₂>0(𝜏₂<0).

5. Numerical Example

In this section, we have performed numerical simulation for both systems (2.2) as well as (2.6) using MATLAB. We are taking parametric values 𝛽=0.3,𝛾=0.2,𝑘=0.3,𝛽₁=0.25,𝜇=0.2,𝜇₁=0.19,𝜉=0.05,𝜉₁=0.01 for both systems. The behavior of the system (2.2) with different values of 𝛼. From Figures 1, 2, 3, 4, it is observed that as 𝛼 increases, the stable system starts oscillating. Again, with these set of values and 𝛼₁=0.12, we get a positive root of (4.12) is 𝛼=𝛼^∗=1.1680. Therefore, one pair of eigenvalues of the characteristic equation (4.8) at 𝛼=𝛼^∗=1.1680 are of the form 𝜆_1,2=±𝑖𝑣, where 𝑣 is positive real number. Here the positive interior equilibrium is 𝐸₃(1.06505,0.3967,0.462745). It follows from Section 4, Theorem 4.3, that 𝛼^∗=1.1680, the positive equilibrium 𝐸₃ is stable when 𝛼<𝛼^∗ (see Figure 7) and the system (2.6) undergoes a Hopf-bifurcation at 𝛼>𝛼^∗ (see Figure 8). Now keeping 𝛼=1.1, further reduction of 𝛼₁ at 0.01 and 0.07 are shown on Figures 5-6. Moreover, we can determine the stability and direction of periodic bifurcation from the positive equilibrium at the critical point 𝛼^∗. For instance, when 𝛼=𝛼^∗=1.1680, 𝑔₁₁=−0.0000105695−0.0000200596𝑖, 𝑔₂₀=0.0000362402−5.6110710⁻6𝑖, 𝑔₀₂=−0.0000151012+0.0000457302𝑖, 𝑔₂₁=−6.90210⁻8−4.0204810⁻8𝑖, 𝐶₁(0)=−4.0826910⁻8+1.6292810⁻9𝑖, 𝑢^(𝛼^∗)=0.00115975. It follows from (4.30) that 𝜇₂>0 and 𝛽₂<0. Therefore, the bifurcation takes place when 𝛼 crosses 𝛼^∗ to the right (𝛼>𝛼^∗), and the corresponding periodic orbits are orbitally asymptotically stable, as depicted in Figure 8.

Figure 1 

Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.1.

Figure 2 

Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.1775.

Figure 3 

Stable phase-space diagram for the system with nutrient, phytoplankton, and zooplankton at 𝛼=1.25.

Figure 4 

Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.3.

Figure 5 

Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.1, 𝛼1=0.01.

Figure 6 

Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.1, 𝛼1=0.07.

Figure 7 

Stable phase-space diagram for the system with nutrient, phytoplankton, and zooplankton at 𝛼=1.1, 𝛼1=0.12.

Figure 8 

Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.1775, 𝛼1=0.12.

6. Conclusion

In this paper, we have studied the role of delay on plankton ecosystem in the presence of a toxic producing phytoplankton. In this system, it has been assumed that the toxic phytoplankton reduces the growth of zooplankton and studied the effect of delay in the digestion of nutrient by phytoplankton biomass. The local stability conditions of all the feasible equilibria of the system are established. The interior equilibrium are locally asymptotically stable under certain conditions. We have shown numerically with the set of parameters that at 𝛼₁=0.01, the system exhibits the chaotic behavior (see Figure 5). Figure 6 exhibits the oscillatory behavior of the system and further it is observed that the increase in 𝛼₁ makes the system stable (see Figure 7). For the comparison of the system (2.6), without delay in the system, a numerical simulation of the system (2.2) is shown in Figures 1, 2, 3, 4. By applying the normal form theory and the center manifold theorem, we define the explicit formulae that determine the stability and direction of the bifurcating periodic solutions. For numerical experiment, it is observed that when the input rate of nutrient, 𝛼, exceeds the critical value 𝛼^∗(1.1680), the system (2.6) leads to oscillatory behaviour shown in Figure 8. Thus, the quantitative level of abundance of system populations depends crucially on the input rate of nutrient. Further from Theorem 4.4, we can determine the direction and stability of Hopf-bifurcation. For these chosen set of parametric values, the Hopf-bifurcation is supercritical and stable. Hence, we conclude that the supply rate (𝛼) of nutrient and delay parameter (𝛼₁) play an important role in changing the dynamical behaviour of the system.