Nonlinear Analysis: Real World Applications 10 (2009) 314–332
www.elsevier.com/locate/nonrwa
Dynamical analysis of toxin producing Phytoplankton–Zooplankton
interactions
Tapan Saha, Malay Bandyopadhyay∗
Department of Mathematics, Scottish Church College, 1 & 3, Urquhart Square, Kolkata-700 006, India
Received 5 April 2007; accepted 13 September 2007
Abstract
In the present paper we consider a toxin producing phytoplankton–zooplankton model in which the toxin liberation by
phytoplankton species follows a discrete time variation. Firstly we consider the elementary dynamical properties of the toxicphytoplankton–zooplankton
interacting model system in absence of time delay. Then we establish the existence of local Hopfbifurcation
as the time delay crosses a threshold value and also prove the existence of stability switching phenomena. Explicit
results are derived for stability and direction of the bifurcating periodic orbit by using normal form theory and center manifold
arguments. Global existence of periodic orbits is also established by using a global Hopf-bifurcation theorem. Finally, the basic
outcomes are mentioned along with numerical results to provide some support to the analytical findings.
c 2007 Elsevier Ltd. All rights reserved.
Keywords: Toxic phytoplankton; Time delay; Hopf-bifurcation; Normal form; Center manifold argument; Global Hopf-bifurcation
1. Introduction
Almost all aquatic life is based upon plankton, which are the most abundant form of life floating freely near the
surfaces of all aquatic environments, namely, lakes, rivers, estuaries and oceans [1,2]. The plant forms of plankton
community are known as phytoplankton, they are capable of photo-synthesis in presence of sunlight and serve as the
basic food source and occupy the first trophic level for all aquatic food chains. The animals in the plankton community
are known as zooplankton. Phytoplankton are consumed by zooplankton which are a most favourable food source for
fish and other aquatic animals. Phytoplankton are not only the basis for all aquatic food chains, they render very
useful service by producing a huge amount of oxygen for human and other living animals after absorbing carbondioxide
from surrounding environments [3]. The most common features of the phytoplankton population is rapid
increase of biomass due to rapid cell proliferation and almost equally rapid decrease in population, separated by some
fixed time period. This type of rapid change in phytoplankton population density is known as ‘bloom’. The excitable
nature of blooms is the main characteristics in the plankton ecosystem. Due to the accumulation of high biomass
or to the presence of toxicity, some of these blooms, more adequately called “harmful algal blooms” (HABs, [4]),
are noxious to marine ecosystems or to human health and can produce great socioeconomic damage. Reduction of
∗ Corresponding author. Tel.: +91 33 2350 3862; fax: +91 33 2350 5207.
E-mail addresses: tsmath@rediffmail.com (T. Saha), mb math scc@yahoo.com (M. Bandyopadhyay).
1468-1218/$ - see front matter c 2007 Elsevier Ltd. All rights reserved.
doi:10.1016/j.nonrwa.2007.09.001
T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332 315
grazing pressure of zooplankton due to release of toxic substances by phytoplankton is one of the most interesting
topics of research in the last few decades [1,5–10]. Chattopadhyay et al. [8] proposed a mathematical model for
toxic phytoplankton (Noctilucca Scintillans belonging to the group Dinoflagellates of the division Dinophyta) and
zooplankton (Paracalanus belonging to the group Copepoda) interaction and explores the role of toxin producing
phytoplankton (TPP) behind harmful algal blooms. The general form of mathematical model they have considered is
the following two nonlinear coupled ordinary differential equations
dP
dt
= r P 1 −
P
k
− β f (P)Z (1a)
dZ
dt
= β1 f (P)Z − δZ − ρg(P)Z (1b)
where P ≡ P(t) is the density of toxin producing phytoplankton (TPP) population and Z ≡ Z(t) is the density of
zooplankton population at any instant of time ‘t’ subject to the non-negative initial condition P(0) = P0 ≥ 0 and
Z(0) = Z0 ≥ 0. In above model system ‘r’ is the intrinsic growth rate and ‘k’ is the environmental carrying capacity
of TPP population. ‘ f (P)’ represents the functional response for the grazing of phytoplankton by zooplankton
and ‘g(P)’ describes the distribution of toxic substance which ultimately contributes to the death of zooplankton
populations. β(>0) is the maximum uptake rate for zooplankton species, β1(>0) denotes the ratio of biomass
conversion (satisfying the obvious restriction 0 < β1 < β) and δ(>0) is the natural death rate of zooplankton. Finally
ρ(>0) denotes the rate of toxic substances produced by per unit biomass of phytoplankton. They have analyzed the
above model system by taking various combinations of functional response terms. In [10], the authors have analyzed
the following type of model system
dP
dt
= r P 1 −
P
k
− β f (P)Z (2a)
dZ
dt
= β1 f (P)Z − δZ − ρg(P(t − τ))Z (2b)
with f (P) = P and g(P) = P/(γ + P) based upon the assumption that the process of toxin liberation follows a
discrete variation and ‘τ’ is the discrete time period required for the maturation of phytoplankton cells to liberate toxic
substances. Based upon the model system (2), the authors studied the oscillation of phytoplankton and zooplankton
populations along with the stability condition for oscillatory behaviour.
Based upon their ideas, we intend to study a model system similar to (2) with the assumption that f (P) and g(P)
are described by same type of function, namely, Holling type-II functions. With this assumption the model system (2)
takes a new form as follows
dP
dt
= r P 1 −
P
k
− β f (P)Z (3a)
dZ
dt
= β1 f (P)Z − δZ − ρ f (P(t − τ))Z (3b)
subject to the initial condition P(θ) ≥ 0, Z(0) = Z0 ≥ 0 where P(θ) is a continuous function in −τ ≤ θ ≤ 0. As
phytoplankton is the most favourable food source for zooplankton within aquatic environments and the Holling type-II
functional form is a reasonable assumption to describe the law of predation we will assume Holling type-II functional
form for ‘ f (P)’ [10,11]. It is quite reasonable to assume that the law of grazing must be same whether it contributes
toward the growth of zooplankton species or it suppresses the rate of grazing due to presence of toxic substances. In
this paper we are mainly interested in the dynamics of the toxic-phytoplankton–zooplankton model system with the
same functional form for f (P) and g(P) in presence of time delay. So far as our knowledge goes, in most of the works
authors have assumed f (P) and g(P) are of different type. Based upon the above model system we intend to study
the periodic oscillatory nature of phytoplankton and zooplankton populations by considering a discrete time delay as
bifurcation parameter. After finding the conditions for stable periodic oscillation we also study the global existence of
periodic solutions for the model system (3).
The paper is organized as follows. In the next section we briefly mention the basic local dynamical behaviour of the
model system near steady states, where we make the assumption that the process of toxin liberation follows no time
316 T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332
lag, that is its effect is instantaneous to the grazers. In Section 3 we consider the delayed phytoplankton–zooplankton
model whose construction is based upon the assumption that process of toxin liberation follows a discrete time variation,
then we deduce the condition for exchange of stability behaviour around the coexisting equilibrium point through
a Hopf-bifurcation by considering the discrete time delay as bifurcation parameter. Section 4 consists of a purely algebraic
technique based on normal form theory and center manifold arguments for retarded functional differential
equations [12] to investigate the stability of small amplitude periodic solutions arising from Hopf-bifurcation. Using
the technique for global Hopf-bifurcation results due to Wu [13], we prove the global existence of periodic solutions
in Section 5. In Section 6, we give some numerical simulation results in support of the analytical findings. Finally, we
mention basic outcomes of our mathematical findings and their ecological significance in the concluding section.
2. Basic mathematical model: Stability analysis
The interaction of toxic-phytoplankton–zooplankton systems and their dynamical behaviour will be considered in
this paper based upon the following nonlinear ordinary differential equation model system (with τ = 0)
dP
dt
= r P 1 −
P
k
−
β P Z
γ + P
(4a)
dZ
dt
=
β1 P Z
γ + P
− δZ −
ρP Z
γ + P
(4b)
subject to the initial condition P(0) = P0 ≥ 0, Z(0) = Z0 ≥ 0. Here ‘r’ is the intrinsic growth rate of phytoplankton
species in absence of zooplankton and ‘γ ’ is the half saturation constant for a Holling type II functional response.
The biological significance of all other parameters have already been mentioned in the introduction. Now we make a
basic assumption that β1 > ρ, that is, the ratio of biomass consumed by zooplankton is greater than the rate of toxic
substance liberation by phytoplankton species. The vector fields associated with the model system (4) are smooth
within R2
+ and thus the existence and uniqueness criteria hold good and it is easy to verify that P(t), Z(t) > 0
whenever P(0), Z(0) > 0 [14]. Regarding the boundedness of the solution for the model system (4) we state the
following lemma:
Lemma 1. All the solutions of the model system (4) with the positive initial conditions P0, Z0 are uniformly bounded
within a region B ⊆ R2
+ where B = {(P, Z) ∈ R2
+ : β1 P + βZ = L
δ + ; L = (r+δ)β1k
2 }.
This lemma is obvious and we omit its proof (interested readers may consult the work [15]).
The equilibrium states for the model system (4) are given by E0(0, 0) (trivial equilibrium), E1(k, 0) (axial
equilibrium) and E∗(P∗, Z∗) (positive interior equilibrium), where P∗ = δγ
β1−ρ−δ and Z∗ = r
β (1 − P∗
k )(γ + P∗).
Thus for existence of E∗ we must have β1 > ρ + δ and P∗ < k. For local stability of these equilibrium states we
have to find the Jacobian matrix at these equilibrium states corresponding to the model system (4). We then have to
calculate eigenvalues of the corresponding Jacobian matrices to characterize the local stability of these equilibrium
states. The Jacobian matrix for the model system (4) at any point (P, Z) within positive quadrant is given by
J(P, Z) =



