International Journal of Differential Equations

International Journal of Differential Equations / 2011 / Article

Research Article | Open Access

Volume 2011 |Article ID 603183 |

Swati Khare, O. P. Misra, Chhatrapal Singh, Joydip Dhar, "Role of Delay on Planktonic Ecosystem in the Presence of a Toxic Producing Phytoplankton", International Journal of Differential Equations, vol. 2011, Article ID 603183, 16 pages, 2011.

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

Academic Editor: Gershon Wolansky
Received11 May 2011
Accepted25 Jul 2011
Published11 Oct 2011


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., 𝑁0) in the system. The loss of nutrient due to leaching is assumed to be given by the term π‘Žπ‘. We take π‘Ž1 as the growth rate of phytoplankton biomass, 𝑀 as the rate of predation of phytoplankton by zooplankton, and 𝑀1 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 𝑐1π‘₯𝑦. We further assume that phytoplankton and zooplankton biomass deplete due to natural mortality at the rate of 𝑏1 and π‘Ž2, 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]: 𝑑𝑁𝑑𝑑=𝑁0βˆ’π‘Žπ‘βˆ’π‘π‘π‘₯+π‘˜π‘1π‘₯,𝑑π‘₯𝑑𝑑=π‘Ž1𝑁π‘₯βˆ’π‘1π‘₯βˆ’π‘€π‘₯𝑦,𝐷+π‘₯𝑑𝑦=𝑀𝑑𝑑1π‘₯𝑦𝐷+π‘₯βˆ’π‘Ž2π‘¦βˆ’π‘1π‘₯𝑦,(2.1) with nonnegative initial conditions 𝑁(0),π‘₯(0),𝑦(0)>0. The above system of equations can be nondimensionalised using the relations π‘₯1=𝑁/𝐷, π‘₯2=π‘₯/𝐷, π‘₯3=𝑦/𝐷, 𝜏=π‘Žπ‘‘, and introducing the new parameters 𝛼=𝑁0/π‘Žπ·, 𝛽=𝑏𝐷/π‘Ž, 𝛾=𝑏1/π‘Ž, 𝛽1=π‘Ž1𝐷/π‘Ž, πœ‡=𝑀/π‘Ž, πœ‡1=𝑀1/π‘Ž, πœ‰=π‘Ž2/π‘Ž, and πœ‰1=𝑐1𝐷/π‘Ž. The non-dimensionalised equations of the above system (2.1) are as follows in which we have replaced 𝜏 by 𝑑 and get 𝑑π‘₯1𝑑𝑑=π›Όβˆ’π‘₯1βˆ’π›½π‘₯1π‘₯2+π‘˜π›Ύπ‘₯2,𝑑π‘₯2𝑑𝑑=𝛽1π‘₯1π‘₯2βˆ’π›Ύπ‘₯2βˆ’πœ‡π‘₯2π‘₯31+π‘₯2,𝑑π‘₯3=πœ‡π‘‘π‘‘1π‘₯2π‘₯31+π‘₯2βˆ’πœ‰π‘₯3βˆ’πœ‰1π‘₯2π‘₯3,(2.2) with initial conditions: π‘₯1(0)>0,π‘₯2(0)>0,π‘₯3(0)>0, where, 𝛼, 𝛽, 𝛾, 𝛽1, πœ‡, πœ‡1, πœ‰, πœ‰1, 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 𝑑π‘₯1𝑑𝑑=π›Όβˆ’π‘₯1βˆ’π›½π‘₯1π‘₯2+π‘˜π›Ύπ‘₯2,𝑑π‘₯2𝑑𝑑=𝛽1π‘₯2ξ€œπ‘‘βˆ’βˆžπ›Ό1ξ€·expβˆ’π›Ό1ξ€Έ(π‘‘βˆ’π‘ )𝑓(𝑠)π‘‘π‘ βˆ’π›Ύπ‘₯2βˆ’πœ‡π‘₯2π‘₯31+π‘₯2,(2.3) where 𝑓(𝑠)=𝑓(π‘₯1)=π‘₯1(𝑠), 𝛼1 is delay parameter, 𝑑π‘₯3=πœ‡π‘‘π‘‘1π‘₯2π‘₯31+π‘₯2βˆ’πœ‰π‘₯3βˆ’πœ‰1π‘₯2π‘₯3,(2.4) putting ξ€œπ‘…(𝑑)=π‘‘βˆ’βˆžπ›Ό1ξ€·expβˆ’π›Ό1(ξ€Έπ‘‘βˆ’π‘ )𝑓(𝑠)𝑑𝑠.(2.5) The above system of delay differential equation can be written as 𝑑π‘₯1𝑑𝑑=π›Όβˆ’π‘₯1βˆ’π›½π‘₯1π‘₯2+π‘˜π›Ύπ‘₯2,𝑑π‘₯2𝑑𝑑=𝛽1π‘₯2π‘…βˆ’π›Ύπ‘₯2βˆ’πœ‡π‘₯2π‘₯31+π‘₯2,𝑑π‘₯3=πœ‡π‘‘π‘‘1π‘₯2π‘₯31+π‘₯2βˆ’πœ‰π‘₯3βˆ’πœ‰1π‘₯2π‘₯3,𝑑𝑅𝑑𝑑=𝛼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 Ξ©1, where Ξ©1=ξ‚»ξ€·π‘₯1,π‘₯2,π‘₯3ξ€Έ,π‘…βˆΆ0≀π‘₯1(𝑑)+π‘₯2(𝑑)+π‘₯3𝛼(𝑑)+𝑅(𝑑)β‰€πœ‚1ξ‚Ό,(3.1)πœ‚1=min{(1βˆ’π›Ό),𝛾(1βˆ’π‘˜π›½1/𝛽),πœ‰,𝛼1}.
Proof. Let us consider a time dependent function: π‘Š1(𝑑)=π‘₯1𝛽(𝑑)+𝛽1π‘₯2(𝑑)+π›½πœ‡π›½1πœ‡1π‘₯3(𝑑)+𝑅(𝑑).(3.2) Clearly, π‘‘π‘Š1=𝑑𝑑𝑑π‘₯1+𝛽𝑑𝑑𝛽1𝑑π‘₯2+π‘‘π‘‘π›½πœ‡π›½1πœ‡1𝑑π‘₯3+𝑑𝑑𝑑𝑅𝑑𝑑.(3.3) Using (2.6) in the above expression we obtain π‘‘π‘Š1=ξ€Ίπ‘‘π‘‘π›Όβˆ’π‘₯1βˆ’π›½π‘₯1π‘₯2+π‘˜π›Ύπ‘₯2ξ€»+𝛽𝛽1𝛽1π‘₯2π‘…βˆ’π›Ύπ‘₯2βˆ’πœ‡π‘₯2π‘₯31+π‘₯2ξ‚Ή+πœ‡π›½πœ‡1𝛽1ξ‚Έπœ‡1π‘₯2π‘₯31+π‘₯2βˆ’πœ‰π‘₯3βˆ’πœ‰1π‘₯2π‘₯3ξ‚Ή+𝛼1ξ€·π‘₯1,ξ€·βˆ’π‘…ξ€Έξ€»=π›Όβˆ’1βˆ’π›Ό1ξ€Έπ‘₯1ξ‚΅π›½βˆ’π›Ύ1βˆ’π‘˜1𝛽𝛽𝛽1π‘₯2βˆ’πœ‰πœ‡π›½πœ‡1𝛽1π‘₯3βˆ’πœ‰1πœ‡π›½πœ‡1𝛽1π‘₯2π‘₯3βˆ’π›Ό1𝑅+π‘…βˆ’π‘₯1𝛽π‘₯2β‰€π›Όβˆ’πœ‚1π‘Š1ξ€·(𝑑)+π‘…βˆ’π‘₯1𝛽π‘₯2,(3.4) where πœ‚1 is chosen as the minimum of {(1βˆ’π›Ό1),𝛾(1βˆ’π‘˜π›½1/𝛽),πœ‰,𝛼1}. Thus, π‘‘π‘Š1𝑑𝑑+πœ‚1π‘Š1≀𝛼+π‘…βˆ’π‘₯1𝛽π‘₯2.(3.5) Now applying the theorem of differential inequalities [22], we obtain 0<π‘Š1(𝑑)β‰€π‘Š1(0)π‘’βˆ’πœ‚1𝑑+π›Όπœ‚1+ξ€œξ€·π‘…βˆ’π‘₯1𝛽π‘₯2π‘’βˆ’πœ‚1𝑑𝑑𝑑,(3.6) as π‘‘β†’βˆž, 𝑅→π‘₯1, which gives 0β‰€π‘Š1β‰€π›Όπœ‚1.(3.7) Hence the solution of the system (2.6) is bounded in Ξ©1.

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: 𝐸1≑(𝛼,0,0,𝛼), (ii)boundary equilibrium: 𝐸2≑(𝛾/𝛽1,(𝛼𝛽1βˆ’π›Ύ)/(π›½βˆ’π‘˜π›½1)𝛾,0,𝛾/𝛽1). The boundary equilibrium 𝐸2 exists if either 𝛾/𝛼<𝛽1<𝛽/π‘˜ or 𝛽/π‘˜<𝛽1<𝛾/𝛼 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: 𝐸3≑(π‘₯βˆ—1,π‘₯βˆ—2,π‘₯βˆ—3,π‘…βˆ—), where π‘₯βˆ—1=𝛼+π‘˜π›Ύπ‘₯βˆ—21+𝛽π‘₯βˆ—2,π‘₯βˆ—3=1πœ‡π›½ξ€Ίξ€·1π‘₯βˆ—1βˆ’π›Ύξ€Έξ€·1+π‘₯βˆ—2ξ€Έξ€»,π‘…βˆ—=𝛼+π‘˜π›Ύπ‘₯βˆ—21+𝛽π‘₯βˆ—2,π‘₯βˆ—2=βˆ’π‘ 2±𝑠22βˆ’4𝑠1𝑠32𝑠1,𝑠1=πœ‰1,𝑠2ξ€·πœ‡=βˆ’1βˆ’ξ€·πœ‰+πœ‰1ξ€Έξ€Έ,𝑠3=πœ‰.(3.8) The positive interior equilibrium 𝐸3 exists if πœ‡1>πœ‰+πœ‰1, 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>𝛾/𝛽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, 𝐸1(𝛼,0,0,𝛼), 𝐸2(𝛾/𝛽,(𝛼𝛽1βˆ’π›Ύ)/(π›½βˆ’π‘˜π›½1)𝛾,0,𝛾/𝛽), and 𝐸3(π‘₯βˆ—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 𝐸1 is𝑉1=⎑⎒⎒⎒⎒⎒⎒⎣0ξ€·βˆ’1βˆ’(π›Όπ›½βˆ’π‘˜π›Ύ)00𝛼𝛽1ξ€Έπ›Όβˆ’π›Ύ0000βˆ’πœ‰0100βˆ’π›Ό1⎀βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦.(4.1)

The eigenvalues of the characteristic equation of 𝑉1 are πœ†1=βˆ’1, πœ†2=(𝛼𝛽1βˆ’π›Ύ), πœ†3=βˆ’πœ‰, and πœ†4=βˆ’π›Ό1. It is seen from these eigenvalues that the equilibrium 𝐸1 is locally asymptotically stable if 𝛽1<𝛾/𝛼, 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 𝐸2 is𝑉2=βŽ‘βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ£βˆ’π›½1(π›Όπ›½βˆ’π‘˜π›Ύ)π›Ύξ€·π›½βˆ’π‘˜π›½1ξ€Έβˆ’ξ€·π›½βˆ’π‘˜π›½1𝛾𝛽1ξ€·0000βˆ’πœ‡π›Όπ›½1ξ€Έβˆ’π›Ύπ›Ύξ€·π›½βˆ’π‘˜π›½1ξ€Έ+𝛼𝛽1ξ€Έπ›½βˆ’π›Ύ1𝛼𝛽1ξ€Έβˆ’π›Ύπ›Ύξ€·π›½βˆ’π‘˜π›½1ξ€Έ00𝑗330𝛼100βˆ’π›Ό1⎀βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦,(4.2) where𝑗33=βˆ’ξ‚Έξ‚€βˆšξ€·πœ‰π›Ύπ›½βˆ’π‘˜π›½1ξ€Έβˆ’βˆšπœ‰1𝛼𝛽1ξ€Έξ‚βˆ’π›Ύ2+ξ‚΅ξ‚€βˆšβˆšπœ‰+πœ‰12βˆ’πœ‡1𝛾𝛼𝛽1βˆ’π›Ύξ€Έξ€·π›½βˆ’π‘˜π›½1ξ€Έξ‚Ήξ€·π›Ύξ€·π›½βˆ’π‘˜π›½1ξ€Έ+𝛼𝛽1π›Ύξ€·βˆ’π›Ύξ€Έξ€Έξ€·π›½βˆ’π‘˜π›½1ξ€Έξ€Έ.(4.3)

The eigenvalues πœ†1, πœ†2, and πœ†3 of the above matrix are the roots of the following cubic equation: πœ†3+𝛽1(π›Όπ›½βˆ’π‘˜π›Ύ)π›Ύξ€·π›½βˆ’π‘˜π›½1ξ€Έ+𝛼1ξƒ­πœ†2+𝛼1𝛽1(π›Όπ›½βˆ’π‘˜π›Ύ)π›Ύξ€·π›½βˆ’π‘˜π›½1ξ€Έπœ†+𝛼1𝛼𝛽1ξ€Έβˆ’π›Ύ=0,(4.4) and the fourth eigenvalue is given by πœ†4=𝑗33.

Clearly, πœ†1, πœ†2, and πœ†3 have negative real parts if (π›Όπ›½βˆ’π‘˜π›Ύ)(π›½βˆ’π‘˜π›½1)>0, 𝛽1>𝛾/𝛼, 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 𝛽1(π›Όπ›½βˆ’π‘˜π›Ύ)(𝛽1(π›Όπ›½βˆ’π‘˜π›Ύ)+𝛼1𝛾(π›½βˆ’π‘˜π›½1))>(π›½βˆ’π‘˜π›½1)2(𝛼𝛽1βˆ’π›Ύ)𝛾2 are satisfied, the eigenvalue πœ†4 is negative when πœ‰+πœ‰1√+2πœ‰βˆšπœ‰1>πœ‡1. Thus, we can state the following theorem.

Theorem 4.1. The second boundary equilibrium 𝐸2 is linearly asymptotically stable if ξ€·(i)(π›Όπ›½βˆ’π‘˜π›Ύ)π›½βˆ’π‘˜π›½1ξ€Έ>0,(ii)𝛽1>𝛾𝛼,(iii)πœ‰+πœ‰1√+2πœ‰βˆšπœ‰1>πœ‡1,(iv)𝛽1𝛽(π›Όπ›½βˆ’π‘˜π›Ύ)1(π›Όπ›½βˆ’π‘˜π›Ύ)+𝛼1π›Ύξ€·π›½βˆ’π‘˜π›½1>ξ€·ξ€Έξ€Έπ›½βˆ’π‘˜π›½1ξ€Έ2𝛼𝛽1ξ€Έπ›Ύβˆ’π›Ύ2,(4.5) are satisfied.

For the sake of convenience, the equilibrium points (π‘₯βˆ—1,π‘₯βˆ—2,π‘₯βˆ—3,π‘…βˆ—) of the system is shifted to new points (𝑛1,𝑛2,𝑛3,𝑛4) through transformations 𝑛1=π‘₯1βˆ’π‘₯βˆ—1, 𝑛2=π‘₯2βˆ’π‘₯βˆ—2, 𝑛3=π‘₯3βˆ’π‘₯βˆ—3, 𝑛4=π‘…βˆ’π‘…βˆ—.

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𝑛2𝑛3𝑛4βŽžβŽŸβŽŸβŽŸβŽŸβŽŸβŽŸβŽ βŽ›βŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽβˆ’ξ€·,𝐴=1+𝛽π‘₯βˆ—2ξ€Έβˆ’ξ€·π›½π‘₯βˆ—1ξ€Έ0βˆ’π‘˜π›Ύ00πœ‡π‘₯βˆ—2π‘₯βˆ—3ξ€·1+π‘₯βˆ—2ξ€Έ2βˆ’πœ‡π‘₯βˆ—21+π‘₯βˆ—2𝛽1π‘₯βˆ—20ξƒ©πœ‡1ξ€·1+π‘₯βˆ—2ξ€Έ2βˆ’πœ‰1ξƒͺπ‘₯βˆ—3𝛼00100βˆ’π›Ό1⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠,βŽ›βŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽœβŽπ΅=βˆ’π›½π‘›1𝑛2𝛽1𝑛2𝑛4+𝑝2𝑝3𝑝24𝑛32βˆ’π‘2𝑝3𝑛22+𝑝2𝑝4𝑛2𝑛3βˆ’π‘2𝑝24𝑛22𝑛3+𝑝2𝑝34𝑛32𝑛3+𝑝3𝑝4πœ‡π‘›22βˆ’πœ‡π‘24𝑝3𝑛32βˆ’πœ‡π‘4𝑛2𝑛3+πœ‡π‘24𝑛22𝑛3𝑝3𝑝5π‘₯βˆ—2𝑛22βˆ’π‘3𝑝4𝑝5π‘₯βˆ—2𝑛32βˆ’π‘5π‘₯βˆ—2𝑛2𝑛3+𝑝4𝑝5π‘₯βˆ—2𝑛22𝑛3βˆ’π‘24𝑝5π‘₯βˆ—2𝑛32𝑛3βˆ’π‘5π‘₯βˆ—3𝑛22+𝑝5𝑝3𝑛32βˆ’πœ‰1𝑛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 πœ†4+𝐴1πœ†3+𝐴2πœ†2+𝐴3πœ†+𝐴4=0,(4.8) where𝐴1=𝑝1βˆ’π‘2𝑝3,𝐴2=𝑔2𝛼1+𝑝2𝑝5π‘₯βˆ—3βˆ’π‘2𝑝1𝑝3+πœ‰1π‘₯βˆ—3ξ€Έ,𝐴3=𝛼1𝛽1𝑝6π‘₯βˆ—2+𝑝1𝑝2𝑝5π‘₯βˆ—3βˆ’π‘2𝑝1πœ‰1π‘₯βˆ—3+𝑔2𝑝3𝛼1ξ€Έ,𝐴4=𝑝2𝑝5βˆ’πœ‰1𝛼1𝑔2π‘₯βˆ—3,𝑝1=1+𝛽π‘₯βˆ—2+𝛼1,𝑝2=πœ‡π‘₯βˆ—21+π‘₯βˆ—2,𝑝3=π‘₯βˆ—31+π‘₯βˆ—2,𝑝4=11+π‘₯βˆ—2𝑔2=1+𝛽π‘₯βˆ—2,𝑝5=πœ‡1ξ€·1+π‘₯βˆ—2ξ€Έ2,𝑝6=π›Όπ›½βˆ’π‘˜π›Ύ1+𝛽π‘₯βˆ—2.(4.9)

Using the Routh-Hurwitz criteria [21, 23], we derive that the equilibrium 𝐸3 is locally asymptotically stable, if 𝐴1>0, 𝐴3>0, 𝐴1𝐴2>𝐴3, and 𝐴1𝐴2𝐴3βˆ’π΄23βˆ’π΄21𝐴4>0. Here the conditions 𝐴1>0, 𝐴3>0, 𝐴1𝐴2>𝐴3, and 𝐴1𝐴2𝐴3βˆ’π΄23βˆ’π΄21𝐴4>0 requires πœ‡1𝑝1>πœ‰1πœ‡π‘₯βˆ—2π‘₯βˆ—3,π‘Ÿ1+π‘Ÿ3>π‘Ÿ2+π‘Ÿ4,𝑝1𝑝4𝛼1+𝐿1>𝑝22𝑝3𝑝5π‘₯βˆ—3+𝐿2,π‘Ÿ1𝐿1+π‘Ÿ2𝐿3+π‘Ÿ3𝐿4+π‘Ÿ4𝐿2>π‘Ÿ1𝐿3+π‘Ÿ2𝐿1+π‘Ÿ3𝐿3+π‘Ÿ4𝐿4,(4.10) where π‘Ÿ1=𝑝1𝑝2𝑝5π‘₯βˆ—3, π‘Ÿ2=𝑝1𝑝2πœ‰1π‘₯βˆ—3, π‘Ÿ3=𝛼1𝛽1𝑝6π‘₯βˆ—2, π‘Ÿ4=𝑝2𝑝3𝑔2𝛼1, 𝐿1=𝑝1𝑝22𝑝3+𝑝22𝑝3πœ‰1π‘₯βˆ—3, 𝐿2=𝑝21𝑝2𝑝3+𝛼1𝛽1𝑝6π‘₯βˆ—2, 𝐿3=𝐿2+𝑝22𝑝3𝑝5π‘₯βˆ—3, 𝐿4=𝐿1+𝑝1𝑔2𝛼1.

Thus, we can state the following theorem.

Theorem 4.2. The interior equilibrium 𝐸3 if it exists is linearly asymptotically stable when (i)πœ‡1𝑝1>πœ‰1πœ‡π‘₯βˆ—2π‘₯βˆ—3,(ii)π‘Ÿ1+π‘Ÿ3>π‘Ÿ2+π‘Ÿ4,(iii)𝑝1𝑔2𝛼1+𝐿1>𝑝22𝑝3𝑝5π‘₯βˆ—3+𝐿2,(iv)π‘Ÿ1𝐿1+π‘Ÿ2𝐿3+π‘Ÿ3𝐿4+π‘Ÿ4𝐿2>π‘Ÿ1𝐿3+π‘Ÿ2𝐿1+π‘Ÿ3𝐿3+π‘Ÿ4𝐿4,(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) 𝐴1(π›Όβˆ—)𝐴2(π›Όβˆ—)>𝐴3(π›Όβˆ—) (iii) 𝐴1(π›Όβˆ—)𝐴2(π›Όβˆ—)𝐴3(π›Όβˆ—)βˆ’π΄23(π›Όβˆ—)βˆ’π΄21(π›Όβˆ—)𝐴4(π›Όβˆ—)=0, and (iv) the eigenvalues of the characteristic equation (4.8) should be of the form πœ†π‘–=𝑒𝑖+𝑖𝑣𝑖, where 𝑑𝑒𝑖/𝑑𝛼≠0, 𝑖=1,2,3,4. The condition 𝐴1𝐴2𝐴3βˆ’π΄23βˆ’π΄21𝐴4=0 becomes βˆ’πΊ1𝛼3+𝐺2𝛼2+𝐺3𝛼+𝐺4=0,(4.12) where 𝐺1=β„Ž1+𝑓1, 𝐺2=β„Ž2βˆ’π‘ž25βˆ’π‘“2, 𝐺3=β„Ž3βˆ’2π‘ž5π‘ž6βˆ’π‘“3, 𝐺4=β„Ž4βˆ’π‘ž26βˆ’π‘“4, 𝑓1=π‘ž22π‘ž7, 𝑓2=π‘ž22π‘ž8βˆ’2π‘ž1π‘ž2π‘ž7, 𝑓3=π‘ž21π‘ž7βˆ’2π‘ž1π‘ž2π‘ž8, 𝑓4=π‘ž21π‘ž8, β„Ž1=𝑙2π‘ž5, β„Ž2=𝑙1π‘ž5βˆ’π‘™2π‘ž6, β„Ž3=𝑙1π‘ž6+𝑙3π‘ž5, β„Ž4=𝑙3π‘ž6, π‘ž1=𝑝1βˆ’π‘˜2𝑔6, π‘ž2=π‘˜2𝑔5, π‘ž3=(𝑝2π‘˜6βˆ’π‘1π‘˜2)𝑔5, π‘ž4=(𝑝2π‘˜6βˆ’π‘1π‘˜2)𝑔6+𝑔2𝛼1, π‘ž5=(𝑝2π‘˜6𝑝1βˆ’π‘”2𝛼1π‘˜2)𝑔5+𝛽𝛽1𝛼1𝑔3π‘₯βˆ—2, π‘ž6=(𝑝2π‘˜6𝑝1βˆ’π‘”2𝛼1π‘˜2)𝑔6+(𝛽𝑔4βˆ’π‘˜π›Ύ)𝛽1𝛼1π‘₯βˆ—2, π‘ž7=𝑝2𝑔2𝛼1π‘˜6𝑔5, π‘ž8=𝑝2π‘˜6𝑔2𝛼1𝑔6, 𝑔1=π‘˜π›Ύπ‘₯βˆ—2, 𝑔3=1/𝑔2, 𝑔4=𝑔1/𝑔2, 𝑔5=(1+π‘₯βˆ—2)𝑔3𝛽1/πœ‡1, 𝑔6=(1+π‘₯βˆ—2)(𝑔4𝛽1βˆ’π›Ύ)/πœ‡, π‘˜2=πœ‡π‘₯βˆ—2/(1+π‘₯βˆ—2)2, π‘˜6=πœ‡1/(1+π‘₯βˆ—2)2βˆ’πœ‰1. 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 (𝑒+𝑖𝑣)4+𝐴1(𝑒+𝑖𝑣)3+𝐴2(𝑒+𝑖𝑣)2+𝐴3(𝑒+𝑖𝑣)+𝐴4=0.(4.13) On separating the real and imaginary parts and eliminating 𝑣 between real and imaginary parts, we have 𝑒4+𝑣4βˆ’6𝑒2𝑣2+𝐴1𝑒3βˆ’3𝐴1𝑒𝑣2+𝐴2𝑒2βˆ’π‘£2ξ€Έ+𝐴3𝑒+𝐴4=0,(4.14)βˆ’4𝑒𝑣3+4𝑒3π‘£βˆ’π΄1𝑣3+3𝐴1𝑒2𝑣+2𝐴2𝑒𝑣+𝐴3𝑣=0.(4.15) Substituting the value of 𝑣2 from (4.15) in (4.14), we get𝑒4ξ€·4𝑒+𝐴1ξ€Έ2+ξ€·4𝑒3+3𝐴1𝑒2+2𝐴2𝑒+𝐴3ξ€Έ2+ξ€·4𝑒3+3𝐴1𝑒2+2𝐴2𝑒+𝐴3ξ€ΈΓ—ξ€·βˆ’6𝑒2βˆ’3𝐴1π‘’βˆ’π΄2ξ€Έξ€·4𝑒+𝐴1ξ€Έ+𝐴1𝑒3+𝐴2𝑒2+𝐴3𝑒+𝐴4ξ€Έξ€·4𝑒+𝐴1ξ€Έ2=0,(4.16) differentiating with respect to 𝛼 and putting 𝛼=π›Όβˆ—,we have 𝑑𝑒𝑑𝛼𝛼=π›Όβˆ—=𝐴(βˆ’π‘‘/𝑑𝛼)1𝐴2𝐴3βˆ’π΄23βˆ’π΄21𝐴4ξ€Έ4𝐴3𝐴2βˆ’π΄3ξ€Έβˆ’2𝐴1𝐴22βˆ’π΄1𝐴3ξ€Έ+8𝐴1𝐴2β‰ 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 𝐴1ξ€·π›Όβˆ—ξ€Έπ΄2ξ€·π›Όβˆ—ξ€Έπ΄3ξ€·π›Όβˆ—ξ€Έβˆ’π΄23ξ€·π›Όβˆ—ξ€Έ2βˆ’π΄21ξ€·π›Όβˆ—ξ€Έπ΄4ξ€·π›Όβˆ—ξ€Έξ‚ƒ=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 𝐴1𝐴2𝐴3βˆ’π΄23βˆ’π΄21𝐴4=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=±𝑣𝑖,πœ†3=βˆ’π‘š1,πœ†4=βˆ’π‘š2,(4.20) where βˆšπ‘£=𝐴3/√𝐴1, π‘š1=(βˆ’π΄21βˆ’βˆšπ΄1𝐴31βˆ’4𝐴1𝐴2+4𝐴3)/2𝐴1, and π‘š2=(βˆ’π΄21+√𝐴1𝐴31βˆ’4𝐴1𝐴2+4𝐴3)/2𝐴1.

Next, we seek a transformation matrix 𝑃 which reduces the matrix 𝐴 to the form π‘ƒβˆ’1βŽ›βŽœβŽœβŽœβŽœβŽœβŽœβŽπ΄π‘ƒ=0βˆ’π‘£00𝑣00000βˆ’π‘š10000βˆ’π‘š2⎞⎟⎟⎟⎟⎟⎟⎠,(4.21) where the nonsingular matrix 𝑃 is given as βŽ›βŽœβŽœβŽœβŽœβŽœβŽœβŽπ‘Žπ‘ƒ=11π‘Ž12π‘Ž13π‘Ž140π‘Ž22π‘Ž23π‘Ž24π‘Ž310π‘Ž31π‘Ž31π‘Ž41π‘Ž42π‘Ž43π‘Ž44⎞⎟⎟⎟⎟⎟⎟⎠,(4.22) where, π‘Ž11=𝑣2π‘βˆ’π‘’π‘‘π›Ό1βˆ’π›Ό1𝑣, π‘Ž12=𝛼1𝑣𝑐+𝑣3+π‘’π‘‘πœ”, π‘Ž13=(𝛼1βˆ’π‘š1)(π‘π‘š1+π‘š21βˆ’π‘’π‘‘), π‘Ž14=(𝛼1βˆ’π‘š2)(π‘π‘š2+π‘š22βˆ’π‘’π‘‘), π‘Ž22=𝛼1𝑣𝛽1π‘₯βˆ—2, π‘Ž23=𝛼1π‘š1𝛽1π‘₯βˆ—2, π‘Ž24=𝛼1π‘š2𝛽1π‘₯βˆ—2, π‘Ž31=𝛼1𝑒𝛽1π‘₯βˆ—2, π‘Ž41=(πœ”2+𝑒𝑑)𝛼1, π‘Ž42=𝑣𝑐𝛼1, π‘Ž43=(π‘’π‘‘βˆ’π‘π‘š1βˆ’π‘š21)𝛼1, π‘Ž44=(π‘’π‘‘βˆ’π‘π‘š2βˆ’π‘š22)𝛼1, 𝑐=πœ‡π‘₯βˆ—2π‘₯βˆ—3/(1+π‘₯βˆ—2)2, 𝑑=βˆ’πœ‡π‘₯βˆ—2/(1+π‘₯βˆ—2), 𝑒=(πœ‡1/(1+π‘₯βˆ—2)2βˆ’πœ‰1)π‘₯βˆ—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 Ξ©=π‘ƒβˆ’1𝐴𝑃 and 𝐹=π‘ƒβˆ’1βŽ›βŽœβŽœβŽœβŽœβŽœβŽœβŽπΉπ‘“=1𝑦1,𝑦2,𝑦3,𝑦4𝐹2𝑦1,𝑦2,𝑦3,𝑦4𝐹3𝑦1,𝑦2,𝑦3,𝑦4𝐹4𝑦1,𝑦2,𝑦3,𝑦4ξ€ΈβŽžβŽŸβŽŸβŽŸβŽŸβŽŸβŽŸβŽ βŽ›βŽœβŽœβŽœβŽœβŽœβŽœβŽπ‘“,𝑓isgivenby𝑓=1𝑦1,𝑦2,𝑦3,𝑦4𝑓2𝑦1,𝑦2,𝑦3,𝑦4𝑓3𝑦1,𝑦2,𝑦3,𝑦4𝑓4𝑦1,𝑦2,𝑦3,𝑦4ξ€ΈβŽžβŽŸβŽŸβŽŸβŽŸβŽŸβŽŸβŽ ,(4.24) where, 𝑓1(𝑦1,𝑦2,𝑦3,𝑦4)=βˆ’π›½(π‘Ž11𝑦1+π‘Ž12𝑦2+π‘Ž13𝑦3+π‘Ž14𝑦4)(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4), 𝑓2(𝑦1,𝑦2,𝑦3,𝑦4) =β€‰βˆ’πœ‡π‘4(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)(π‘Ž31𝑦1+π‘Ž31𝑦3+π‘Ž31𝑦4)+πœ‡π‘24(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)2(π‘Ž31y1+π‘Ž31𝑦3+π‘Ž31𝑦4𝑝)βˆ’2𝑝24(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)2(π‘Ž31𝑦1+π‘Ž31𝑦3+π‘Ž31𝑦4)+𝑝2𝑝34(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)3(π‘Ž31𝑦1+π‘Ž31𝑦3+π‘Ž31𝑦4)+(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)(π‘Ž41𝑦1+π‘Ž42𝑦2+π‘Ž43𝑦3+π‘Ž44𝑦4)𝛽1+𝑝2𝑝4(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)π‘₯βˆ—3+πœ‡π‘24(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)2π‘₯βˆ—3βˆ’π‘2𝑝24(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)2π‘₯βˆ—3βˆ’πœ‡π‘34(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)3π‘₯βˆ—3+𝑝2𝑝34(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)3π‘₯βˆ—3, 𝑓3(𝑦1,𝑦2,𝑦3,𝑦4)=βˆ’π‘2𝑝4(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)(π‘Ž31𝑦1+π‘Ž31𝑦3+π‘Ž31𝑦4)+𝑝2𝑝24(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)2(π‘Ž31𝑦1+π‘Ž31𝑦3+π‘Ž31𝑦4)βˆ’(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)(π‘Ž31𝑦1+π‘Ž31𝑦3+π‘Ž31𝑦4)πœ‰1βˆ’π‘2𝑝34(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)2(π‘Ž31𝑦1+π‘Ž31𝑦3+π‘Ž31𝑦4)π‘₯βˆ—2βˆ’π‘2𝑝4(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)2π‘₯βˆ—3+𝑝2𝑝24(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)2π‘₯βˆ—3+𝑝2𝑝24(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)3π‘₯βˆ—3βˆ’π‘2𝑝34(π‘Ž22𝑦2+π‘Ž23𝑦3+π‘Ž24𝑦4)3π‘₯βˆ—3, 𝑓4(𝑦1,𝑦2,𝑦3,𝑦4)=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 𝑔11=14ξƒ©πœ•2𝐹1πœ•π‘¦21+πœ•2𝐹1πœ•π‘¦22ξƒ©πœ•+𝑖2𝐹2πœ•π‘¦21+πœ•2𝐹2πœ•π‘¦22ξƒͺξƒͺ,(4.25)𝑔02=14ξƒ©πœ•2𝐹1πœ•π‘¦21βˆ’πœ•2𝐹1πœ•π‘¦22πœ•βˆ’22𝐹2πœ•π‘¦1πœ•π‘¦2ξƒ©πœ•+𝑖2𝐹2πœ•π‘¦21βˆ’πœ•2𝐹2πœ•π‘¦22πœ•+22𝐹1πœ•π‘¦1πœ•π‘¦2ξƒͺξƒͺ,(4.26)𝑔20=14ξƒ©πœ•2𝐹1πœ•π‘¦21βˆ’πœ•2𝐹1πœ•π‘¦22πœ•+22𝐹2πœ•π‘¦1πœ•π‘¦2ξƒ©πœ•+𝑖2𝐹2πœ•π‘¦21βˆ’πœ•2𝐹2πœ•π‘¦22πœ•βˆ’22𝐹1πœ•π‘¦1πœ•π‘¦2ξƒͺξƒͺ,(4.27)𝐺21=18ξƒ©πœ•3𝐹1πœ•π‘¦31+πœ•3𝐹2πœ•π‘¦32+πœ•3𝐹2πœ•π‘¦21πœ•π‘¦2+πœ•3𝐹1πœ•π‘¦1πœ•π‘¦22ξƒ©πœ•+𝑖3𝐹2πœ•π‘¦31+πœ•3𝐹1πœ•π‘¦32+πœ•3𝐹1πœ•π‘¦21πœ•π‘¦2+πœ•3𝐹2πœ•π‘¦1πœ•π‘¦22,β„Žξƒͺξƒͺ111=14ξƒ©πœ•2𝐹3πœ•π‘¦21+πœ•2𝐹3πœ•π‘¦22ξƒͺ,β„Ž211=14ξƒ©πœ•2𝐹4πœ•π‘¦21+πœ•2𝐹4πœ•π‘¦22ξƒͺ,β„Ž120=14ξƒ©πœ•2𝐹3πœ•π‘¦21βˆ’πœ•2𝐹3πœ•π‘¦22πœ•βˆ’2𝑖2𝐹3πœ•π‘¦1πœ•π‘¦2ξƒͺ,β„Ž220=14ξƒ©πœ•2𝐹4πœ•π‘¦21βˆ’πœ•2𝐹4πœ•π‘¦22πœ•βˆ’2𝑖2𝐹4πœ•π‘¦1πœ•π‘¦2ξƒͺ,𝑀111=β„Ž111π‘š1,𝑀211=β„Ž211π‘š2,𝑀120=β„Ž120π‘š1+2𝑖𝑀,𝑀220=β„Ž220π‘š2,𝐺+2𝑖𝑀1110=12ξ‚΅πœ•2𝐹1πœ•π‘¦1πœ•π‘¦3+πœ•2𝐹2πœ•π‘¦2πœ•π‘¦3ξ‚΅πœ•+𝑖2𝐹2πœ•π‘¦1πœ•π‘¦3βˆ’πœ•2𝐹1πœ•π‘¦2πœ•π‘¦3,𝐺2110=12ξ‚΅πœ•2𝐹1πœ•π‘¦1πœ•π‘¦4+πœ•2𝐹2πœ•π‘¦2πœ•π‘¦4ξ‚΅πœ•+𝑖2𝐹2πœ•π‘¦1πœ•π‘¦4βˆ’πœ•2𝐹1πœ•π‘¦2πœ•π‘¦4,𝐺1101=12ξ‚΅πœ•2𝐹1πœ•π‘¦1πœ•π‘¦3βˆ’πœ•2𝐹2πœ•π‘¦2πœ•π‘¦3ξ‚΅πœ•+𝑖2𝐹2πœ•π‘¦1πœ•π‘¦3+πœ•2𝐹1πœ•π‘¦2πœ•π‘¦3,𝐺2101=12ξ‚΅πœ•2𝐹1πœ•π‘¦1πœ•π‘¦4βˆ’πœ•2𝐹2πœ•π‘¦2πœ•π‘¦4ξ‚΅πœ•+𝑖2𝐹2πœ•π‘¦1πœ•π‘¦4+πœ•2𝐹1πœ•π‘¦2πœ•π‘¦4,ξ‚Άξ‚Ά(4.28)𝑔21=𝐺21+ξ€·2𝐺1110𝑀111+𝐺1101𝑀120+2𝐺2110𝑀211+𝐺2101𝑀220ξ€Έ.(4.29) Thus, we can determine 𝑔11, 𝑔20, 𝑔02, 𝑔21 from (4.25), (4.26), (4.27), and (4.29), respectively. Thus, we can compute the following values: 𝐢1𝑖(0)=𝑔2𝑀20𝑔11βˆ’||𝑔11||2βˆ’13||𝑔02||2ξ‚„+𝑔212,πœ‡2=βˆ’Re𝐢1(0)Reπœ†ξ…ž(π›Όβˆ—),𝜏2ξ€·=βˆ’Im𝐢1(0)+πœ‡2Imπœ†ξ…žξ€·π›Όβˆ—ξ€Έξ€Έπ‘£π›Όβˆ—,𝛽2=2Re𝐢1(0),(4.30) which determine the qualities of bifurcation periodic solution in the center manifold at the critical value π›Όβˆ—.

Theorem 4.4. The parameter πœ‡2 determines the direction of the Hopf-bifurcation if πœ‡2>0(πœ‡2<0), then the Hopf-bifurcation is supercritical (subcritical) and the bifurcation periodic solutions exist for 𝛼>π›Όβˆ—(𝛼>π›Όβˆ—); 𝛽2 determines the stability of bifurcating periodic solution; the bifurcation periodic solutions are orbitally asymptotically stable (unstable) if 𝛽2<0(𝛽2>0); 𝜏2 determines the periodic of the bifurcating periodic solution; the period increases (decreases) if 𝜏2>0(𝜏2<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,𝛽1=0.25,πœ‡=0.2,πœ‡1=0.19,πœ‰=0.05,πœ‰1=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 𝛼1=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 𝐸3(1.06505,0.3967,0.462745). It follows from Section 4, Theorem 4.3, that π›Όβˆ—=1.1680, the positive equilibrium 𝐸3 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 𝛼1 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, 𝑔11=βˆ’0.0000105695βˆ’0.0000200596𝑖, 𝑔20=0.0000362402βˆ’5.6110710βˆ’6𝑖, 𝑔02=βˆ’0.0000151012+0.0000457302𝑖, 𝑔21=βˆ’6.90210βˆ’8βˆ’4.0204810βˆ’8𝑖, 𝐢1(0)=βˆ’4.0826910βˆ’8+1.6292810βˆ’9𝑖, π‘’ξ…ž(π›Όβˆ—)=0.00115975. It follows from (4.30) that πœ‡2>0 and 𝛽2<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.

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 𝛼1=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 𝛼1 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 (𝛼1) play an important role in changing the dynamical behaviour of the system.


  1. J. Duinker and G. Wefer, β€œDas CO2-problem und die rolle des ozeans,,” Naturwissenschaften, vol. 81, no. 6, pp. 237–242, 1994. View at: Publisher Site | Google Scholar
  2. T. M. Work et al., β€œDomic acid in toxication of brown pelicans and cormorants in Santa Cruz California,” in Toxic Phytoplankton Blooms in the Sea, T. J. Smayda and Y. Shimuza, Eds., vol. 3, pp. 643–649, Elsevier, Amsterdam, The Netherlands, 1993. View at: Google Scholar
  3. K. A. Steidinger et al., β€œPfiesteria piscicida a new toxic dinoflagellate genus and species of the order Dinamoebales,” Journal of Phycology, vol. 32, pp. 157–164, 1996. View at: Publisher Site | Google Scholar
  4. T. G. Nielsen, T. Kiorboe, and P. K. Bjornsen, β€œEffects of a Chrysochromulina polylepis subsurface bloom on the planktonic community,” Marine Ecology Progress Series, vol. 62, pp. 21–35, 1990. View at: Publisher Site | Google Scholar
  5. S. A. Levin and L. Segel, β€œHypothesis for the origin plankton patchiness,” Nature, vol. 259, no. 5545, article 659, 1976. View at: Publisher Site | Google Scholar
  6. S. Uye, β€œImpact of copepod grazing on the red tide flagellate chattonella antiqua,” Marine Biology, vol. 92, pp. 35–43, 1986. View at: Publisher Site | Google Scholar
  7. R. K. Upadhyay and J. Chattopadhyay, β€œChaos to order: role of toxin producing phytoplankton in aquatic systems,” Nonlinear Analysis. Modelling and Control, vol. 10, no. 4, pp. 383–396, 2005. View at: Google Scholar | Zentralblatt MATH
  8. S. R.-J. Jang, J. Baglama, and J. Rick, β€œNutrient-phytoplankton-zooplankton models with a toxin,” Mathematical and Computer Modelling, vol. 43, no. 1-2, pp. 105–118, 2006. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  9. T. Saha and M. Bandyopadhyay, β€œDynamical analysis of toxin producing phytoplankton-zooplankton interactions,” Nonlinear Analysis. Real World Applications, vol. 10, no. 1, pp. 314–332, 2009. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  10. E. J. Buskey and C. J. Hyatt, β€œEffects of the Texas (USA) “brown tide” alga on planktonic grazers,” Marine Ecology Progress Series, vol. 126, no. 1–3, pp. 285–292, 1995. View at: Publisher Site | Google Scholar
  11. F. C. Hansen, β€œTrophic interactions between zooplankton and Phaeocystis cf. globosa,” Helgoländer Meeresuntersuchungen, vol. 49, no. 1–4, pp. 283–293, 1995. View at: Publisher Site | Google Scholar
  12. R. Xu, Q. Gan, and Z. Ma, β€œStability and bifurcation analysis on a ratio-dependent predator-prey model with time delay,” Journal of Computational and Applied Mathematics, vol. 230, no. 1, pp. 187–203, 2009. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  13. N. Bairagi and D. Jana, β€œOn the stability and Hopf bifurcation of a delay-induced predator-prey system with habitat complexity,” Applied Mathematical Modelling, vol. 35, no. 7, pp. 3255–3267, 2011. View at: Publisher Site | Google Scholar
  14. S. Khare, O. P. Misra, and J. Dhar, β€œRole of toxin producing phytoplankton on a plankton ecosystem,” Nonlinear Analysis. Hybrid Systems, vol. 4, no. 3, pp. 496–502, 2010. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  15. J. Dhar, A. K. Sharma, and S. Tegar, β€œThe role of delay in digestion of plankton by fish population: a fishery model,” Journal of Nonlinear Science and Its Applications, vol. 1, no. 1, pp. 13–19, 2008. View at: Google Scholar | Zentralblatt MATH
  16. J. Chattopadhyay, R. R. Sarkar, and A. El Abdllaoui, β€œA delay differential equation model on harmful algal blooms in the presence of toxic substances,” IMA Journal of Mathematics Applied in Medicine and Biology, vol. 19, no. 2, pp. 137–161, 2002. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  17. N. MacDonald, Time Lags in Biological Models, vol. 27 of Lecture Notes in Biomathematics, Springer, Berlin, Germany, 1978.
  18. Y. Tang and L. Zhou, β€œGreat time delay in a system with material cycling and delayed biomass growth,” IMA Journal of Applied Mathematics, vol. 70, no. 2, pp. 191–200, 2005. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  19. Y. Tang and L. Zhou, β€œStability switch and Hopf bifurcation for a diffusive prey-predator system with delay,” Journal of Mathematical Analysis and Applications, vol. 334, no. 2, pp. 1290–1307, 2007. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  20. M. Kot, Elements of Mathematical Ecology, Cambridge University Press, Cambridge, UK, 2001. View at: Publisher Site
  21. R. M. May, Stability and Complexity in Model Ecosystem, Princeton University Press, Princeton, NJ, USA, 2001.
  22. G. Birkhoff and G. C. Rota, Ordinary Differential Equation, Ginn, Boston, Mass, USA, 1982.
  23. J. D. Murray, Mathematical Biology, Springer, Berlin, Germany, 2001.
  24. Brian D. Hassard, N. D. Kazarinoff, and Y. H. Wan, Theory and Applications of Hopf Bifurcation, vol. 41 of London Mathematical Society Lecture Note Series, Cambridge University Press, Cambridge, UK, 1981. View at: Zentralblatt MATH
  25. J. E. Marsden and M. McCracken, The Hopf Bifurcation and Its Applications, Springer, New York, NY, USA, 1976.

Copyright © 2011 Swati Khare et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

We are committed to sharing findings related to COVID-19 as quickly as possible. We will be providing unlimited waivers of publication charges for accepted research articles as well as case reports and case series related to COVID-19. Review articles are excluded from this waiver policy. Sign up here as a reviewer to help fast-track new submissions.