Stability and bifurcation in a two harmful phytoplankton–zooplankton system

Chaos, Solitons and Fractals 39 (2009) 1395–1409 www.elsevier.com/locate/chaos

Stability and bifurcation in a two harmful phytoplankton–zooplankton system q Jiantao Zhao, Junjie Wei

*

Department of Mathematics, Harbin Institute of Technology, Harbin 150001, PR China Accepted 29 May 2007

Abstract In this paper, a mathematical model consisting of two harmful phytoplankton and zooplankton with discrete time delays is considered. We prove that a sequence of Hopf bifurcations occur at the interior equilibrium as the delay increases. Meanwhile, the phenomenon of stability switches is found under certain conditions. The direction of the Hopf bifurcations and the stability of the bifurcating periodic solutions are determined by using the theory of normal form and center manifold. Numerical simulations are given to support the theoretical results.  2007 Elsevier Ltd. All rights reserved.

1. Introduction There has been a global increase in harmful plankton blooms in last three decades [3,12,27]. Recently, there has been increasing interest in investigating the dynamics of harmful algal bloom (HAB) and its control [5,14,26,29]. Plankton model is a kind of ecologic model and different plankton models have been established and studied, previously [2,9,31,37]. A broad classification of HAB species distinguishes two groups: the toxin producers, which can contaminate seafood or kill fish, and the high-biomass producers, which can cause anoxia and indiscriminate mortalities of marine life after reaching dense concentrations. Some HAB species have characteristics of both groups. Phytoplankton blooms may be dispersed by physical factors or removed from the water column due to sinking, while mortality of individual cells within blooms may be effected by autolysis, viruses, predatory bacteria, or grazing zooplankton [6,17,30,32]. Effect of environmental factors and existence of harmful species in a community play an important role in generating hypothesis about the nature of interactions between species in an assemblage. Microzooplankton grazing of harmful algal bloom species (HABs) is important in two ways: first, grazing may serve to prevent blooms or lessen their magnitude; second, for toxin producing species, grazers may accumulate and subsequently transfer toxins [19,33]. There are several papers attempt to establish the role of different hydrological parameters in the formation of plankton blooms and to look for a suitable form of functional response to describe the reduction of zooplankton population due to toxin producing phytoplankton (TPP), for example, see [3,8,12]. The adverse of harmful algal bloom is clear, but the control of such problem is under investigation. Hence studies regarding the pattern of bloom and its control are necessary towards q *

This research is supported by the NNSF of China. Corresponding author. E-mail address: [email protected] (J. Wei).

0960-0779/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2007.05.019

1396

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

its serious ecological problem. Recent studies reveal that some times, bloom of certain harmful species leads to release of both toxins and allelopathic substances, which affect the growth of other microorganisms [4]. Among marine algae, allelopathy has been observed both in vitro and in situ [10], but the chemical nature and the role of allelopathic compounds are still not known [28]. The current study has been originated from the theoretical as well as experimental results on the interaction of harmful algal bloom and different types of phytoplankton–zooplankton interactions [7,8,23–25]. In this paper, we shall discuss a two harmful phytoplankton–zooplankton system proposed by Sarkar et al. [26]. This model is described by the system of differential equations   dP 1 P1 ¼ r1 P 1 1   a1 P 1 P 2  q1 P 1 Z; dt k   dP 2 P2 ð1Þ ¼ r2 P 2 1   a2 P 1 P 2  q2 P 2 Z; dt k dZ ¼ ðc1 P 1 þ c2 P 2 ÞZ  dZ  h1 P 1 ðt  s1 ÞZ  h2 P 2 ðt  s2 ÞZ; dt where P1(t) and P2(t) are the concentrations of the two kind of harmful phytoplankton and Z(t) is the concentration of zooplankton at time t; r1 and r2 are the growth rates of the harmful phytoplankton, respectively, and k is the environmental carrying capacity, which is common for both the plankton species and can be interpreted as competition between the plankton due to limited resource; q1 and q2 are the maximum zooplankton ingestion rates for both the harmful phytoplankton species; c1 and c2 are the maximum zooplankton conversion rates; d is the natural death rate of zooplankton; a1 and a2 are the inhibitory effects of the two competing harmful phytoplankton; h1(c1 > h1) and h2(c2 > h2) are the rates of toxin liberation by the harmful phytoplankton, respectively, which reduces the growth of zooplankton. In [26], the authors have studied the stability of the feasible steady states of system (1). Some numerical results are given. The purpose of the present paper is study system (1) in view of bifurcation. By analyzing the distribution of the characteristic equation associated with system (1) in detail, the phenomenon of the stability switch is found. We not only investigate the stability of the fixed point and the existence of the Hopf bifurcations, but also the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions. The remainder of this paper is organized as following. In Section 2, we employ the method due to Ruan and Wei [21,22] to analyze the distribution of the characteristic equation associated with this model, and obtain the existence of the local Hopf bifurcation. The phenomenon of stability switch is found when the time delay varies. In Section 3, the direction and stability of periodic solutions bifurcating from the Hopf bifurcation are determined by using the normal form theory and center manifold argument presented in [13]. Finally, some numerical simulation examples are given to illustrate the results found. We would like to mention that there are several excellent articles on the study of stability and bifurcation for ecologic model with delay or without delay, we refer to [1,11,15,16,18,20,34–36] and the references therein. 2. Stability of the interior equilibrium and local Hopf bifurcations From [26] we know that this system has seven equilibriums with certain conditions. In this section, we shall study the stability of the interior equilibrium and the existence of local Hopf bifurcations. Let E ðP 1 ; P 2 ; Z  Þ denote the interior equilibrium, where dðq1 r2  a1 q2 kÞ þ kðq2 r1  q1 r2 Þðc2  h2 Þ ; ðq1 r2  a1 q2 kÞðc1  h1 Þ þ ðq2 r1  q1 r2 Þðc2  h2 Þ dðq2 r1  a2 q1 kÞ  kðq2 r1  q1 r2 Þðc1  h1 Þ P 2 ¼ ; ðq1 r2  a1 q2 kÞðc1  h1 Þ þ ðq2 r1  q1 r2 Þðc2  h2 Þ P 1 ¼

and r1 k  r1 P 1  a2 P 2 : q1 k From Section 1, we know that c1 > h1,c2 > h2. And hence we get that the positive interior equilibrium E* exists if and only if either  q r2 ðq r1  q1 r2 kÞðc2  h2 Þ ðq2 r1  a2 q1 kÞðc2  h2 Þ þ q1 r2 ðc1  h1 Þ a1 < min 1 þ 2 ; ; q2 k dq2 q2 ðc1  h1 Þ  ½r1 r2 q1 þ a2 kðq2 r1  q1 r2 Þðc1  h1 Þ ðr2 q1  a2 q1 kÞðc2  h2 Þ ½r1 r2 q1 þ a2 kðq2 r1  a2 q1 kÞd   r1 q2 kðc1  h1 Þ  r1 q2 d q2 kðc1  h1 Þ  q2 d r1 q2 k 2 ðc1  h1 Þ  r1 q2 kd Z ¼

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