r −
2r P
k
−
βγ Z
(γ + P)2
−
β P
γ + P
(β1 − ρ)γ Z
(γ + P)2
(β1 − ρ)P
γ + P
− δ


 . (5)
It is easy to check that E0 is always a saddle point, E1 is locally asymptotically stable whenever E∗ does not exist and
a saddle point whenever E∗ exists. At E∗ the Jacobian matrix takes the form
J∗ =




r −
2r P∗
k
−
βγ Z∗
(γ + P∗)2
−
β P∗
γ + P∗
(β1 − ρ)γ Z∗
(γ + P∗)2
0



 .
Thus by Routh–Hurwitz criteria E∗ is locally asymptotically stable whenever k < 2P∗ +γ . In this case E∗ is a global
attractor also. The stability behaviour changes as k passes through the critical value k = k∗ = 2P∗ +γ and the system
exhibits Hopf-bifurcation. Thus the existence criteria for Hopf-bifurcation can be stated as follows:
T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332 317
Lemma 2. If k = k∗ = 2P∗ + γ with β1 > ρ + δ, then the model system (4) exhibits a Hopf-bifurcation around E∗.
This implies that a small amplitude periodic orbit bifurcates near E∗ whenever k passes through its critical value k∗.
A detailed analysis in this direction can be found in [8].
3. Mathematical model with delay: Local stability analysis
The delayed model system is developed based upon the fact that the liberation of toxic substances by phytoplankton
species is not an instantaneous phenomena rather it must be mediated by some time lag which is required for the
maturity of toxic-phytoplankton. If we assume ‘τ’ is the discrete time period required for the maturity of toxicphytoplankton
then model system (4) can be extended as a delay differential equation model system as follows
dP
dt
= r P 1 −
P
k
−
β P Z
γ + P
(6a)
dZ
dt
=
β1 P Z
γ + P
− δZ −
ρP(t − τ)Z
γ + P(t − τ)
(6b)
subject to the initial condition P(θ) ≥ 0, Z(0) = Z0 ≥ 0 where P(θ) for −τ ≤ θ ≤ 0. Under these assumptions the
model system (6) must posses an unique solution [16,17]. We will consider here the time delay as a control parameter
starting with the assumption that E∗ is locally asymptotically stable whenever τ = 0. For local asymptotic stability
analysis of E∗(P∗, Z∗) for the model system (6), we linearize the model system (6) with help of the transformations
P(t) = P∗ + x1(t), Z(t) = Z∗ + x2(t) (x1 and x2 are small perturbation terms whose second and higher order terms
can be neglected) and which results in the following system of linear delay differential equations :
dx1
dt
= ax1 + bx2 (7a)
dx2
dt
= cx1 + dx1(t − τ) (7b)
where the coefficients are given by
a =
r P∗
γ + P∗
1 −
2P∗
k
−
γ
k
, b = −
β P∗
γ + P∗
, c =
β1γ Z∗
(γ + P∗)2
, d = −
ργ Z∗
(γ + P∗)2
.
The characteristic equation corresponding to the linearized system (7) is given by
G(λ, τ) ≡ λ2
− aλ − b(c + de−λτ
) = 0. (8)
At τ = 0, E∗ is locally asymptotically stable whenever k < 2P∗ +γ . Now E∗ will be locally asymptotically stable
for τ ≥ 0 if the real part of G(λ, 0) is negative and G(iω, τ) = 0 for every real ω and τ ≥ 0. We assume that for
some values of τ(= 0), there exists a real number ω > 0 such that λ = iω is a root of the characteristic Eq. (8). Then
separating real and imaginary parts of G(iω, τ) = 0 we get
ω2
+ bc + bd cos ωτ = 0 (9)
bd sin ωτ − aω = 0. (10)
The above two equations can be combined as
ω4
+ ω2
(a2
+ 2bc) + b2
(c2
− d2
) = 0. (11)
Let us denote
∆ = (a2
+ 2bc)2
− 4b2
(c2
− d2
)
Then the roots of biquadratic equation (11) are given by
ω2
± =
1
2
−(a2
+ 2bc) ±
√
∆ .
318 T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332
Now as we need ω(>0) as a real quantity, we have to consider the following cases
(i) ∆ < 0 imply no purely imaginary roots of the form iω,
(ii) ∆ > 0, c2 > d2, a2 + 2bc > 0 imply no purely imaginary roots of the form iω,
(iii) ∆ > 0, c2 < d2, a2 + 2bc > 0 imply one purely imaginary root iω+,
(iv) ∆ > 0, c2 < d2, a2 + 2bc < 0 implies one purely imaginary root iω+,
(v) ∆ > 0, c2 > d2, a2 + 2bc < 0 imply two purely imaginary roots iω±.
Since the existence of the positive interior equilibrium E∗ demands c2 > d2, the cases (iii) and (iv) can be excluded.
For the cases (i) and (ii) we don’t have any purely imaginary roots of the characteristic Eq. (8). This shows that
the positive interior equilibrium point E∗ is absolutely stable (locally asymptotically stable for all τ ≥ 0) under the
parametric restrictions P∗ < k < 2P∗ + γ , β1 > ρ + δ with the conditions (i) and (ii).
We now consider the case (v), which implies c2 > d2, a2 + 2bc < 0; in this case there are two purely imaginary
roots given by ω± = 1
2 [−(a2 + 2bc) ±
√
∆] where 0 < ω− < ω+. From the Eqs. (9) and (10) we have
sin(ω±τ) =
aω±
bd
< 0 (12)
cos(ω±τ) =
(−a2 ±
√
∆)
2bd
. (13)
Since in this case (−a2 ±
√
∆) < 0, cos(ω±τ) < 0. Thus τ±
k = 1
ω±
[arcsin(aω±
bd ) + 2kπ], k = 0, 1, 2, . . .. Let us
define θ± = arcsin(aω±
bd ). Then π < θ− < θ+ < 3π
2 , 0 < ω− < ω+ and arcsin is decreasing. We have thus two
sequences of delays τ+
k = θ++2kπ
ω+
and τ−
k = θ−+2kπ
ω−
for which there are two purely imaginary roots of the Eq. (8).
We will now study how the real parts of the roots of (8) vary as τ varies in a small neighborhood of τ+
k as well as τ−
k .
Let λ = u + iω be a root of (8), then substituting λ = u + iω in the characteristic equation and separating real and
imaginary parts we have
H1(u, ω, τ) ≡ u2
− au − bc − ω2
− bde−uτ
cos ωτ = 0
H2(u, ω, τ) ≡ 2uω − aω + bde−uτ
sin ωτ = 0.
Now we have H1(0, ω±, τ±
k ) = H2(0, ω±, τk)± = 0. Also we have | J |(0,ω±,τ±
k ) > 0 where
J =



