Table of Contents Author Guidelines Submit a Manuscript
International Journal of Differential Equations
Volume 2011, Article ID 603183, 16 pages
http://dx.doi.org/10.1155/2011/603183
Research Article

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

1School of Mathematics and Allied Sciences, Jiwaji University, MP Gwalior 474011, India
2Department of Applied Sciences, ABV-Indian Institute of Information Technology and Management, MP Gwalior 474010, India

Received 11 May 2011; Accepted 25 July 2011

Academic Editor: Gershon Wolansky

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.

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 [79]. 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 [1517]. 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𝑡𝛼1exp𝛼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 𝑅(𝑡)=𝑡𝛼1exp𝛼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±𝑠224𝑠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=01(𝛼𝛽𝑘𝛾)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𝛾𝛽10000𝜇𝛼𝛽1𝛾𝛾𝛽𝑘𝛽1+𝛼𝛽1𝛽𝛾1𝛼𝛽1𝛾𝛾𝛽𝑘𝛽100𝑗330𝛼100𝛼1,(4.2) where𝑗33=𝜉𝛾𝛽𝑘𝛽1𝜉1𝛼𝛽1𝛾2+𝜉+𝜉12𝜇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>𝛽𝑘𝛽12𝛼𝛽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𝛽𝑥10𝑘𝛾00𝜇𝑥2𝑥31+𝑥22𝜇𝑥21+𝑥2𝛽1𝑥20𝜇11+𝑥22𝜉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𝑛30.(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=𝜇11+𝑥22,𝑝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=32𝑞5𝑞6𝑓3, 𝐺4=4𝑞26𝑓4, 𝑓1=𝑞22𝑞7, 𝑓2=𝑞22𝑞82𝑞1𝑞2𝑞7, 𝑓3=𝑞21𝑞72𝑞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+𝑣46𝑢2𝑣2+𝐴1𝑢33𝐴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𝑢44𝑢+𝐴12+4𝑢3+3𝐴1𝑢2+2𝐴2𝑢+𝐴32+4𝑢3+3𝐴1𝑢2+2𝐴2𝑢+𝐴3×6𝑢23𝐴1𝑢𝐴24𝑢+𝐴1+𝐴1𝑢3+𝐴2𝑢2+𝐴3𝑢+𝐴44𝑢+𝐴12=0,(4.16) differentiating with respect to 𝛼 and putting 𝛼=𝛼,we have 𝑑𝑢𝑑𝛼𝛼=𝛼=𝐴(𝑑/𝑑𝛼)1𝐴2𝐴3𝐴23𝐴21𝐴44𝐴3𝐴2𝐴32𝐴1𝐴22𝐴1𝐴3+8𝐴1𝐴20.(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𝐴314𝐴1𝐴2+4𝐴3)/2𝐴1, and 𝑚2=(𝐴21+𝐴1𝐴314𝐴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 forṁ𝑌=Ω𝑌+𝐹,(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||213||𝑔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.00001056950.0000200596𝑖, 𝑔20=0.00003624025.61107106𝑖, 𝑔02=0.0000151012+0.0000457302𝑖, 𝑔21=6.9021084.02048108𝑖, 𝐶1(0)=4.08269108+1.62928109𝑖, 𝑢(𝛼)=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.

603183.fig.001
Figure 1: Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.1.
603183.fig.002
Figure 2: Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.1775.
603183.fig.003
Figure 3: Stable phase-space diagram for the system with nutrient, phytoplankton, and zooplankton at 𝛼=1.25.
603183.fig.004
Figure 4: Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.3.
603183.fig.005
Figure 5: Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.1, 𝛼1=0.01.
603183.fig.006
Figure 6: Phase-space diagram for nutrient, phytoplankton, and zooplankton at 𝛼=1.1, 𝛼1=0.07.
603183.fig.007
Figure 7: Stable phase-space diagram for the system with nutrient, phytoplankton, and zooplankton at 𝛼=1.1, 𝛼1=0.12.
603183.fig.008
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 𝛼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.

References

  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 · View at 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 · View at 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 · View at 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 · View at 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 · View at 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 · View at 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 · View at Google Scholar · View at 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 · View at Google Scholar · View at 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 · View at 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 · View at 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 · View at Google Scholar · View at 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 · View at 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 · View at Google Scholar · View at 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 · View at 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 · View at Google Scholar · View at 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 · View at Google Scholar · View at 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 · View at Google Scholar · View at Zentralblatt MATH
  20. M. Kot, Elements of Mathematical Ecology, Cambridge University Press, Cambridge, UK, 2001. View at Publisher · View at Google Scholar
  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.