1397

and a2 <

q2 r1 ðq2 r1  q1 r2 Þðc1  h1 Þ  ; q1 k dq1

or  a1 > max

q1 r2 ðq2 r1  q1 r2 kÞðc2  h2 Þ ðq2 r1  a2 q1 kÞðc2  h2 Þ þ q1 r2 ðc1  h1 Þ ; þ ; q2 k dq2 q2 ðc1  h1 Þ  ½r1 r2 q1 þ a2 kðq2 r1  q1 r2 Þðc1  h1 Þ ðr2 q1  a2 q1 kÞðc2  h2 Þ ½r1 r2 q1 þ a2 kðq2 r1  a2 q1 kÞd   r1 q2 kðc1  h1 Þ  r1 q2 d q2 kðc1  h1 Þ  q2 d r1 q2 k 2 ðc1  h1 Þ  r1 q2 kd

and a2 >

q2 r1 ðq2 r1  q1 r2 Þðc1  h1 Þ :  q1 k dq1

Let u1 ¼ P 1  P 1 ; u2 ¼ P 2  P 2 and u3 ¼ Z  Z  . Then system (1) can be written in the following form du1 r1 P 1 r1 ¼ u1  a1 P 1 u2  q1 P 1 u3  u21  a1 u1 u2  q1 u1 u3 ; dt k k du2 r2 P 2 r2 2   ¼ a2 P 2 u1  u2  q2 P 2 u3  u2  a2 u1 u2  q2 u2 u3 ; dt k k du3  ¼ ðc1 u1 þ c2 u2 ÞZ  ½h1 u1 ðt  s1 Þ þ h2 u2 ðt  s2 ÞZ  þ ðc1 u1 þ c2 u2 Þu3  ½h1 u1 ðt  s1 Þ þ h2 u2 ðt  s2 Þu3 : dt

ð2Þ

The characteristic equation of system (2) at the equilibrium E(0, 0, 0) is k3 þ a2 k2 þ a1 k þ a0 þ ðb1 k þ b0 Þeks1 þ ðc1 k þ c0 Þeks2 ¼ 0;

ð3Þ

where a2 ¼

r1 P 1 þr2 P 2 k

;

   a1 a2 P 1 P 2 þ q1 c1 P 1 þ q2 r2 P 2 ;

   a0 ¼ r2kq1  a1 q2 c1 þ r1kq2  a2 q1 c2 P 1 P 2 Z  ;

a1 ¼



r1 r2 k2

b1 ¼ q1 h1 P 1 Z  ;  b0 ¼ h1 r2kq1  a1 q2 P 1 P 2 Z  ; c1 ¼ q2 h2 P 2 Z  ;  c0 ¼ h2 r1kq2  a2 q1 P 1 P 2 Z  : Now we use the method due to Ruan and Wei [21,22] to investigate the distribution of roots of Eq. (3). For s1 = 0 and s2 = 0, Eq. (3) becomes k3 þ a2 k2 þ ða1 þ b1 þ c1 Þk þ a0 þ b0 þ c0 ¼ 0:

ð4Þ

By Routh-Hurwitz criterion we know that if ðH 1 Þ ða0 þ b0 þ c0 Þ > 0;

a2 ða1 þ b1 þ c1 Þ > a0 þ b0 þ c0 ;

then all roots of Eq. (4) have negative real parts. Obviously, im (s1) (m > 0) is a root of Eq. (3) with s2 = 0 if and only if m3 i  a2 m2 þ ða1 þ c1 Þmi þ a0 þ c0 þ ðib1 m þ b0 Þðcos ms1  i sin ms1 Þ ¼ 0: Separating the real and imaginary parts gives  3 m þ ða1 þ c1 Þm ¼ b0 sin ms1  b1 m cos ms1 ; a2 m  ða0 þ c0 Þ ¼ b0 cos ms1 þ b1 m sin ms1 ;

ð5Þ

which leads to m6 þ pm4 þ qm2 þ r ¼ 0;

ð6Þ

1398

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

where p ¼ a22  2ða1 þ c1 Þ; q ¼ ða1 þ c1 Þ2  2a2 ða0 þ c0 Þ  b21 ; r ¼ ða0 þ c0 Þ2  b20 : Let x = m2, then Eq. (6) becomes x3 þ px2 þ qx þ r ¼ 0:

ð7Þ

Denote hðxÞ ¼ x3 þ px2 þ qx þ r:

ð8Þ

Since limx!1h(x) = + 1, we conclude that Eq. (8) has at least one positive root if r < 0. From Eq. (8) it follows that dhðxÞ ¼ 3x2 þ 2px þ q: dx Obviously, if D = p2  3q 6 0, dh(x)/dx P 0, then h(x) is monotonically increasing in x 2 [0,1). Thus for r > 0 and D 6 0, Eq. (8) has no positive roots for x 2 [0,1). When r P 0 and D > 0, the following equation 3x2 þ 2px þ q ¼ 0 has two real roots pffiffiffiffi pffiffiffiffi p þ D p  D   and x2 ¼ : x1 ¼ 3 3 p ffiffiffi ffi p ffiffiffiffi Clearly, h00 ðx1 Þ ¼ 2 D > 0 and h00 ðx2 Þ ¼ 2 D < 0, so we know that x1 and x2 are the local minimum and the local maximum of h(x), respectively. We obtain the followings. Lemma 2.1. For the polynomial Eq. (7), we have the following results. (i) If r < 0, then Eq. (7) has at least one positive root. (ii) If r P 0 and D 6 0, then Eq. (7) has no positive roots. (iii) If r P 0 and D > 0, then Eq. (7) has positive roots if and only if x1 > 0 and hðx1 Þ 6 0.

Proof. Summarizing the discussions above, we only need to prove (iii). Noticing that r P 0, x1 is the local minimum of h(x) and limx!+1h(x) = + 1, we know that the sufficiency is true. In what follows, we need to prove the necessity of iii. Suppose that either x1 6 0 or x1 > 0 and hðx1 Þ > 0. Since h(x) is increasing for x P 0 and h(0) = r P 0, we know that h(x) has no positive real roots for x1 6 0. If x1 > 0 and hðx1 Þ > 0, since hðx1 Þ is the local minimum value, it follows that h(0) = r > 0 that when x1 > 0 and hðx1 Þ > 0, h(x) has no positive real roots, too. This proves the lemma. Without loss of generality, we assume that Eq. (7) has three positive roots, denoted by x1, x2 and x3, respectively. pffiffiffiffiffi Then Eq. (5) has three positive roots mk ¼ xk ; k ¼ 1; 2; 3. By (4), we have cos mk s1 ¼ Denote ðjÞ s1k