∂ H1
∂u
∂ H1
∂ω
∂ H2
∂u
∂ H2
∂ω


 .
Hence by the implicit function theorem H1(u, ω, τ) = 0 = H2(u, ω, τ) defines u, ω as a function of τ in a
neighborhood of (0, ω±, τ±
k ) such that u+(τ+
k ) = 0, ω+(τ+
k ) = ω+, du+
dτ (τ+
k ) > 0 and u−(τ−
k ) = 0ω−(τ−
k ) = ω−,
du−
dτ (τ−
k ) < 0. Now we are in a position to state the following theorem regarding the Hopf-bifurcation:
Theorem 1. For the delayed model system (6),
(i) if P∗ < k < 2P∗ + δ with β1 > ρ + δ, either (a2 + 2bc)2 < 4b2(c2 − d2) or a2 + 2bc > 0, and
(a2 + 2bc)2 − 4b2(c2 − d2) > 0 then E∗ is locally asymptotically stable for all τ ≥ 0
(ii) if P∗ < k < 2P∗ + γ with β1 > ρ + δ, and a2 + 2bc < 0, (a2 + 2bc)2 − 4b2(c2 − d2) > 0 then there exists a
positive integer n, such that the equilibrium E∗ switches n times from stability to instability to stability and so on
such that E∗ is locally asymptotically stable whenever τ ∈ [0, τ+
0 ) ∪ (τ−
0 , τ+
1 ) ∪ − − − − − − ∪(τ−
n−1, τ+
n ) and
is unstable whenever τ ∈ (τ+
0 , τ−
0 ) ∪ (τ+
1 , τ−
1 ) ∪ − − − − − − ∪(τ+
n−1, τ−
n−1) and τ > τ+
n . The model system (6)
undergoes a Hopf-bifurcation around E∗ for every τ = τ±
k .
4. Stability of periodic solutions
In the previous section we have obtained the conditions for which the model system (6) undergoes a Hopf
bifurcation by considering ‘τ’ as bifurcation parameter. In this section we determine the stability of periodic solutions
T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332 319
bifurcating from E∗ at τ = τ±
k . For this purpose we deduce the normal form of the model system (6) by following
the way as introduced by Faria and Magalhaes [12], choosing ‘τ’ as bifurcation parameter in a small neighborhood of
τ = τ±
k . Introducing the time scaling t → t
τ in Eq. (6) we get the following transformed system of equations
dP
dt
= τ r P 1 −
P
k
−
β P Z
γ + P
(14a)
dZ
dt
= τ
β1 P Z
γ + P
− δZ −
ρP(t − 1)Z(t)
γ + P(t − 1)
. (14b)
We now apply the transformation x1(t) = P(t) − P∗, x2(t) = Z(t) − Z∗, thus (14) reduces to
dx1(t)
dt
= τ ax1(t) + bx2(t) +
m≥2
1
m!
h(11)
m xm
1 (t) +
i+ j≥2
1
i! j!
h
(11)
i j xi
1(t)x
j
2 (t) (15a)
dx2(t)
dt
= τ cx1(t) + dx1(t − 1) +
m≥2
1
m!
h(2)
m xm
1 (t) +
i+ j≥2
1
i! j!
h
(21)
i j xi
1(t)x
j
2 (t)
+
i+ j≥2
h
(22)
i j xi
1(t − 1)x
j
2 (t) (15b)
where a, b, c, d have the same expression as given in Section 3 and the h’s are given as follows
h(11)
= r P 1 −
P
k
, h(2)
=
β1 P Z
γ + P
, h(21)
= Z
β1 P
γ + P
− δ , h(22)
= −Z δ +
ρP
γ + P
h(11)
m =
∂mh(11)
∂ Pm
(P∗,Z∗)
, h
(11)
i j =
∂i+ j h(11)
∂ Pi ∂ Z j
(P∗,Z∗)
, h(2)
m =
∂mh(2)
∂ Pm
(P∗,Z∗)
h
(21)
i j =
∂i+ j h(21)
∂ Pi ∂ Z j
(P∗,Z∗)
, h
(22)
i j =
∂i+ j h(22)
∂ Pi ∂ Z j
(P∗,Z∗)
.
The system of Eq. (15) in a phase space C := C([−1, 0], R2
+) can be written in matrix form as follows
dX
dt
= τ L(Xt ) + τ F(Xt ) (16)
where X ≡ X(t) is a column vector (x1(t), x2(t))T
, Xt ∈ C satisfying the condition Xt (θ) = X(t + θ), −1 ≤ θ ≤ 0.
L is a continuous linear function, mapping C into R2
+ and is defined by
L(φ) =
aφ1(0) + bφ2(0)
cφ1(0) + dφ1(−1)
, φ ∈ C.
This hypothesis implies that there exist an 2×2 matrix η(θ) for −1 ≤ θ ≤ 0 whose elements are of bounded variation
such that
L(φ) =
0
−1
[dη(θ)]φ(θ) (17)
and
dη(θ) =
a b
c 0
δ(θ) +
0 0
0 d
δ(θ + 1) dθ (18)
320 T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332
where δ denotes the Dirac-delta function. Also F : C → R2
+ is given by
F(φ) =





m≥2
1
m!
h(11)
m φm
1 (0) +
i+ j≥2
1
i! j!
h
(11)
i j φi
1(0)φ
j
2 (0)
m≥2
1
m!
h(2)
m φm
1 (0) +
i+ j≥2
1
i! j!
h
(21)
i j φi
1(0)φ
j
2 (0) +
i+ j≥2
h
(22)
i j φi
1(−1)φ
j
2 (0)





.
We now expand F(φ) in a Taylor series about φ
F(φ) = F0(φ) + F1(φ) +
1
2!
F2(φ) + · · · .
Comparing we get F0(φ) = F1(φ) = (0, 0)T and
1
n!
Fn(φ) =





1
n!
h(11)
n φn
1 (0) +
i+ j=n
1
i! j!
h
(11)
i j φi
1(0)φ
j
2 (0)
1
n!
h(2)
n φn
1 (0) +
i+ j=n
1
i! j!
h
(21)
i j φi
1(0)φ
j
2 (0) +
i+ j=n
h
(22)
i j φi
1(−1)φ
j
2 (0)





. (19)
Introducing a new parameter κ = τ − τk we rewrite (16) as
dX
dt
= τk L(Xt ) + F(Xt , κ) (20)
where F(Xt , κ) = κL(Xt )+(κ+τk)F(Xt ) and τk denotes any one of the critical values τ±
k (k = 0, 1, 2, . . .) at which
Eq. (8) has a pair of purely imaginary roots ±iω. The linear operator defined by (17) generates a strong continuous
semigroup of bounded linear operators with the infinitesimal generator
Aφ =



dφ
dθ
, −1 ≤ θ < 0
0
−1
[dη(θ)]φ(θ), θ = 0.
(21)
At τ = τk, the infinitesimal generator A has a pair of imaginary eigenvalues ±iσk where σk = ωτk. Let us consider
∧ = {−iσk, iσk}. Let P be the eigen-space associated with ∧ and we decompose C by ∧ as C = P ⊕Q where Q is the
complementary space of P in C and dimP = 2. Let us define C∗ = C([0, 1], R2
+
∗
), where R2
+
∗
is the two dimensional
vector space of row vectors and for any ψ in C∗ consider the adjoint bilinear form (., .) on C∗ × C associated with
the linear equation
(ψ, φ) = ψ(0)φ(0) −
0
−1
θ
0
ψ(ξ − θ)dη(θ)φ(ξ)dξ (22)
and the adjoint operator A∗ of A is defined as
A∗
ψ(s) =



−
dψ(s)
ds
0 ≤ s ≤ 1
0
−1
dη(s)ψ(−s) s = 0.
(23)
Let Φ = (Φ1, Φ2) be a basis for P and Ψ = Ψ1
Ψ2
be a basis for the dual space P∗ in C∗ associated with the
eigenvalues ±iσk of the adjoint equations. Then it can be normalized so that (Ψ, Φ) = I where I is a 2 × 2
identity matrix. We consider (20) in C([−1, 0] : C) and still denote it by C. Thus we write the 2 × 2 matrices Φ
and Ψ of the form Φ(θ) = Φ1(θ) Φ2(θ) , Φ1(θ) = eiσkθ v, Φ2(θ) = Φ1(θ), −1 ≤ θ ≤ 0 and Ψ(s) = Ψ1(s)
Ψ2(s) ,
Ψ1(s) = e−iσksuT, Ψ2(s) = Ψ1(s), 0 ≤ s ≤ 1 where v = v1
v2
and u = u1
u2
are vectors in C2 and also ˙Φ = ΦB,
T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332 321
where B is a 2 × 2 diagonal matrix B = iσk 0
0 −iσk
. With the help of (22) after some computation, we have
v =


1
iσk − aτk
bτk

 , u =


u1
bτku1
iσk

 (24a)
where
u1 =
iσk
(2 + dτke−iσk )iσk − aτk(1 + dτke−iσk )
. (24b)
Enlarge the phase space C to the space BC of the functions from [−1, 0] to R2
+ uniformly continuous on [−1, 0) and
with a jump discontinuity at 0. The elements of BC are of the form φ + X0α where φ ∈ C, α ∈ R2
+ that is we can
identify BC with C × R2
+. The projection of C upon P, φ → Φ(Ψ, φ) associated with decomposition C = P ⊕ Q
is now replaced by π : BC → P such that π(φ + X0α) = Φ[(Ψ, φ) + Ψ(0)α] leading to the decomposition
BC = P ⊕ ker π with the property Q ⊂ ker π. We decompose Xt = Φn(t) + st , n(t) ∈ C2, st ∈ Q = ker π ∩ C1
and Eq. (20) is equivalent to
˙n = Bn(t) + ψ(0)F(Φn + s, κ) (25a)
ds
dt
= AQ s + (I − π)X0 F(Φn + s, κ) (25b)
where X0 = X0(θ) is given by
X0(θ) =
I θ = 0
0 −1 ≤ θ < 0
(26)
and AQ : Q → ker π. We refer to [12] for detailed explanation of the notations described above. We consider the
Taylor expansion of the second terms of the right hand sides of Eq. (14)
ψ(0)F(Φn + s, κ) =
1
2!
f 1
2 (n, s, κ) +
1
3!
f 1
3 (n, s, κ) + h.o.t (27a)
(I − π)X0 F(Φn + s, κ) =
1
2!
f 2
2 (n, s, κ) +
1
3!
f 2
3 (n, s, κ) + h.o.t (27b)
where f 1
j (n, s, κ) and f 2
j (n, s, κ) are homogeneous polynomials in (n, s, κ) of degree j, j = 2, 3 with coefficients
in C2 and ker π, h.o.t. stands for higher order terms. The flow of (25) on the center manifold tangent to the invariant
subspace P at n = 0, κ = 0 is given by
˙n = Bn +
1
2
h1
2(n, 0, κ) +
1
3!
h1
3(n, 0, κ) + h.o.t. (28)
where h1
j , j = 2, 3 are homogeneous polynomials of degree j in (n, κ), which is in normal form. We consider
V
(m+p)
j (X) as the linear space of homogeneous polynomials of degree j in m + p variables, n = (n1, n2, . . . ., nm),
κ = (κ1, κ2, . . . , κp) with coefficients in X. For j ≥ 2, let Mj denote operator defined in V
m+p
j (Rm × ker π) with
values in the same space by
Mj (p, h) = (M1
j p, M2
j h)
(M1
j p)(n, κ) = [B, p(., κ)](n)
(M2
j h)(n, κ) = Dx h(n, κ)Bn − AQ (h(n, κ))
where [B, p(., κ)] denotes the Lie bracket defined by
[B, p(., κ)](n) = Dn p(n, κ)Bn − Bp(n, κ)
322 T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332
Here the operator M1
j acts on V 3
j (C2) and we consider the decomposition V 3
j (C2) = (Im M1
j ) ⊕ (ker M1
j ). We have
ker(M1
j ) = span{nqκlek : (q, λ) = λk, k = 1, 2, q ∈ N,l ∈ N, |(q,l)| = j} where {e1, e2} is the canonical basis for
C2. We have now
(Mj p)(n, κ) = iσk



n1
∂p1
∂n1
− n2
∂p1
∂n2
− p1
n1
∂p2
∂n1
− n2
∂p2
∂n2
+ p2


 .
Thus
ker(M1
2 ) = span
n1κ
0
,
0
n2κ
and
ker(M1
3 ) = span
n2
1n2
0
,
n1κ2
0
,
0
n1n2
2
,
0
n2κ2 .
From (27a)
f 1
2 (n, s, κ) = ψ(0)[2L(κ)(Φn + s) + F2(Φn + s, τk)].
We have now
Φn =
eiσkθ
v1n1 + e−iσkθ
v1n2
eiσkθ
v2n1 + e−iσkθ
v2n2
with the help of (19) we get
f
(1)
2 (n, 0, κ) =
2A1n1κ + 2A2n2κ + a20n2
1 + 2a11n1n2 + a02n2
2
2A1n2κ + 2A2n1κ + a02n2
1 + 2a11n1n2 + a20n2
2
where A1 = iσk
τk
uTv, A2 = −iσk
τk
uTv and bar denotes the complex conjugate.
a20 = τk
u1h
(11)
2 v2
1 + 2u1h
(11)
11 v1v2 + u2h
(2)
2 v2
1 + 2u2h
(21)
11 v1v2
+2u2h
(22)
11 e−iσk v1v2 + u2h
(22)
20 e−2iσk v2
1
a11 = τk
u1h
(11)
2 v1v1 + u1h
(11)
11 (v1v2 + v2v1) + u2h
(2)
2 v1v1 + u2h
(21)
11 (v1v2 + v2v1)
+u2h
(22)
11 (e−iσk v1v2 + eiσk v2v1) + u2h
(22)
20 v1v1
a02 = τk
u1h
(11)
2 v1
2
+ 2u1h
(11)
11 v1v2 + u2h
(2)
2 v1
2
+ 2u2h
(21)
11 v1v2
+2u2h
(22)
11 eiσk v1v2 + u2h
(22)
20 e2iσk v1
2 .
Following [12], the second order terms in (n, κ) of the normal form on the center manifold are given by
h1
2(n, 0, κ) = Projker(M1
2 ) f 1
2 (n, o, κ).
Thus we get
h1
2(n, 0, κ) =
2A1n1κ
2A1n2κ
.
For the cubic terms in the normal form, we have
h1
3(n, 0, κ) ∈ ker(M1
3 ).
In order to study the local stability of the periodic orbit arising from Hopf-bifurcation the terms o(|n|κ2) are irrelevant.
Hence we only need to compute the coefficients of n2
1n2
0
and
0
n1n2
2
. The canonical basis for V 3
2 (C2) consists of
T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332 323
following twelve elements
n2
1
0
,
n1n2
0
,
n1κ
0
,
n2
2
0
,
n2κ
0
,
κ2
0
0
n2
1
,
0
n1n2
,
0
n1κ
,
0
n2
2
,
0
n2κ
,
0
κ2 .
The image of each one of the elements of this basis under M1
2 /iσk are respectively
n2
1
0
, −
n1n2
0
,
0
0
, −3
n2
2
0
, −2
n2κ
0
, −
κ2
0
3
0
n2
1
,
0
n1n2
, 2
0
n1κ
, −
0
n2
2
,
0
0
,
0
κ2 .
We have
1
3!
h1
3(n, 0, κ) =
1
3!
Projker(M1
3 ) f 1
3 (n, 0, κ)
=
1
3!
Projker (M1
3 ) f 1
3 (n, 0, 0) + o(| n | κ2
)
where ker (M1
3 ) = span n2
1n2
0
,
0
n1n2
2
and f 1
3 (n, 0, 0) is the third order terms of the equation after computing the
normal form of the equation up to second order. Following the computation scheme given in [12] we have
f 1
3 (n, 0, 0) = f 1
3 (n, 0, 0) +
3
2
(Dn f 1
2 )U2 − (DnU2)h1
2
(n,0,0)
+
3
2
(Ds f 1
2 )g
(n,0,0)
We now compute 1
3! h1
3(n, 0, κ) step by step
Step 1: Computation for Projker (M1
3 ) (Dn f 1
2 )U2 (n,0,0)
Following the procedure in [12] we get
U2(n, 0, 0) = (M1
2 )−1
P1
I,2 f 1
2 (n, 0, 0)
where (M1
2 )−1 is the right inverse of (M1
2 ) and P1
I,2 f 1
2 (n, 0, 0) is the projection of f 1
2 (n, 0, 0) on the image space of
(M1
2 ). Now
f 1
2 (n, 0, 0) =
a20n2
1 + 2a11n1n2 + a02n2
2
a02n2
1 + 2a11n1n2 + a20n2
2
.
Hence
U2(n, 0, 0) =