ða2 m2k  a0  c0 Þb0 þ ½m3k  ða1 þ c1 Þmk b1 mk : b20 þ b21 m2k

" # 1 ða2 m2k  a0  c0 Þb0 þ ½m3k  ða1 þ c1 Þmk b1 mk ¼ arccos þ 2jp ; mk b20 þ b21 m2k ðjÞ

where k = 1, 2, 3; j = 0, 1, . . . Then ±imk is a pair of purely imaginary roots of Eq. (3) with s1 ¼ s1k ; s2 ¼ 0. For convenience, we state a result due to Ruan and Wei [22] as a lemma in the followings. Lemma 2.2. Consider the exponential polynomial

h i ð0Þ ð0Þ ð1Þ n1 ð1Þ eks1 þ    þ pn1 k þ pð1Þ pðk; eks1 ; . . . ; eksm Þ ¼ kn þ p1 kn1 þ    þ pn1 k þ pð0Þ n þ p1 k n h i ðmÞ ðmÞ eksm ; þ    þ p1 kn1 þ    þ pn1 k þ pðmÞ n

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

1399

ðiÞ

where si P 0 ði ¼ 1; 2; . . . ; mÞ and pj ði ¼ 0; 1; . . . ; m; j ¼ 1; 2; . . . ; nÞ are constants. As ðs1 ; s2 ; . . . ; sm Þ vary, the sum of the zeros of pðk; eks1 ; . . . ; eksm Þ on the open right half plane can change only if a zero appears on or crosses the imaginary axis. Let k(s1) = a(s1) + im(s1) be the root of Eq. (3) satisfying aðsj1k Þ ¼ 0, mðsj1k Þ ¼ mk . Then the following transversality condition holds. ðjÞ

Lemma 2.3. Suppose that xk ¼ m2k . Then the sign of a0ðs1k Þ is coincident with that of h 0 (xk). Proof. Let k = k(s1) be the root of Eq. (3). Substituting k(s1) into Eq. (3) and differentiating the both sides of Eq. (3) with respect to s1, it follows that f3k2 þ 2a2 k þ a1 þ c1 þ ½b1  s1 ðb1 k þ b0 Þeks1 g

dk ¼ kðb1 k þ b0 Þeks1 : ds1

Thus, 

dk ds1

1 ¼

ð3k2 þ 2a2 k þ a1 þ c1 Þeks1 b1 s1 þ  kðb1 k þ b0 Þ kðb1 k þ b0 Þ k

From (5)–(7), we have