1
iσk
(a20n2
1 − 2a11n1n2 −
1
3
a02n2
2)
1
iσk
1
3
a02n2
1 + 2a11n1n2 − a20n2
2


 .
Again
Dn f 1
2 =
2a20n1 + 2a11n2 2a11n1 + 2a02n2
2a02n1 + 2a11n2 2a11n1 + 2a20n2
and after some algebraic calculations we get
Projker (M1
3 ) (Dn f 1
2 )U2
(n,0,0)
=




2i
σk
(a20a11 − 2|a11|2
−
1
3
|a02|2
)n2
1n2
−
2i
σk
(a20a11 − 2|a11|2
−
1
3
|a02|2
)n1n2
2



 .
324 T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332
Step 2: Computation for Projker (M1
3 ) (DnU2)h1
2 (n,0,0)
We have h1
2(n, 0, 0) = 0 and so
Projker (M1
3 ) (DnU2)h1
2
(n,0,0)
= 0.
Step 3: Computation for Projker (M1
3 ) f 1
3 (n, 0, 0)
We have f 1
3 (n, 0, 0) = Ψ(0)F3(Φn, τk) and from (19)
1
3!
F3(Φn, τk) = τk



1
3!
h
(11)
3 φ3
1(0) +
1
2!
h
(11)
21 φ2
1(0)φ2(0)
1
3!
h
(2)
3 φ3
1(0) +
1
2!
h
(21)
21 φ2
1(0)φ2(0) +
1
2!
h
(22)
21 φ2
1(−1)φ2(0) +
1
3!
h
(22)
30 φ3
1(−1)


 .
After computing we have
Projker (M1
3 ) f 1
3 (n, 0, 0) =
3a21n2
1n2
3a21n1n2
2
where
a21 = τk



u1h
(11)
3 v2
1v1 + u1h
(11)
21 (v2
1v2 + 2v1v1v2) + u2h
(2)
3 v2
1v1
+ u2h
(21)
21 (v2
1v2 + 2v1v1v2) + u2h
(22)
21 (e−2iσk v2
1v2 + 2v1v1v2)
+ u2h
(22)
30 e−2iσk v2
1v1


 .
Step 4: Computation for Projker (M1
3 ) (Ds f 1
2 )g (n,0,0)
Here g is a second order homogeneous polynomial in (n1, n2, κ) with coefficients in Q . Let
g(n1, n2, κ) = g110n1n2 + g101n1κ + g011n2κ + g200n2
1 + g020n2
2 + g002κ2
where g = g1
g2
is the unique solution of
(M2
2 g)(n, κ) = (I − π) [2L(κ)(Φn) + F2(Φn, τk)]
which acts in V 3
2 (Q ). Following [12] we have
(M2
2 g)(n, κ) = Dng(n, κ)Bn − AQ (g(n, κ))
= Dng(n, κ)Bn − ˙g(n, κ) − X0 [L(τk)(g(n, κ)) − ˙g(n, κ)(0)]
= (I − π)X0 [2l(κ)(Φn) + F2(Φn, τk)] .
Thus from the above it follows that h = h(n, 0)(θ) can be evaluated by the system
˙g(n) − Dng(n)Bn = ΦΨ(0) [2L(0)(Φn) + F2(Φn, τk)] (29)
˙g(n)(0) − L(τk)(g(n)) = 2L(0)(Φn) + F2(Φn, τk) (30)
where we have used π(φ + X0κ) = Φ[(Ψ, φ) + Ψ(0)κ] and ˙g denotes the derivative of g(n)(θ) with respect to θ.
We have now
f 1
2 (n, s, 0) = Ψ(0) [2L(0)(Φn + s) + F2(Φn + s, τk)]
= τk
u1 u2
u1 u2
h
(11)
2 c2
1 + 2h
(11)
11 c1d1
h
(2)
2 c2
1 + 2h
(21)
11 c1d1 + 2h
(22)
11 c2d1 + h
(22)
20 c2
2
where
c1 = v1n1 + v1n2 + s1(0)
c2 = e−iσk v1n1 + eiσk v1n2 + s1(−1)
d1 = v2n1 + v2n2 + s2(0).
T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332 325
After some computation, we have the second order column matrix
Ds f 1
2 g
(n,0,0)
= 2τk
u1 u2
u1 u2
w1
w2
where
w1 = b0(h
(11)
2 g1
(0) + h
(11)
11 g2
(0) + h
(2)
2 g1
(0) + h
(11)
11 g2
(0)) + c0(h
(22)
11 g2
(0) + h
(22)
20 g1
(−1))
+ d0(h
(21)
11 g1
(0) + h
(11)
11 g1
(0) + h
(22)
11 g1
(−1))
w2 = [b0(h
(11)
2 g1
(0) + h
(11)
11 g2
(0) + h
(2)
2 g1
(0) + h
(11)
11 g2
(0)) + c0(h
(22)
11 g2
(0) + h
(22)
20 g1
(−1))
+ d0(h
(21)
11 g1
(0) + h
(22)
11 g1
(−1) + h
(11)
11 g1
(0))]
and
b0 = v1n1 + v1n2, c0 = e−iσk v1n1 + eiσk v1n2, d0 = v2n1 + v2n2
with
g1
(0) = g1 ∂s1(0)
∂s1
, g2
(0) = g2 ∂s2(0)
∂s2
, g1
(−1) = g1 ∂s1(−1)
∂s1
.
Also we may write
g1
(0) = g1
110(0)n1n2 + g1
200(0) + g1
020(0)n2
2
g1
(−1) = g1
110(−1)n1n2 + g1
200(−1) + g1
020(−1)n2
2
g2
(0) = g2
110(0)n1n2 + g2
200(0) + g2
020(0)n2
2.
Thus
Projker (M1
3 ) (Ds f 1
2 )g
(n,0,0)
=
2c3n2
1n2
2c3n1n2
2
where c3 is given by
c3 = uT
τk






[h
(11)
2 g1
110(0)v1 + h
(11)
2 g1
200(0)v1 + h
(11)
11 g1
110(0)v2 + h
(11)
11 g1
200(0)v2 + h
(11)
11 g2
110(0)v1]
[h
(11)
11 g2
200(0)v1 + h
(2)
2 g1
110(0)v1 + h
(2)
2 g1
200(0)v1 + h
(21)
11 g1
110(0)v2 + h
(21)
11 g1
200(0)v2
+h
(11)
11 g2
110(0)v1 + h
(11)
11 g2
200(0)v1 + h
(22)
11 g1
110(−1)v2 + h
(22)
11 g1
200(−1)v2
+h
(22)
11 e−iσk g2
110(0)v1 + h
(22)
11 eiσk g2
200(0)v1 + h
(22)
20 e−iσk g1
110(−1)v1 + h
(22)
20 eiσk g1
200(−1)v1]






.
Again from (29) and (30) we obtain the differential equation for g110(θ) and g200(θ) as
˙g110(θ) = φ1 φ2
2a11
2a11
(31)
˙g200(θ) − 2iσk g200(θ) = φ1 φ2
a20
a02
(32)
with the conditions
˙g110(0) − L(τk)(g110) = τk
a1
b1
(33)
˙g200(0) − L(τk)(g200) = τk
a2
b2
(34)
where
g110(θ) =
g1
110(θ)
g2
110(θ)
, g200(θ) =
g1
200(θ)
g2
100(θ)
326 T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332
and
a1 = 2h
(11)
2 v1v1 + 2h11
11(v1v2 + v2v1)
b1 = 2h
(2)
2 v1v1 + 2h
(21)
11 (v1v2 + v2v1) + 2h
(22)
11 (e−iσk v1v2 + eiσk v1v2) + 2h
(22)
20 v1v1
a2 = h
(11)
2 v2
1 + 2h
(11)
11 v1v2
b2 = h
(2)
2 v2
1 + 2h
(21)
11 v1v2 + 2h
(22)
11 e−iσk v1v2 + h
(22)
20 e−2iσk v2
1.
Solving (31) and (32) subject to the conditions (33) and (34) we get
g110(θ) =
2
iσk
a11eiσkθ
v − a11e−iσkθ
v + C1 (35)
g200(θ) = −
1
iσk
a20eiσkθ
v +
1
3
a02e−iσkθ
v + C2e2iσkθ
(36)
where C1 = c
(1)
1 c
(2)
1
and C2 = c
(1)
2 c
(2)
2
are as follows
c
(1)
1 = −
b1
c + d
c
(2)
1 =
1
b
ab1
c + d
− a1
c
(1)
2 =
bb2τ2
k + 2iσka2τk
2iσk(2iσk − aτk) + bτ2
k (c + de−2iσk )
c
(2)
2 =
a2τ2
k (c + de−2iσk ) + b2τk(2iσk − aτk)
2iσk(2iσk − aτk) + bτ2
k (c + de−2iσk )
.
Thus the above four steps together give the result
1
3!
h1
3(n, 0, κ) =
A3n2
1n2
A3n1n2
2
where A3 = i
2σk
(a20a11 − 2|a11|2 − 1
3 |a02|2) + 1
2 (a21 + c3). Hence the normal form (28) becomes
˙n = Bn +
A1n1κ
A1n2κ
+
A3n2
1n2
A3n1n2
2
+ o(|n|κ2
+ |n|4
). (37)
This normal form relative to P can be written in terms of the co-ordinates (x, y) by using the change of variables
n1 = x − iy, n2 = x + iy, where i =
√
−1. Finally using the polar co-ordinates x = r cos θ, y = r sin θ, the normal
form (37) becomes
˙r = k1rκ + k2r3
+ o(κ2
ρ + |(r, κ)|4
) (38a)
˙θ = −σk + o(|(r, k)|) (38b)
where k1 = Re(A1), k2 = Re(A3). Summarizing the results obtained above leads us to the following theorem,
Theorem 2. The flow of the Eq. (20) on the center manifold of the origin at κ = 0 is given by (38). Moreover if
k2 = 0 (generic Hopf bifurcation) the periodic orbits of Eq. (38) bifurcating from the origin r = 0, κ = 0 are given by
r(t, κ) = −
k1κ
k2
1
2
+ o(κ)
θ(t, κ) = −σk + o |κ|
1
2
and the stability of periodic orbits depends upon the following conditions
T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332 327
(i) if k1k2 < 0 (respectively k1k2 > 0) there exists a unique non trivial periodic orbit in the neighborhood of
(r, κ) = (0, 0) for κ > 0 (respectively κ < 0) and no periodic orbit for κ < 0 (respectively κ > 0),
(ii) the non trivial periodic orbit is stable if k2 < 0 and unstable if k2 > 0.
5. Global existence of periodic solutions
In the earlier sections we have established that the model system (6) undergoes a Hopf bifurcation at E∗ for τ = τ±
k
and also investigated the stability of bifurcating periodic solutions. In this section by using a global Hopf bifurcation
theorem due to Wu [13] we study the global continuation of periodic solutions bifurcating from the point (E∗, τ+
k ),
k = 0, 1, 2, . . . for the model system (6). For convenience we write the model system (6) as
˙Γ = F(Γt , τ, T ) (39)
where Γt (θ) = Γ(t + θ) ∈ C([−τ, 0], R2
+) and Γt = (Pt , Zt ). We use the following notations
= C([−τ, 0], R2
+)
= cl{(Γ, τ, T ) ∈ × R+ × R+ : Γ is a T -periodic solution of (39)}
N = {(Γ, τ, T ) : F(Γ, τ, T ) = 0}.
Let µ(E∗,τ+
k , 2π
ω+
) denote the connected component of (E∗, τ+
k , 2π
ω+
) in where the expression for τ+
k and ω+ are given
in Section 3.
Lemma 3. Whenever E∗ exists, all the nontrivial periodic solution of (39) are uniformly bounded within R2
+.
Proof. Let ξ1, ξ2, η1 and η2 be the real quantities defined by
P(ξ1) = min{P(t)} P(η1) = max{P(t)}
Z(ξ2) = min{Z(t)} Z(η2) = max{Z(t)}
where P(t) and Z(t) are nonconstant periodic solutions to (39). As our state space is R2
+, we will look only for
positive periodic solutions. From (6a) we get
0 = r −
r
k
P(η1) −
βZ(η1)
γ + P(η1)
< r 1 −
P(η1)
k
.
This implies
P(η1) < k.
From (6b)
0 =
β1 P(η2)
γ + P(η2)
− δ −
ρP(η2 − τ)
γ + P(η2 − τ)
<
β1 P(η2)
γ + P(η2)
− δ.
This shows
P(η2) >
δγ
β1 − δ
.
Hence
δγ
β1 − δ
< P(η2) < P(η1) < k. (40)
Again from (6a)
0 = r 1 −
Pη2
k
−
βZ(η2)
γ + P(η2)
< r −
βZ(η2)
γ + P(η2)
.
This implies
βZ(η2) < r (γ + P(η2)) < r(γ + k).
328 T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332
That is
Z(η2) <
r(γ + k)
β
. (41)
Hence by the relation (40) and (41) the lemma is confirmed.
Lemma 4. Suppose E∗ exists and k < 2P∗ +γ . Then there is no any nontrivial τ-periodic solution for the system (6).
Proof. Let us assume that the model system (6) has a nontrivial τ-periodic solution. Then the model system (3) has
also a periodic solution. For the model system (3) we note that the P-axis and Z-axis are invariant and the orbits of
the model system (3) do not intersect each other. There is no solution which will cross the co-ordinate axes. If there
be any periodic solution within the first quadrant, then there must be equilibrium in its interior and obviously E∗
is there. Now if k < 2P∗ + γ holds with the existence of E∗ then E∗ is globally asymptotically stable. Hence by
the Benedixon-Dulac criterion the system (3) does not have any periodic orbit under the restriction mentioned in the
lemma and the lemma is confirmed.
Theorem 3. Suppose that E∗ exists with k < 2P∗ + γ and a2 + 2bc < 0. Then for each τ > τ+
k (k = 0, 1, . . .),
system (6) has at least k + 1 periodic solutions.
Proof. The characteristic matrix of Eq. (6) at an equilibrium Γ = (P, Z) is given by
(Γ, τ, T )(λ) =
λ − a −b
−c − de−λτ
λ − e
(42)
where
a = r 1 −
2P
k
−
βγ Z
(γ + P)2
, b = −
β P
γ + P
, c =
β1γ Z
(γ + P)2
,
d = −
ργ Z
(γ + P)2
and e =
(β1 − ρ)P
γ + P
− δ
(Γ, τ, T ) is called a center if F(Γ, τ, T ) = 0 and det( (Γ, τ, T )(2π
T i)) = 0. A center is said to be isolated if it is
the only center in some neighbourhood of (Γ, τ, T ). The model system (6) has three equilibria E0(0, 0), E1(k, 0) and
E∗(P∗, Z∗). It is clear from Section 2 that the equilibria E0(0, 0) and E1(k, 0) are not centers. From the discusssion
about the local Hopf-bifurcation theorem in Section 3, it follows that (E∗, τ+
k , 2π
ω+
) is an isolated center and there exist
> 0, δ > 0 and a smooth curve λ : (τ+
k − δ, τ+
k + δ) → C such that det( (λ(τ))) = 0, | λ(τ) − iω+ |< for all
τ ∈ [τ+
k − δ, τ+
k + δ] and
λ(τ+
k ) = iω+,
dReλ(τ)
dτ τ=τ+
k
> 0.
Let
Ω , 2π
ω+
= (u, T ) : 0 < u < , T −
2π
ω+
< .
Now if | τ − τ+
k |≤ δ and (u, T ) ∈ ∂Ω , 2π
ω+
, then det( (E∗, τ, T )(u + 2π
T i)) = 0 iff u = 0, τ = τ+
k and T = 2π
ω+
.
Hence all the assumptions (A1) ∼ (A4) of Theorem 8 in [13] are satisfied for m = 1. Moreover if we define
H±(E∗, τ+
k , 2π
ω+
)(u, T ) = det( (E∗, τ+
k ± δ, T )(u + 2π
T i)) then we have the crossing number of isolated centers
(E∗, τ+
k , 2π
ω+
) as follows
γ1 E∗, τ+
k ,
2π
ω+
= degB H−
(E∗, τ+
k ,
2π
ω+
), Ω , 2π
ω+
− degB H+
E∗, τ+
k ,
2π
ω+
, Ω , 2π
ω+
= −1 = 0.
T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332 329
Thus by Theorem 8 of [13], we conclude that the connected component µ(E∗,τ+
k , 2π
ω+
) through (E∗, τ+
k , 2π
ω+
) in Σ is
unbounded. We will now show that the projection of µ(E∗,τ+
k , 2π
ω+
) onto τ-space is [τ+
k , ∞]. With the help of Lemma 4
the projection of µ(E∗,τ+
k , 2π
ω+
) onto τ-space is away from zero. Let us suppose that the projection of µ(E∗,τ+
k , 2π
ω+
) onto
τ-space is bounded. This implies that the projection of µ(E∗,τ+
k , 2π
ω+
) onto τ space must include an interval of the form
(0, τ). With the help of Lemma 4 this shows that the projection of µ(E∗,τ+
k , 2π
ω+
) onto T -space is bounded. Hence by
Lemma 3 the connected component µ(E∗,τ+
k , 2π
ω+
) is bounded and we have a contradiction and consequently the theorem
is confirmed.
6. Numerical results
In this section, we present some numerical results for some particular values of the parameters associated with
the model system (6). For ecological justification behind the choice of numerical values and related information,
interested readers may consult the literature by Chattopadhyay et al. [8]. We consider the following model system
dP
dt
= 2P 1 −
P
108
−
0.7P Z
5 + P
(43a)
dZ
dt
=
0.6P Z
5 + P
− .3654Z −
0.2P(t − τ)Z
5 + P(t − τ)
(43b)
which has an interior equilibrium point E∗(52.80346820809254, 84.40611033414001) and this set of parameter
values satisfies the local asymptotic stability conditions of E∗ in absence of time delay. Substituting P(t) =
52.80346820809254 + x1(t) and Z(t) = 84.40611033414001 + x2(t) in (43) and then linearizing we get
dx1
dt
= −0.0441006743738x1 − 0.63945x2,
dx2
dt
= 0.07578571428571x1 − 0.0252619047619x1(t − τ).
(44)
As the conditions ∆ > 0, c2 > d2 and a2 + 2bc < 0 are satisfied, we find two purely imaginary roots iω± with
ω+ = 0.24583783672861, ω− = 0.18585273350794
After some usual algebraic calculations one can find the minimum value of the delay parameter ‘τ’ for the model
system (43) for which the stability behaviour changes and the first critical values are given by
τ+
0 = 15.77200313837373, and τ−
0 = 19.7669696683186
such that, E∗ is locally stable for τ ∈ [0, 15.77200313837373) and is unstable for τ ∈ (15.77200313837373,
19.7669696683186). The eigenvectors as defined in (24) are given by
u =
0.42399156898581 + 0.09428853634629i
−0.2452543732444 + 1.1028465446809i
v =
1
−0.06896657185675 − 0.38445200833312i
.
Now one can calculate the relevant quantities to obtain the normal form by using any mathematical software
(e.g. MATLAB, MAPLE), which results in the desired normal form (see Eq. (38)) as follows
˙r = −0.02766106548628rκ − 0.00483827117482r3
+ o(κ2
ρ + |(r, κ)|4
) (45a)
˙θ = −3.87735513241468 + o(|(r, k)|). (45b)
Finally the stability determining quantities for Hopf-bifurcating periodic solutions are given by
k1 = −0.02766106548628,
k2 = Re(−0.00483827117482 − 0.00917806049283i) = −0.00483827117482.
330 T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332
Fig. 1. Solution of (43) for τ = 14.5 < τ+
0 showing E∗ = (52.80, 84.41) is stable.
Fig. 2. Hopf bifurcating periodic solution of (43) for τ = 16.2 > τ+
0 .
As k2 < 0, the periodic solution is stable. Interested readers may verify the results with other sets of parametric
values or with some realistic field data.
From our analytical findings we have seen that E∗ is locally asymptotically stable for τ < τ+
0 , Fig. 1 shows the
simulation result for the model system (43) with τ = 15.4 < τ+
0 . Interior equilibrium point looses its stability as τ
passes through its critical value τ = τ+
0 and a Hopf bifurcation occurs, a stable Hopf-bifurcating periodic solution is
depicted in Fig. 2. E∗ remains unstable whenever τ+
0 < τ < τ−
0 and regains its stability for τ > τ−
0 . The equilibrium
point E∗ remains locally asymptotically stable whenever the delay parameter lies in the range τ−
0 < τ < τ+
1 (see
Fig. 3). E∗ again switches from stability to instability as τ passes through τ = τ+
1 and an unstable solution for the
model system (43) is shown in Fig. 4. The numerical simulations we have done here illustrate the stable periodic
solution arising from Hopf bifurcation at τ = τ+
0 and the switching of stability that occurs as the magnitude of the
delay parameter increases gradually.
7. Discussion
Aquatic environments are not only a common habitat for phytoplankton as well as zooplankton but they constitute
the aquatic ecosystem as it is. Global increase of harmful phytoplankton blooms in the last two decades have received a
great deal of attention from several researchers, but the pattern and mechanism behind the planktonic blooms as well as
its possible control have not yet taken a proper shape and hence it emerges as an interesting subject area for theoretical
ecologists and their co-researchers. It is established by many of the present day researchers that one of the key factors
behind phytoplankton bloom and succession is the toxic substance released by some phytoplankton species. The
main objective of the present paper is to study the qualitative behaviour of an interacting phytoplankton–zooplankton
system in presence of toxic substances with the help of a delay differential equation model system with same type
T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332 331
Fig. 3. E∗ regains its stability for τ−
0 < τ = 20.5.
Fig. 4. E∗ switches from stability to instability at τ = 42.5 > τ+
1 .
of functional response associated with the growth of zooplankton due to grazing of phytoplankton and the term
which has a negative effect on the zooplankton growth due to presence of toxic substances. Analysis of the delayed
model system is based upon the assumption that the toxin liberation by the phytoplankton species follows a discrete
time variation and we have established that when the time delay ‘τ’ crosses a threshold value ‘τ±
k ’, the delayed
toxic-phytoplankton–zooplankton model system enters into a Hopf bifurcation and we have a periodic orbit around
a coexisting equilibrium point E∗ for k = 0, 1, 2, . . . .. We have shown that there are two possible cases for local
asymptotic stability of E∗ for the delayed model system, in the first possible case the positive interior equilibrium
E∗ is absolutely stable under certain parametric restrictions mentioned in the theorem of Section 3, (in this situation
bloom phenomena never arise for the model system under consideration) and the second possible case includes a
sequence of two time delays τ+
k (k = 0, 1, 2, . . . .) and τ−
k (k = 0, 1, 2, . . . .) under certain parametric restrictions
for which the interior equilibrium point E∗ is stable whenever τ ∈ [0, τ+
0 ) ∪ (τ−
0 , τ+
1 ) ∪ · · · and unstable for
τ ∈ (τ+
0 , τ−
0 )∪(τ+
1 , τ−
1 )∪· · ·. This phenomenon is known as switching of stability which arises for our model system.
The most interesting as well as mathematically important results we have presented in this paper is the stability criteria
for the Hopf-bifurcating periodic solution by considering the discrete time lag ‘τ’ as bifurcation parameter. In order
to understand the importance of the delay differential model system in studying the toxic-phytoplankton–zooplankton
interaction compared to the non-delayed version and the effect of discrete time delays on dynamical behaviour of the
model system we have considered ‘τ’ as bifurcation parameter. The local stability analysis of bifurcating periodic
solutions is based on normal form theory and we have obtained the criterion for stable and unstable periodic solutions.
According to the analytical results obtained in Section 4 our numerical example shows that the model system does
not have any periodic orbit whenever 0 ≤ τ < 15.77200313837373 and all trajectories eventually approach the
coexisting equilibrium point E∗ starting from any point within the positive quadrant of state space. A unique periodic
orbit emerges in the vicinity of E∗ when ‘τ’ crosses the critical value τ = 15.7720 and for τ < 19.7669696683186
332 T. Saha, M. Bandyopadhyay / Nonlinear Analysis: Real World Applications 10 (2009) 314–332
the periodic orbit is stable (see numerical results). This result indicates that there is a threshold limit of toxin liberation
by the phytoplankton species below which the system does not have any excitable nature and above which the system
shows excitability. These results are consistent with the dynamical behaviour of the toxic-phytoplankton–zooplankton
interaction model with discrete time lag ‘τ’, and agree well with some earlier works. We have established the global
existence of the bifurcating periodic solutions indicating the global oscillatory nature for the solution of the delayed
model system when it exhibits oscillatory behaviour. Numerical simulations which we have performed show the
switching of stability behaviour for the coexisting equilibrium point as the time delay crosses its different threshold
values; interested readers may check the stability of oscillatory solutions obtained for τ > τ+
1 .
There is still a tremendous amount of work to do with the phytoplankton–zooplankton model, especially in presence
of toxic substances. It would be interesting to study the emerging dynamical behaviour with various types of functional
response terms to describe the grazing of phytoplankton by zooplankton and the effect of toxic substances liberated
by the phytoplankton on their grazers in presence of discrete and continuous time delay in various terms associated
with the model system. On the other hand, the study of dynamical behaviour of planktonic interaction models with
varying magnitude of different parameters involved with the model system keeping the discrete time delay as a
constant and verification of these results with some practical data set may suggest some control mechanism for
harmful phytoplankton blooms. We hope, these issues will be well addressed in the near future and we leave them
for the subject matter of our future research. Finally we would like to remark that the subject area, namely, ‘harmful
phytoplankton blooms’ demands more in-depth and consistent research to understand the underlying mechanism
behind recurrent bloom formations of different phytoplankton species along with possible control mechanisms.
Acknowledgement
The authors are thankful to Prof. C.G. Chakrabarti, S.N. Bose Professor of Theoretical Physics, Department of
Applied Mathematics, University of Calcutta, for his continuous help and guidance throughout the preparation of this
paper.
References
[1] A.E. Abdllaoui, J. Chattopadhyay, O. Arino, Comparisons, by models, of some basic mechanisms acting on the dynamics of the zooplanktontoxic
phytoplankton systems, M3AS 12 (10) (2002) 1421–1451.
[2] E.P. Odum, Fundamentals of Ecology, Saunders, Philadelphia, 1971.
[3] J. Duinker, G. Wefer, Das CO2 Und Die Rolle Des Ozeans, Naturwissenschahten 81 (1994) 237–242.
[4] T. Smayda, What is a bloom? A commentary, Limnol. Oceonogr. 42 (1997) 1132–1136.
[5] Anderson, et al., Marine biotoxins and harmful algae: A national plan, Woods Hole Oceanographic Institution Technical Report WHOI
93-102, 1993.
[6] E.J. Buskey, C.J. Hyatt, Effects of the Texas (USA) brown tide alga on planktonic grazers, Marine Ecol. Prog. Ser. 126 (1995) 285–292.
[7] E.J. Buskey, D.A. Stockwell, Effects of a persistent brown tide on zooplankton population in the Laguno Madre of Southern Texas,
in: T.J. Smayda, Shimuzu (Eds.), Toxic Phytoplankton Blooms in the Sea, Elsevier, Amsterdam, 1993, pp. 659–666.
[8] J. Chattopadhyay, R.R. Sarkar, S. Mandal, Toxin producing plankton may act as a biological control for planktonic blooms: A field study and
mathematical modelling, J. Theoret. Biol. 215 (2002) 333–344.
[9] R.R. Sarkar, J. Chattopadhyay, The role of environmental stochasticity in a toxic phytoplankton-non-toxic phytoplankton zooplankton system,
Environmetrics 14 (2003) 775–792.
[10] J. Chattopadhyay, R.R. Sarkar, A. Abdllaoui, A delay differential equation model on harmful algal blooms in the presence of toxic substances,
IMA J. Math. Appl. Med. Biol. 19 (2002) 137–161.
[11] D. Ludwig, et al., Qualitative analysis of an insect outbreak system: The Spruce budworm and forest, J. Animal Ecol. 47 (1978) 315–328.
[12] T. Faria, L.T. Magalhaes, Normal forms for retarded functional differential equations with parameters and applications to Hopf bifurcations,
J. Differential Equations 122 (1995) 181–200.
[13] J. Wu, Symmetric functional differential equations and neural networks with memory, Trans. Amer. Math. Soc. 350 (1998) 4799–4838.
[14] G. Birkoff, G.C. Rota, Ordinary Differential Equations, Wiley, Boston, 1982.
[15] M. Bandyopadhyay, C.G. Chakrabarti, Deterministic and stochastic analysis of a nonlinear prey-predator system, J. Biol. Sys. 11 (2) (2003)
161–172.
[16] J.K. Hale, Theory of Functional Differential Equations, Springer-Verlag, New York, 1977.
[17] Y. Kuang, Delay Differential Equation With Applications in Population Dynamics, Academic, New York, 1993.