2 

 ð3k þ 2a2 k þ a1 þ cÞeks1 b1 0 ðjÞ 1 þ Re a ðsk Þ ¼ Re kðb1 k þ b0 Þ kðb1 k þ b0 Þ n o 1 1 xk ¼ 3m6k þ 2½a22  2ða1 þ c1 Þm4k þ ½ða1 þ c1 Þ2  2a2 ða0 þ c0 Þ  b21 m2k ¼ ð3m6k þ 2pm4k þ qm2k Þ ¼ h0 ðxk Þ; K K K where K ¼ b21 m4k þ b20 m2k . Notice that K > 0 and xk > 0, we conclude that   sign a0 ðsj1k Þ ¼ signfh0 ðxk Þg This proves the lemma. From Lemmas 2.1, 2.2 and 2.3 we know that if Eq. (6) has more than one positive roots. Then the stability switch may exist. Summarizing the above discussions we can ensure the stability interval. Let I be the stability interval. We have the followings. Theorem 2.1. Suppose that (H1) is satisfied and s2 = 0. (i) when r P 0 and D 6 0, all roots of Eq. (3) have negative real parts for all s1 P 0, and hence the equilibrium E of system (2) is asymptotically stable for all s1 P 0. (ii) If either r < 0 or r P 0, D > 0, x1 > 0 and hðx1 Þ 6 0 hold, then h(x) has at least a positive root xk, the stability switch may exist, all roots of Eq. (3) have negative real parts for s1 2 I. Hence the equilibrium E of system (2) is asymptotically stable for s1 2 I. (iii) If all conditions as stated in ii and h 0 (xk)50 hold, then system (2) undergoes a Hopf bifurcation at the equilibriumE, ðjÞ when s1 ¼ s1k ðj ¼ 0; 1; 2; . . .Þ. Now let s1 ¼ s1 2 I, s2 > 0, and k = ix(s2)(x > 0) be a root of Eq. (3). Then we get ( x3 þ a1 x þ b1 x cos xs1  b0 sin xs1 ¼ c0 sin xs1  c1 x cos xs2 ; a2 x2  a0  b0 cos xs1  b1 sin xs1 ¼ c1 x sin xs2 þ c0 cos xs1 :

ð9Þ

From Eq. (9) it follows that  x6 þ ða22  2a1 Þx4 þ a21  2a0 a2 þ b21  c21 x2 þ a20 þ b20 þ c20 þ 2½b0 ðx3  a1 xÞ  b1 xða2 x2  a0 Þ sin xs1 þ 2½b1 xða1 x  x3 Þ  b0 ða2 x2  a0 Þ cos xs1 ¼ 0:

ð10Þ

We know that Eq. (10) has at most finite positive roots. If Eq. (10) has positive root, without loss of generality, we assume that Eq. (10) has N positive roots, denoted by xi ; ði ¼ 1; 2; . . . ; N Þ. Notice Eq. (9) we get ðjÞ s2i ; i ¼ 1; 2; . . . ; N; j ¼ 0; 1; 2; . . .. Define

1400

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409 ð0Þ

s02 ¼ s2i0 ¼

n

min

i2f1;2;...;N g

ð0Þ

s2i

o ; x0 ¼ xi0 : ðjÞ

ðjÞ

Let k(s2) = a(s2) + ix(s2) be the root of Eq. (3) satisfying aðs2i Þ ¼ 0, xðs2i Þ ¼ xi . By computation, we get  a0 ðs02 Þ ¼ K1 ðb1 c1 s1 x30  b0 c0 s1 x0 þ b1 c0 x0 Þ sin x0 ðs1 þ s02 Þ þ ðb0 c1 s1 þ b1 c0 s1  b1 c1 Þx20 cos x0 ðs1 þ s02 Þ



 1 þ c0 ða1 x0 þ 3x30 Þ þ 2a2 b1 x30 sin x0 s02  c1 x20 þ 2a2 c0 x20  b1 x0 ða1 x0 þ 3x30 Þ cos x0 s02 ; where K1 ¼

c21 x40

þ

c20 x20 .

ð11Þ

Summarizing the discussions above, we have the following conclusions.

Theorem 2.2. Suppose that (H1), s1 ¼ s1 2 I hold and Eq. (10) has positive roots. s02 and a0 ðs02 Þ have the same meaning as last definition. We have (i) All roots of Eq. (3) have negative real parts for s2 2 ½0; s02 Þ, and the equilibrium E of system (2) is asymptotically stable for s2 2 ½0; s02 Þ. (ii) If a0 ðs02 Þ–0 hold, then system (2) undergoes a Hopf bifurcation at the equilibrium E, when s2 ¼ s02 .

3. Stability and direction of the Hopf bifurcation In the previous section, we obtained conditions for Hopf bifurcation to occur when s2 ¼ s02 . In this section we study the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions when s2 ¼ s02 , using techniques from normal form and center manifold theory (see e.g. Hassard et al. [13]). We assume s1 > s02 . Letting s2 ¼ s02 þ l and dropping the bars for simplification of notations, system (2) is transformed into an FDE in C ¼ Cð½s1 ; 0; R3 Þ as _ ¼ Ll ðut Þ þ F ðl; ut Þ; uðtÞ

ð12Þ 3

3

where ut(h) = u(t + h) 2 C, and Ll:C ! R , F:R · C ! R are given, respectively, by  Ll u ¼ A1 uð0Þ þ B1 u s1 þ B2 uðs02 Þ;

ð13Þ

where 0



r1 P 1 k

B B A1 ¼ B a2 P  1 @  c1 Z 0 0 B B B1 ¼ @ 0 h1 Z 

a1 P 1 

r2 P 2 k

c2 Z  1 0 0 C 0 0C A;

q1 P 1

1

C C ; q2 P 2 C A 0 0

1

0

0

0

B B2 ¼ B @0

0

C 0C A;

0 0

0

h2 Z 

0

T

uðtÞ ¼ ðu1 ðtÞ; u2 ðtÞ; u3 ðtÞÞ ; and

0

 rk1 u21 ð0Þ  a1 u1 ð0Þu2 ð0Þ  q1 u1 ð0Þu3 ð0Þ

B r2 2 F ðl; uÞ ¼ B @ a2 u1 ð0Þu2 ð0Þ  k u2 ð0Þ  q2 u2 ð0Þu3 ð0Þ

1 C C: A

ð14Þ

½c1 u1 ð0Þ þ c2 u2 ð0Þ  h1 u1 ðs1 Þ  h2 u2 ðs02 Þu3 ð0Þ From the discussion in Section 2, we know that system (10) undergoes a Hopf bifurcation at (0,0,0) when l = 0, and the associated characteristic equation of system (10) with l = 0 has a pair of simple imaginary roots ±ix0. By the Riesz representation theorem, there exists a function g(h,l) of bounded variation for h 2 ½s1 ; 0, such that Z 0 Ll u ¼ dgðh; lÞuðhÞ 8u 2 C: ð15Þ s1

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

In fact, we can choose 8 h ¼ 0; > < A1 þ B2 ; h 2 ½s02 ; 0Þ; gðh; lÞ ¼ B2 ; > :  B1 dðh þ s1 Þ; h 2 ½s1 ; s02 Þ;  0; h–0; where dðhÞ ¼ 1; h ¼ 0:

1401

ð16Þ

For u 2 C 1 ð½s1 ; 0; R3 Þ, we set 8 < duðhÞ ; h 2 ½s1 ; 0Þ; dh AðlÞu ¼ R 0 : s dgðn; lÞuðnÞ; h ¼ 0;

ð17Þ

1

and

 RðlÞu ¼

h 2 ½s1 ; 0Þ; h ¼ 0:

0; F ðl; uÞ;

ð18Þ

Then system (10) can be rewritten as u_ t ¼ AðlÞut þ RðlÞut ;

ð19Þ ½s1 ; 0.

For w 2 C where ut(h) = u(t + h) for h 2 8 duðsÞ < ; s 2 ð0; s1 ; ds A u ¼ R 0 : s dgT ðt; 0ÞuðtÞ; s ¼ 0:

1

([0,1],(R3)*),

define

1

and w 2 C 1 ð½0; s1 ; ðR3 Þ Þ, using the bilinear form Z 0 Z h  T ðn  hÞdgðhÞuðnÞdn;  T ð0Þuð0Þ  w hw; ui ¼ w

For u 2 C

1

ð½s1 ; 0; R3 Þ

s1

ð20Þ

n¼0

we know that A* and A = A(0) are adjoint operators. By the discussion in Section 2, we know that ±ix0 are eigenvalues of A(0). Thus they are eigenvalues of A*. Suppose that qðhÞ ¼ qð0Þeix0 h is an eigenvector of A(0) corresponding to the eigenvalue ix0. Then, A(0) = ix0q(h). When h = 0, we obtain " # Z 0 ix0 h qð0Þ ¼ 0; ð21Þ dgðhÞe ix0 I  s1

which yields q(0) = (1, a, b)T, where a¼

ðix0 k þ r1 P 1 Þq2 P 2  a2 q1 kP 1 P 2 ; ðix0 k þ r2 P 2 Þq1 P 1  a1 q2 kP 1 P 2

b¼

a1 a2 k 2 P 1 P 2  ðix0 k þ r1 P 1 Þðix0 k þ r2 P 2 Þ : ðix0 k þ r2 P 2 Þq1 kP 1  a1 q2 k 2 P 1 P 2

Similarly, it can be verified that q ðsÞ ¼ Dð1; a ; b Þeix0 s is the eigenvector of A* corresponding to ix0, where a ¼

ðr1 P 1  ix0 kÞq2 P 2  a2 q1 kP 1 P 2 ; ðr2 P 2  ix0 kÞq1 P 1  a1 q2 kP 1 P 2

b ¼

a1 a2 k 2 P 1 P 2  ðr1 P 1  ix0 kÞðr2 P 2  ix0 kÞ : ðix0 k þ r2 P 2 Þq1 kP 1  a1 q2 k 2 P 1 P 2

By (20), we get  hq ðsÞ; qðhÞi ¼ D 1; a ; b ð1; a; bÞT  ¼ 1 þ aa þ bb 

Z

Z

0

s1 0



Z

   ix ðnhÞ 1;  a ; b e 0 dgðhÞð1; a; bÞT eix0 n dn

h

n¼0

1; a ; b heix0 h dgðhÞð1; a; bÞT

s1 



¼ 1 þ aa þ bb þ

 s1 h1 b Z  eix0 s1

0

þ s02 h1 a b Z  eix0 s2 :

1402

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

We can choose ¼ D

1 

0

1 þ aa þ bb þ s1 h1 b Z  eix0 s1 þ s02 h1 ab Z  eix0 s2

;

so that hq*,qi = 1. Following the algorithms in [13] and using the same notations as there to compute the coordinates describing the center manifold C0 at l = 0. Let ut be the solution of Eq. (10) when l = 0. Define zðtÞ ¼ hq ðsÞ; ut ðhÞi;

W ðt; hÞ ¼ ut ðhÞ  2RezðtÞqðhÞ:

ð22Þ

On the center manifold C0 we have W ðt; hÞ ¼ W ðzðtÞ; zðtÞ; hÞ; where W ðz; z; hÞ ¼ W 20 ðhÞ

z2 z2 þ W 11 ðhÞzz þ W 02 ðhÞ þ    ; 2 2

 . Note that W is real if ut is real. We z and z are local coordinates for center manifold C0 in the direction of q* and q consider only real solution. For solution ut 2 C0 of (13), since l = 0, z_ ðtÞ ¼ ix0 z þ hq ðhÞ; F ð0; W ðz; z; hÞ þ 2RefzðtÞqðhÞgÞi q ð0ÞF 0 ðz; zÞ: ¼ ix0 z þ q ð0ÞF ð0; W ðz; z; 0Þ þ 2RefzðtÞqðhÞgÞ,ix0 z þ  We rewrite this equation as z_ ðtÞ ¼ ix0 zðtÞ þ gðz; zÞ;

ð23Þ

where gðz; zÞ ¼ q ð0ÞF 0 ðz; zÞ ¼ g20

z2 z2 z2z þ g11 zz þ g02 þ g21 þ  2 2 2

By (20), we have ut ðhÞ ¼ ðu1t ; u2t ; u3t Þ ¼ W ðt; hÞ þ zqðhÞ þ zqðhÞ; and then z2 z2 ð1Þ ð1Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ    ; 2 2 z2 z2 ð1Þ ð1Þ ð1Þ u1t ð0Þ ¼ z þ z þ W 20 ð0Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ    ; 2 2 2 z z2 ð3Þ ð3Þ ð3Þ u3t ð0Þ ¼ az þ az þ W 20 ð0Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ    ; 2 2 2 z z2   ð1Þ ð1Þ ð1Þ u1t ðs1 Þ ¼ aeix0 s1 z þ aeix0 s1 z þ W 20 ðs1 Þ þ W 11 ðs1 Þzz þ W 02 ðs1 Þ þ    ; 2 2 2 z z2 0 0 ð1Þ ð1Þ ð1Þ u2t ðs02 Þ ¼ aeix0 s2 z þ aeix0 s2 z þ W 20 ðs02 Þ þ W 11 ðs02 Þzz þ W 02 ðs02 Þ þ    2 2 ð1Þ

u1t ð0Þ ¼ z þ z þ W 20 ð0Þ

Substitute ut(h) into (12), and comparing the coefficients with (22), we have h r   i r2  0 1 b ðc1 b þ c2 ab  h1 beix0 s1  h2 abeix0 s2 Þ ; g20 ¼ 2D   a1 a  q1 b  a a2 a þ a2 þ q2 ab þ  k  k 2r1 2r2 aÞ þ bþ abÞþ b ½c1 ðb þ  bÞ  a1 ða þ aÞ  q1 ðb þ bÞ  a ½a2 ða þ  a a þ q2 ða g11 ¼ D  k k   0 0 beix0 s2 Þg; þc2 ðab þ abÞ  h1 ðbeix0 s1 þ beix0 s1 Þ  h2 ðabeix0 s2 þ a   h r i r2  0 1 b þ b ðc1  a b  h1  a beix0 s2 Þ ; b þ c2  beix0 s1  h2  g20 ¼ 2D   a1 a  q1 b  a a2 a þ a2 þ q2 a k k

ð24Þ

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

1403

n r 1 ð1Þ ð1Þ ð2Þ ð2Þ ð1Þ ð1Þ ð3Þ g21 ¼ D  ð4W 11 ð0Þ þ 2W 20 ð0ÞÞ  a1 ð2W 11 ð0Þ þ W 20 ð0Þ þ 2aW 11 ð0Þ þ  aW 20 ð0ÞÞ  q1 ð2W 11 ð0Þ k hr 2 ð3Þ ð1Þ ð1Þ ð2Þ ð2Þ ð2Þ þ W 20 ð0Þ þ 2bW 11 ð0Þ þ bW 20 ð0ÞÞ  a ð4W 11 ð0Þ þ 2W 20 ð0ÞÞ þ a2 ð2W 11 ð0Þ k i ð2Þ ð1Þ ð1Þ ð3Þ ð3Þ ð2Þ ð2Þ aW 20 ð0Þ þ 2bW 11 ð0Þ þ  bW 20 ð0ÞÞ þW 20 ð0Þ þ 2aW 11 ð0Þ þ aW 20 ð0ÞÞþq2 ð2aW 11 ð0Þ þ  ð3Þ ð3Þ ð1Þ ð3Þ ð2Þ ð2Þ  ð1Þ ð0ÞÞ þ c2 ðð2aW ð3Þ ð0Þ þ  þ b ½c1 ð2W 11 ð0Þ þ W 20 ð0Þ þ 2bW 11 ð0Þ þ bW aW 20 ð0Þ þ 2bW 11 ð0Þ þ  bW 20 ð0ÞÞÞ 20 11       s s    ð3Þ ð3Þ ð1Þ ð1Þ ð3Þ  h1 2eix0 s1 W 11 ð0Þ þ eix0 s1 W 20 ð0Þ þ 2bW 11  10 þ  bW 20  10 h2 ð2aeix0 s1 W 11 ð0Þ s2 s2 o  ð3Þ ð2Þ ð2Þ þaeix0 s1 W 20 ð0Þ þ 2bW 11 ð1Þ þ bW 20 ð1ÞÞ :

We still need to compute W20(h) and W11(h). From (19) and (22), we have  AW  2Refq ð0ÞF 0 qðhÞg; h 2 ½s1 ; 0Þ _ ,AW þ H ðz; z; hÞ: W ¼ u_ t  z_ q  z_ q ¼ AW  2Refq ð0ÞF 0 qð0Þg þ F 0 ; h ¼ 0

ð25Þ

Here H ðz; z; hÞ ¼ H 20 ðhÞ

z2 z2 þ H 11 ðhÞzz þ H 02 ðhÞ þ    2 2

ð26Þ

Expanding the above series and comparing the coefficients, we obtain ½A  2ix0 IW 20 ðhÞ ¼ H 20 ðhÞ; AW 11 ðhÞ ¼ H 11 ðhÞ; . . . :

ð27Þ

By (23), we know that for h 2 ½s1 ; 0Þ, qðhÞ: H ðz; z; hÞ ¼ q ð0ÞF 0 qðhÞ  q ð0ÞF 0 qðhÞ ¼ gqðhÞ  g

ð28Þ

Comparing the coefficients with (26) gives that H 20 ðhÞ ¼ g20 qðhÞ  g02 qðhÞ; H 11 ðhÞ ¼ g11 qðhÞ  g11 qðhÞ:

ð29Þ

From (25), (27) and the definition of A, we obtain W 20 ðhÞ ¼ 2ix0 W 20 ðhÞ þ g20 qðhÞ þ g20 qðhÞ:

ð30Þ

Solving for W20(h), we obtain W 20 ðhÞ ¼

ig20 i g02 qð0Þeix0 h þ E1 e2ix0 h ; qð0Þeix0 h þ x0 3x0

ð31Þ

ig11 i g11 qð0Þeix0 h þ E2 ; qð0Þeix0 h þ x0 x0

ð32Þ

and similarly W 11 ðhÞ ¼ 

where E1 and E2 are both three-dimensional vectors, and can be determined by setting h = 0 in H. In fact, since   H ðz; z; 0Þ ¼ 2Re q ð0ÞF 0 qð0Þ þ F 0 ; we have H 20 ð0Þ ¼ g20 qð0Þ  g02 qð0Þ þ F z2 ;

ð33Þ

and H 11 ð0Þ ¼ g11 qð0Þ  g11 qð0Þ þ F zz ; where F 0 ¼ F z2

z2 z2 þ F zz zz þ F z2 þ    2 2

ð34Þ

1404

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

Hence combining the definition of A, we get Z 0 dgðhÞW 20 ðhÞ ¼ 2ix0 W 20 ðhÞ þ g20 qð0Þ þ g02 qð0Þ  F z2 ;

ð35Þ

s1

and Z

0 s1

dgðhÞW 11 ðhÞ ¼ g11 qð0Þ þ g11 qð0Þ  F zz :

Notice that "

Z

ix0 I 

#

0 ix0 h

e

"

dgðhÞ qð0Þ ¼ 0;

ix0 I 

ð36Þ

Z

s1

we have " 2ix0 I 

Z

#

0 ix0 h

e

dgðhÞ  qð0Þ ¼ 0;

s1

#

0 2ix0 h

dgðhÞ E1 ¼ F z2 :

e s1

Similarly, we have "Z # 0



dgðhÞ E2 ¼ F zz :

s1

Hence, we get 0

2ix0 þ

B B @

r1 P 1 k

a1 P 1

a2 P 2

2ix0 þ 

c1 Z  þ h1 Z  e2ix0 s1

q1 P 1 r2 P 2 k  2ix0 s02

c2 Z  þ h2 Z e

1

0

C B q2 P 2 C AE1 ¼ 2@ 2ix0

1

 rk1  a1 a  q1 b a2 a  rk2 a2  q2 ab 

0

C A

c1 b þ c2 ab  h1 beix0 s1  h2 abeix0 s2

and 0 B B @

r1 P 1 k

a1 P 1

a2 P 2

r2 P 2 k

ðh1  c1 ÞZ  ðh2  c2 ÞZ 

1

0 1 aÞ  q1 ðb þ  bÞ  2rk1  a1 ða þ  C B C C a2 ða þ  aÞ  2rk2 a a  q2 ða bþ abÞ A: q2 P 2 AE2 ¼ @ 0 0   ix0 s2 ix0 s2 ix0 s1 ix0 s1     c1 ðb þ bÞ þ c2 ðab þ  abÞ  h1 ðbe þ be Þ  h2 ð abe þ abe Þ 0 q1 P 1

Then g21 can be expressed by the parameters. Based on the above analysis, we can see that each gij can be determined by the parameters. Thus we can compute the following quantities: i C 1 ð0Þ ¼ 2x0 l2 ¼ 

jg j2 g11 g20  2jg11 j  02 3 2

ReðC 1 ð0ÞÞ ; Reðk0 ðs02 ÞÞ

! þ

g21 ; 2 ð37Þ

b2 ¼ 2ReðC 1 ð0ÞÞ; T2 ¼ 

ImðC 1 ð0ÞÞ þ l2 Imðk0 ðs02 ÞÞ : x0

Theorem 3.1. l2 determines the direction of the Hopf bifurcation: if l2 > 0(l2 < 0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for s2 > s02 ðs2 < s02 Þ; b2 determines the stability of bifurcating periodic solutions: the bifurcating periodic solutions are orbitally asymptotically stable (unstable) if s2 > s02 ðs2 < s02 Þ; and T2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2 > 0(T2 < 0).

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409 1

t-P1 plane

P1 0.5 0

1405

1

0

200

400

600

800

1

1000

t-P2 plane

P2 0.5 P2 0.5 0

0

200

400

600

800

1000

t-Z plane

6

0 1 7

Z 4 2

P1 0

200

400

600

800

6

0.5

5 4

Z

3

1000

0

t

2

Fig. 1. Behavior and phase portrait of system (1) with s1 = 6.0, s2 = 0. The interior equilibrium is stable. The initial value is (1, 1, 4).

1

t-P1 plane

P1 0.5 0

1 0

200

400

600

1

t-P2 plane

P2 0.5

P2 0.5 0

0

200

400

600

t-Z plane

6

0 1 7

Z 4 2

P1 0

200

400

6

0.5

5 4

Z

3

600

2

0

t

Fig. 2. Behavior and phase portrait of system (1) with s1 = 8.0, s2 = 0. Hopf bifurcation occurs from the interior equilibrium. The initial value is (1, 1, 4). 1

t-P1 plane

P2

P1 0.5 0

1 0

200

400

600

800

1

1000

t-P2 plane

P2 0.5

P2 0.5 0

0

200

400

600

800

1000

t-Z plane

6

0 1 7

Z 4 2

P1 0

200

400

600

t

800

1000

P1

6

0.5

5 4 0

3

Z

2

Fig. 3. Behavior and phase portrait of system (1) with s1 = 15.0, s2 = 0. The interior equilibrium become stable again. The initial value is (1, 1, 4).

1406

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409 1

t-P1 plane

P1 0.5 0

1 0

200

400

600

1

t-P2 plane

P2 0.5

P2 0.5 0

0

200

400

600

t-Z plane

6

Z

0 1 7

4

6

P1 0.5 2

0

200

400

5 4

Z

3

600

2

0

t

Fig. 4. Behavior and phase portrait of system (1) with s1 = 20.0, s2 = 0. Hopf bifurcation occurs again from the interior equilibrium. The initial value is (1, 1, 4). 1

t-P1 plane

P1 0.5 0

1 0

200

400

600

800

1

1000

t-P2 plane

P2 0.5

P2 0.5 0

0

200

400

600

800

6

Z

1000

t-Z plane

0 1 7

4

6

P1 0.5 2

0

200

400

600

800

5 4 3

1000

0

t

Z

2

Fig. 5. Behavior and phase portrait of system (1) with s1 = 6.0, s2 = 5.0. The interior equilibrium is stable. The initial value is (1, 1, 4). 1

t-P1 plane

P1 0.5 0

1 0

200

400

1

600

t-P2 plane

P2 0.5

P2 0.5 0

0

200

400

600

t-Z plane

6

0 1 7

Z 4

P1 2

0

200

400

t

600

6

0.5

5 4 3 0

Z

2

Fig. 6. Behavior and phase portrait of system (1) with s1 = 6.0, s2 = 5.9. Hopf bifurcation occurs from the interior equilibrium. The initial value is (1, 1, 4).

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409 1

1407

t-P1 plane

P1 0.5 1 0

0

200

400

600

800

1

1000

t-P2 plane

P2 0.5

P2 0.5 0

0

200

400

600

800

1000

t-Z plane

6

0 1 7

Z 4 2

6

P1 0.5 0

200

400

600

800

5 4 3

1000

0

t

Z

2

Fig. 7. Behavior and phase portrait of system (1) with s1 = 15.0, s2 = 6.0. The interior equilibrium is stable. The initial value is (1, 1, 4).

2

t-P1 plane

P1 1 2 0

0

200

400

600

t-P2 plane

2

P2 1

P2 1 0

0

200

400

10

600

t-Z plane

0 2

Z 5

10 8

P1 1 0

0

200

400

t

600

6 4 0

2

Z

0

Fig. 8. Behavior and phase portrait of system (1) with s1 = 15.0, s2 = 10.0. Hopf bifurcation occurs from the interior equilibrium. The initial value is (1, 1, 4).

4. Numerical examples In this section, we shall use MATLAB to perform some numerical simulations on system (1). We take the same parameter values as the literature [26]: r1 = 2.5, r2 = 2.3, a1 = 0.01, a2 = 0.02, q1 = 0.6, q2 = 0.55, c1 = 0.4, c2 = 0.3, h1 = 0.06, h2 = 0.07, k = 20, d = 0.1. For this set of parameter values we observe that the positive interior equilibrium is E*(0.1501, 0.2129, 4.1283). We will now try to find the effect of single and multiple delays on the dynamics of the system (1). With these parameters, we find that (H1) hold. When s2 = 0, by direct computation we get Eq. (6) has two ðjÞ : positive roots m1 G 0.5278 and m2 G 0.4190. Substituting these parameters into Eq. (5) we get s11 ¼7:0371 þ 11:9045j ðjÞ : ðjÞ ðjÞ and s12 ¼11:6274 þ 12:7969jðj ¼ 0; 1; . . .Þ. Eq. (3) has pure imaginary roots when s1 ¼ s11 or s1 ¼ s12 . Further 0 ðjÞ 0 ðjÞ a ðs11 Þ > 0 and a ðs12 Þ < 0. By Theorem 2.1 we know that the stability switches exist and the equilibrium E* is asymptotically stable when s1 2 ½0; 7:0371Þ [ ð11:6247; 18:9416Þ [    [ ð114:0010; 114:1776Þ. The results are illustrated by Figs. 1–4. : First, let s1 = 6.0 2 [0, 7.0371), we get s02 ¼5:5695. By Theorem 2.2 we know that the interior equilibrium E* is asymptotically stable for s2 2 [0, 5.5695). Furthermore by (37) and direct computation, we have C1(0) G 0.9734  1.0325i, : b2 G 1.9468 < 0, and l2 G 13.5851 > 0. By Theorem 3.1 we know that, at s02 ¼5:5695, the bifurcating periodic solution is orbitally asymptotically stable, and the direction of the Hopf bifurcation is supercritical. The results are illustrated by : Figs. 5 and 6. Second, let s1 = 15.0 2 (11.6274, 18.9416). Similarly, we have s02 ¼8:8025, C1(0) G 1.4622  0.0365i,

1408

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

b2 G 2.9244 < 0, l2 G 22.9389 > 0. By Theorems 2.2 and 3.1 we know that the interior equilibrium E* is asymptot: ically stable for s2 2 [0, 8.8025), and the direction of the Hopf bifurcation at s02 ¼8:8025 is supercritical, and the bifurcating periodic solution is orbitally asymptotically stable. The results are illustrated by Figs. 7 and 8.

5. Conclusion The factors effecting plankton dynamics and the abundance of harmful phytoplankton species have been the subject of intensive research recently. As mentioned in the Introduction, researchers have established with the help of experimental results and mathematical modeling, that toxin-producing phytoplankton may be used as a controlling agent for the termination of planktonic bloom. But few of those studies contain the presence of two harmful phytoplankton in such situation. Moreover, the effects of time lags due to liberation of toxic substances can not be ignored in this context. The reason of occurrence of planktonic bloom and its possible control mechanism is still in infancy, hence the progress of such important areas urgently requires special attention both from experimental and mathematical ecologists. In this paper we have analyzed a model consisting of two harmful phytoplankton and zooplankton. First, in the case of s2 = 0, by analyzing the distribution of the roots of the corresponding characteristic equation, we find out that there are stability switches for the interior equilibrium when s1 varies. Then for s1 in a stability interval, regarding the delay s2 as parameter, we show that there exists a first critical value of s2 at which the interior equilibrium loses its stability and the Hopf bifurcation occurs. We also investigate the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions. From the previous studies we know that the positive impact of toxin-producing phytoplankton for the control of algal bloom. In this paper we also arrive at the same conclusion. Thus we may finally conclude that harmful phytoplankton may be used as a bio-control agent for the HAB problems. The behavior of environmental fluctuations in the two harmful phytoplankton-zooplankton dynamics may give some interesting results and needs further investigations.

Acknowledgement The authors are grateful to Professor M.S. El Naschie for his helpful comments and valuable suggestions.

References [1] Abdusalam H, Fahmy E. Cross-diffusional effect in a telegraph reaction diffusion Lotka-Volterra two competitive system. Chaos, Solitons & Fractals 2003;18:259–66. [2] Abdallah S. Stability and persistence in plankton models with distributed delays. Chaos, Solitons & Fractals 2003;17:879–84. [3] Anderson D. Toxic algae blooms and red tides: a global perspective. In: Red tides: biology, environmental science and toxicology. New York: Elsevier; 1989. p. 11–21. [4] Arzul G, Siguel M, Guzman M, Denn E. Comparison of allelopathic properties in three toxic Alexandrium species. J Exp Mar Biol Ecol 1999;232:285–95. [5] Blaxter J, Southward A. Advances in marine biology. London: Academic Press; 1997. [6] Buskey E, Montagna P, Amos A, Whitledge T. Disruption of grazer populations as a contributing factor to the initiation of the Texas brown tide algal bloom. Limnol Oceanogr 1997;42:1215–22. [7] Chattopadhyay J, Sarkar R, Abdllaoui A. A delay differential equation model on harmful algal blooms in the presence of toxic substances. IMA J Math Appl Med Biol 2002;19:137–61. [8] Chattopadhyay J, Sarkar R, Mandal S. Toxin-producing plankton may act as a biological control for planktonic blooms-field study and mathematical modelling. J Theor Biol 2002;215:333–44. [9] Edwards M, Bees M. Generic dynamics of a simple plankton population model with a non-integer exponent of closure. Chaos, Solitons & Fractals 2001;12:289–300. [10] Fistarol G, Legrand C, Graneli E. Allelopathic effect of Primnesium parvum on a natural plankton community. Mar Ecol Prog Ser 2003;255:115–25. [11] Gakkhar S, Naji R. Chaos in seasonally perturbed ratio-dependent prey–predator system. Chaos Solitons & Fractals 2003;15:107–18. [12] Hallegraeff G. A review of harmful algae blooms and the apparent global increase. Phycologia 1993;32:79–99. [13] Hassard B, Kazarinoff N, Wan Y. Theory and application of Hopf bifurcation. Cambridge: Cambridge University Press; 1981. [14] Huang D, Wang H, Feng J, Zhu Z. Hopf bifurcation of the stochastic model on HAB nonlinear stochastic dynamics. Chaos, Solitons & Fractals 2006;27:1072–9.

J. Zhao, J. Wei / Chaos, Solitons and Fractals 39 (2009) 1395–1409

1409

[15] Hui J, Chen L. Dynamic complexities in a periodically pulsed ratio-dependent predator-prey ecosystem modeled on a chemostat. Chaos, Solitons & Fractals 2006;29:407–16. [16] Jing Z, Yang J. Bifurcation and chaos in discrete-time predator–prey system. Chaos, Solitons & Fractals 2006;27:259–77. [17] Lawrence J, Chan A, Suttle C. A novel virus (HaNIV) causes lysis of the toxic bloom-forming alga Heterosigma akashiwo (Raphidophyceae). J Phycol 2001;37:216–22. [18] Li Z, Wang W, Wang H. The dynamics of a Beddington-type system with impulsive control strategy. Chaos, Solitons & Fractals 2006;29:1229–39. [19] Maneiro I, Frango´pulos M, Guisande C, Ferna´ndez M, Reguera B, Riveiro I. Zooplankton as a potential vector of diarrhetic shellfish poisoning toxins through the food web. Mar Ecol Prog Ser 2000;201:155–63. [20] Meng X, Wei J. Stability and bifurcation of mutual system with time delay. Chaos, Solitons & Fractals 2004;21:729–40. [21] Ruan S, Wei J. On the zeros of a third degree exponential polynomial with applications to a delayed model for the control of testosterone secretion. IMA J Math Appl Med Biol 2001;18:41–52. [22] Ruan S, Wei J. On the zeros of transcendental functions with applications to stability of delay differential equations with two delays. Dyn Contin Discrete Impuls Syst Ser A Math Anal 2003;10:863–74. [23] Sarkar R, Chattopadhyay J. Occurrence of planktonic blooms under environmental fluctuations and its possible control mechanism -mathematical models and experimental observations. J Theor Biol 2003;224:501–16. [24] Sarkar R, Malchow H. Nutrient and toxin producing phytoplankton control algal blooms – a spatiotemporal study in a noisy environment. J Biosci 2005;30:749–60. [25] Sarkar R, Petrovskii S, Biswas M, Gupta A, Chattopadhyay J. An ecological study of a marine plankton community based on the field data collected from Bay of Bengal. Ecol Model 2006;193:589–601. [26] Sarkar R, Mukhopadhyay B, Bhattacharyya R, Sandip Banerjee. Time lags can control algal bloom in two harmful phytoplankton–zooplankton system. Appl Math Comp 2006. doi:10.1016/j.amc.2006.07.113. [27] Smayda T. Novel and nuisance phytoplankton blooms in the sea: evidence for a global epidemic. Toxic marine phytoplankton. New York: Elsevier; 1990. [28] Sole J, Garcia-Ladona E, Ruardij P, Estrada M. Modelling allelopathy among marine algae. Ecol Model 2005;183:373–84. [29] Stoermer E, Smol J. The diatoms. Cambridge: Cambridge University Press; 1999. [30] Stoecker D, Guillard R, Kavee R. Selective predation by Favella ehrenbergii (Tintinnia) on and among dinoflagellates. Biol Bull 1981;160:36–145. [31] Tikhonov D, Enderlein J, Malchow H, Medvinsky A. Chaos and fractals in fish school motion. Chaos, Solitons & Fractals 2001;12:277–88. [32] Turner J, Tester P. Toxic marine phytoplankton, zooplankton grazers and pelagic food webs. Limnol Oceanogr 1997;42:1203–14. [33] Turner J, Doucette G, Powell C, Kulis D, Keafer B, Anderson D. Accumulation of red tide toxins in larger size fractions of zooplankton assemblages from Massachusetts Bay, USA. Mar Ecol Prog Ser 2000;203:95–107. [34] Xiang Z, Song X. Extinction and permanence of a two-prey two-predator system with impulsive on the predator. Chaos, Solitons & Fractals 2006;29:1121–36. [35] Zhang S, Tan D, Chen L. Chaos in periodically forced Holling type II predator–prey system with impulsive perturbations. Chaos, Solitons & Fractals 2006;28:367–76. [36] Zhang S, Chen L. A study of predator–prey models with the Beddington–DeAnglis functional response and impulsive effect. Chaos, Solitons & Fractals 2006;27:237–48. [37] Zhang Y, Xiu Z, Chen L. Chaos in a food chain chemostat with pulsed input and washout. Chaos, Solitons & Fractals 2005;26:159–66.