\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2011 (2011), No. 148, pp. 1--8.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2011 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2011/148\hfil Bifurcations for a phytoplankton model]
{Bifurcations for a phytoplankton model with time delay}

\author[C. Xu \hfil EJDE-2011/148\hfilneg]
{Changjin Xu} 

\address{Changjin Xu \newline
Guizhou Key Laboratory of Economics System Simulation,
School of Mathematics and Statistics \\
Guizhou College of Finance and Economics\\
Guiyang 550004, China}
\email{xcj403@126.com}

\thanks{Submitted April 13, 2011. Published November 3, 2011.}
\thanks{Supported by grant 10961008 from the NNSF and
the Doctoral Foundation of Guizhou \hfill\break\indent
College of Finance and Economics (2010), China}
\subjclass[2000]{34K20}
\keywords{Phytoplankton model; Hopf bifurcation; stability; frequency domain;
\hfill\break\indent Nyquist criterion}

\begin{abstract}
 Applying a frequency domain approach, we investigate a
 phytoplankton model with time delay. We use the delay as a
 bifurcation parameter; as it passes through a sequence of
 critical values, Hopf  bifurcation occurs.
 A family of periodic solutions bifurcate from
 the equilibrium when the bifurcation parameter exceeds
 a critical value. Some numerical simulations illustrate
 our theoretical results.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\allowdisplaybreaks

\section{Introduction}

 Phytoplankton and Zooplankton are a single celled organisms that
drift with the currents on the surface of open oceans. They play an
important role in stabilizing the environment; for example, they
are the staple item for the food web and they are recyclers of most of
energy that flows through the ocean ecosystem. Phytoplankton
systems have received much attention from biologists and
mathematicians \cite{c1,h1}. In 2010, Dhar and Sharma \cite{d1}
investigated the stability of the  phytoplankton system
\begin{equation} \label{e1.1}
\begin{gathered}
   \frac{dP_s(t)}{dt}=rP_s[1-\frac{P_s}{K}]-{\alpha}P_sP_i+{\gamma}P_i,\\
   \frac{dP_i(t)}{dt}={\alpha}P_sP_i-{\beta}P_i,
   \end{gathered}
   \end{equation}
where $P_s,P_i$ are the population densities of susceptible and
infected phytoplankton at any instant of time $t$. $r$ is the
intrinsic growth rate of the population of susceptible
phytoplankton, $K$ is the carrying capacity of the population of
susceptible phytoplankton, $\alpha$ is the disease contact rate of
the disease phytoplankton population, $\beta$ is the removal rate of
the disease phytoplankton population, out of which $\gamma$ fraction
of infected phytoplankton rejoin the susceptible phytoplankton
population. In details, one can see \cite{d1}. Taking into account that
there is a certain time delay during the process of reproduction(for
example, egg formation will take $\tau$ units of time before
hatching), the dynamic behavior of the system not only is affected
by the current state of the system, but also the past state of the
system, i.e., there exists inherent lag in the system. Based on this
view of point, in this paper, we will revise system \eqref{e1.1} as
follows:
\begin{equation}\label{e1.2}
\begin{gathered}
   \frac{dP_s(t)}{dt}=rP_s[1-\frac{P_s(t-\tau)}{K}]
-{\alpha}P_sP_i+{\gamma}P_i,\\
   \frac{dP_i(t)}{dt}={\alpha}P_sP_i-{\beta}P_i.
   \end{gathered}
\end{equation}

In this article, we investigate Hopf bifurcation for
 system \eqref{e1.2}.
It is worth pointing out that many early work on Hopf bifurcation of
the delayed differential equations is used the state-space
formulation for delayed differential equations, known as the ``time
domain'' approach. But there exists another approach that comes from
the theory of feedback systems known as frequency domain method
which was initiated and developed by Allwright \cite{a1},
Mees and Chua \cite{m1}
and Moiola and Chen \cite{m2,m3} and is familiar to control engineers.
This
alternative representation applies the engineering feedback systems
theory and methodology: an approach described in the ``frequency
domain''--the complex domain after the standard Laplace transforms
having been taken on the state-space system in the time domain. This
new methodology has some advantages over the classical time-domain
methods \cite{l1,l2}. A typical one is its pictorial characteristic that
utilizes advanced computer graphical capabilities thereby bypassing
quite a lot of profound and difficult mathematical analysis.

In this paper, we will devote our attention to finding the Hopf
bifurcation point for models \eqref{e1.2} by means of the frequency-domain
approach. We found that if the coefficient $\tau$ is used as a
bifurcation parameter, then Hopf bifurcation occurs for the model
\eqref{e1.2}. That is, a family of periodic solutions bifurcates from the
equilibrium when the bifurcation parameter exceeds a critical value.
Some numerical simulations are carried out to illustrate the
theoretical analysis. We believe that it is the first time to
investigate Hopf bifurcation of the model \eqref{e1.2} using the
frequency-domain approach.

  The remainder of the paper is organized as follows: in Section 2, applying the frequency-domain
approach formulated by Moiola and Chen \cite{m3}, the existence of Hopf
bifurcation parameter is determined and shown that Hopf bifurcation
occurs when the bifurcation parameter exceeds a critical value. In
Section 3, some numerical simulation are carried out to verify the
correctness of theoretical analysis result. Finally, some
conclusions and discussions are included in Section 4.


\section{Stability of the equilibrium and local Hopf
bifurcations}


In model \eqref{e1.2}, we assume that the  condition
\begin{itemize}
\item[(H1)] $(K\alpha-\beta)(\beta-\gamma)>0$.
\end{itemize}

It is obvious that system \eqref{e1.2} has a unique positive
equilibrium  $E_*(P_s^*, P_i^*)$, where,
$$
P_s^*=\frac{\beta}{\alpha}, \quad
P_i^*=\frac{r\beta(K\alpha-\beta)}{K\alpha^2(\beta-\gamma)}.
$$
We can rewrite the nonlinear system \eqref{e1.2} as a matrix
form
\begin{equation} \label{e2.1}
 \frac{dx(t)}{dt}=Ax(t)+H(x),
\end{equation}
where $x=( P_s(t), P_i(t))^T$,
 $$
A=\begin{pmatrix}
 r & 1  \\
 0 & -\beta  \end{pmatrix},
\quad
H(x)= \begin{pmatrix}
         -\frac{rP_sP_s(t-\tau)}{K}-{\alpha}P_sP_i \\
         {\alpha}P_sP_i
       \end{pmatrix}.
$$
Choosing the coefficient $\tau$ as a bifurcation and introducing a
``state-feedback control'' $u=g[y(t-\tau); \tau]$, where
$y(t)=(y_1(t), y_2(t))^T$, we obtain a linear system with a
non-linear feedback as follows
\begin{equation} \label{e2.2}
\begin{gathered}
   \frac{dx}{dt}=Ax+Bu,\\
   y=-Cx,\\
   u=g[y(t-\tau); \tau],
   \end{gathered}
\end{equation}
 where
$$
B=C= \begin{pmatrix}
               1 & 0 \\
               0 & 1
             \end{pmatrix}, \quad
u=g[y(t-\tau),\tau]
= \begin{pmatrix}
 -\frac{ry_1y_1(t-\tau)}{K}-{\alpha}y_1y_2 \\
                           {\alpha}y_1y_2
                         \end{pmatrix}.
$$
Next, taking Laplace transform on \eqref{e2.2}, we obtain
the standard transfer matrix of the linear part of the system:
$$
G(s;\tau)=C[sI-A]^{-1}B.
$$
Then
\begin{equation} \label{e2.3}
 G(s;\tau)= \begin{pmatrix}
                \frac{1}{s-r} & \frac{r}{(s-r)(s+\beta)} \\
                0 & \frac{1}{s+\beta}
              \end{pmatrix}.
\end{equation}
If this feedback system is linearized about the equilibrium
$y=-C(P_s^*, P_i^*)^T$, then the Jacobian of \eqref{e2.3} is
\[
J(\tau)=\frac{\partial{g}}{\partial{y}}\Big|_{y=\tilde{y}=-C(P_s^*,
P_i^*)^T}
=   \begin{pmatrix}
    \frac{r(P_s^*+P_s^*e^{-s\tau})}{K}+\alpha{x_2^*} & \alpha{x_1^*}\\
    -\alpha{x_2^*} & -\alpha{x_1^*}
  \end{pmatrix}.
\]
Let
\begin{align*}
h(\lambda,s;\tau)&= \det|\lambda{I}-G(s;\tau)J(\tau)| \\
&= \lambda^2+\big[\frac{\alpha{P_s^*}}{s+\beta}
 +\frac{r\alpha{P_i^*}}{(s-r)(s+\beta)}
 -\frac{1}{s-r}\big(\frac{rP_s^*(1+e^{-s\tau})}{K}
 +\alpha{P_i^*}\big)\big]\lambda \\
&\quad +\frac{\alpha{P_i^*}}{s+\beta}
 \big[\frac{\alpha{P_s^*}}{s-r}-\frac{r\alpha{P_s^*}}{(s-r)
 (s+\beta)}\big]=0.
\end{align*}
Applying the generalized   Nyquist stability criterion with
$s=i\omega$, we obtain the following  results.

\begin{lemma}[Moiola, Chen, 1996] \label{lem2.1}
If an eigenvalue of the corresponding Jacobian of the nonlinear
system, in the time domain, assumes a purely imaginary value
$i\omega_0$ at a particular $\tau=\tau_0$, then
the corresponding eigenvalue of the constant matrix
$G(i\omega_0;\tau_0)J(\tau_0)$ in the frequency domain must assume
the value $-1+i0$ at $\tau=\tau_0$.
\end{lemma}

To apply Lemma \ref{lem2.1}, let
$\hat{\lambda}=~\hat{\lambda}(i\omega; \tau)$ be the eigenvalue of
$G(i\omega;\tau)J(\tau)$ that satisfies $~\hat{\lambda}(i\omega_0;
\tau_0)=-1+0i$. Then
$$
h(-1,i\omega_0;\tau_0)=0.
$$
Separating the real and imaginary parts and rearranging, we obtain
\begin{gather} \label{e2.4}
E\cos\omega_0\tau_0+F\sin\omega_0\tau_0=S,\\
F\cos\omega_0\tau_0+E\sin\omega_0\tau_0=T,\label{e2.5}
\end{gather}
where
\begin{equation} \label{e2.6}
\begin{gathered}
E=rP_s^*(\beta^2-\omega_0^2), \quad
F=-2r\beta\omega_0{P_s^*},\\
\begin{aligned}
S&= Kr(\omega_0^2-\beta^2)-2K\beta\omega_0^2+K\alpha{P_s^*}
 (\beta{r}+\omega_0^2)\\
&\quad -rK\alpha\beta{P_i^*}
 -(\beta^2-\omega_0^2)(rP_s^*+K\alpha{P_i^*}),
\end{aligned}\\
\begin{aligned}
T&= Kr\alpha\omega_0{P_i^*}-K\omega_0(\beta^2-\omega_0^2)
 +2Kr\beta\omega_0\\
&\quad +K\alpha{P_s^*}\omega_0(\beta-r)
 +2\beta\omega_0(rP_s^*+\alpha{K}P_i^*).
\end{aligned}
\end{gathered}
\end{equation}
It follows from \eqref{e2.4} and \eqref{e2.5} that
\begin{equation} \label{e2.7}
E^2+F^2=S^2+T^2.
\end{equation}
Then
\begin{equation} \label{e2.8}
K^2\omega_0^6+\theta_1\omega_0^4+\theta_2\omega_0^3+\theta_3\omega_0^2
+\theta_4=0,
\end{equation}
where
\begin{equation}
\begin{gathered}
\theta_1= (Kr-2K\beta+K{\alpha}P_s^*+rP_s^*
 +K{\alpha}P_i^*)^2-r^2(P_s^*)^2, \\
\theta_2= 2K[Kr{\alpha}P_i^*-K\beta^2+2Kr\beta+K{\alpha}
 P_s^*(\beta-r)+2\beta(rP_s^*+K{\alpha}P_i^*)], \\
\begin{aligned}
\theta_3&= 2(Kr-2K\beta+K{\alpha}P_s^*+rP_s^*+K{\alpha}P_i^*)
 (Kr\alpha\beta{P_s^*}-Kr\beta^2-Kr\alpha\beta{P_i^*}), \\
      &\quad -2r^2\beta^2(P_s^*)^2,
\end{aligned}\\
\begin{aligned}
 \theta_4&= (Kr\alpha\beta{P_s^*}-Kr\beta^2-Kr\alpha\beta{P_i^*})^2
 +\theta_3^2-r^2\beta^4(P_s^*)^2
 +\big[Kr{\alpha}P_i^*-K\beta^2\\
&\quad +2Kr\beta+K{\alpha}P_s^*(\beta-r)+2\beta(rP_s^*
 +K{\alpha}P_i^*)]^2 -r^2\beta^*(P_s^*)^2.
\end{aligned}
\end{gathered}\label{e2.9}
\end{equation}
It is easy to that that if $\theta_4<0$, then \eqref{e2.8} has
at least one positive root (say $\omega_0$). By \eqref{e2.8},
we can compute the value of $\omega_0$ by means of Matlab software.
Then from \eqref{e2.4} and \eqref{e2.5}, we obtain
\begin{equation} \label{e2.10}
\tau_0=\frac{1}{\omega_0}[\arccos\frac{ES-FT}
{E^2-F^2}+2k\pi](k=0, 1, 2,\dots).
\end{equation}
In the sequel, we will consider the transversality condition for
Hopf bifurcation of system \eqref{e1.2}.

In view of the definition of $h(\lambda,s;\tau)$, we have
\begin{align*}
[\frac{d\lambda}{d\tau}]^{-1}
&= \frac{2K^3(s-r)\lambda+rP_s^*(1+e^{-s\tau})(s+\beta)+K\alpha{P_i^*}(s+\beta)}
{K^2rP_s^*(e^{-s\tau})\lambda(s+\beta)} \\
&\quad -\frac{K^3(s-r)\alpha{P_s^*}+Kr\alpha{P_i^*}}
{K^2rP_s^*(e^{-s\tau})\lambda(s+\beta)}.
\end{align*}
Thus we obtain
$$
[\frac{d\lambda}{d\tau}]^{-1}=\frac{C+iD}{A+iB},
$$
which leads to
$$
\operatorname{Re}[\frac{d\lambda}{d\tau}]_{\tau=\tau_0}^{-1}
=\frac{AC+BD}{A^2+B^2},
$$
where
\begin{gather}
A = -K^2rP_s^*(\beta\omega_0\sin\omega_0\tau_0
 -\omega_0^2\cos\omega_0\tau_0),\label{e2.11}\\
B = -K^2rP_s^*(\omega_0^2\sin\omega_0\tau_0
 +\beta\omega_0\cos\omega_0\tau_0),\label{e2.12}\\
\begin{aligned}
C &=  2K^3(r\beta-\omega_0^2)-K^3r\alpha{P_s^*}
 +Kr\alpha{P_i^*}-K\alpha\beta{P_i^*} \\
&\quad -rP_s^*[\beta(1+\cos\omega_0\tau_0)
 +\omega_0\sin\omega_0\tau_0],
\end{aligned} \label{e2.13} \\
D= 2K^3(\beta+r)+K^3\omega_0-rP_s^*[\omega_0(1+\cos\omega_0\tau_0)
 -\beta\sin\omega_0\tau_0]
  -K\omega_0\alpha {P_i^*}.\label{e2.14}
\end{gather}
To obtain our main result, we assume
\begin{itemize}
\item[(H2)] $AC+BD\neq 0$.
\end{itemize}


\begin{theorem}[Existence of Hopf bifurcation parameter] \label{thm2.1}
 Let $\theta_4, A, B, C, D$ be  defined by
\eqref{e2.9}, \eqref{e2.11}, \eqref{e2.12}, \eqref{e2.13},
 \eqref{e2.14}, respectively.
For \eqref{e1.2}, if $\theta_4<0$ and conditions {\rm (H1)--(H2)} hold,
then Hopf bifurcation point of system \eqref{e1.2} is
$$
\tau_0=\frac{1}{\omega_0}[\arccos\frac{ES-FT}
{E^2-F^2}+2k\pi](k=0, 1, 2,\dots),
$$
where $\omega_0$  is positive real roots of
\eqref{e2.8}, and  $E, F, S, T $ are defined by \eqref{e2.6}.
\end{theorem}


\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig1}
\includegraphics[width=0.6\textwidth]{fig2}
\includegraphics[width=0.6\textwidth]{fig3}
\end{center}
\caption{ 
Trajectories  and phase for \eqref{e3.1} with
$\tau=2.1<\tau_0\approx 2.17$. The positive equilibrium is
asymptotically stable. The initial value is (0.5, 0.5).
}\label{fig1-3}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig4}
\includegraphics[width=0.6\textwidth]{fig5}
\includegraphics[width=0.6\textwidth]{fig6}
\end{center}
\caption{Trajectories  and phase for \eqref{e3.1}
with $\tau=2.4>\tau_0\approx 2.17$. Hopf bifurcation occurs from the
positive equilibrium. The initial value is (0.5, 0.5)}\label{fig4-6}
\end{figure}


\section{Numerical Examples}

In this section, we shall carry out numerical simulations for
supporting our theoretical analysis. As an example, We consider
system \eqref{e1.2} with $r=0.5$, $K=2$, $\alpha=2$,
$\gamma=0.4 \beta=2$; that is,
\begin{equation} \label{e3.1}
\begin{gathered}
\frac{dP_s(t)}{dt}=0.5P_s[1-\frac{P_s(t-\tau)}{2}]-2P_sP_i+0.4P_i,\\
\frac{dP_i(t)}{dt}=2P_sP_i-2P_i.
\end{gathered}
\end{equation}
By Theorem \ref{thm2.1}, we obtain $\tau_0\approx 2.17$.
Numerical simulations for $\tau=2.1$ are shown in
Figure \ref{fig1-3}. Thus we conclude
that when $\tau<\tau_0\approx 2.17$, system \eqref{e3.1} is
asymptotically stable.
Numerical simulations for $\tau=2.4$ are shown in Figure
\ref{fig4-6}. Thus we conclude that when $\tau>\tau_0\approx 2.17$,
system \eqref{e3.1} undergoes a Hopf bifurcation that occurs near
the positive equilibrium. Therefore $\tau_0\approx 2.17$
is a supercritical Hopf bifurcation point.

\subsection*{Conclusions and discussions}
In this paper, we
investigated a class of phytoplankton model with time delay. By
choosing the coefficient $\tau$  as a bifurcating parameter and
analyzing the associating characteristic equation.
It is found that a Hopf bifurcation occurs when the bifurcating
parameter $\tau$ passes through a critical value.
Considering computational complexity, the direction and the
stability of the bifurcating periodic orbits for system
\eqref{e1.2} have not been investigated. It is beyond the scope
of the present paper and will be further investigated elsewhere
in the near future.


\begin{thebibliography}{10}

\bibitem{a1} D. J. Allwright;
\emph{Harmonic balance and the Hopf bifurcation
theorem}, Math. Proc. Combridge Philos. Soc. 82 (1977) 453-467.

\bibitem{c1} J. Chattopadhayay, R. R. Sarkar, S. Mandal;
\emph{Toxin producting Plankton may act
as a biological control for Planktonic bloom-Field study and
mathematical modelling}, J. Biol. Theor. 215(3) (2002) 333-344.

\bibitem{d1} J. Dhar, A. K. Sharma;
\emph{ The role of viral infection in phytoplankton dynamics
with the inclusion of incubation class}, Nonlinear Anal.: Hybird syst. 4 (2010) 9-15. %%%%

\bibitem{h1} M. Horst, F. M. Hilker, V. Sergei, S. V. Petrovskii,
B. Klaus;
\emph{Oscillations and waves
in avirally infected plankton system},
Ecol. Complex. 1 (3) (2004) 211-223.

\bibitem{l1} X. F. Liao, S. W. Li;
\emph{Hopf bifurcation on a two-neuron system with distributed delays:
a frequency domain approach}, Nonlinear Dynam. 31 (2003) 299-326.

\bibitem{l2} X. F. Liao, S. W. Li, G. R. Chen;
\emph{Bifurcation Analysis on  a Two-neuron system with distributed
delays in the frequency domain},  Neural Netw. 17 (2004) 545-561.

\bibitem{m1} A. I. Mees, L. O. Chua;
 \emph{The Hopf bifurcation theorem and its
applications to nonlinear oscillations in circuits and systems},
IEEE, Trans. Circuits. Syst., 26 (1979) 235-254

\bibitem{m2} J. L. Moiola, G. R. Chen;
\emph{Hopf bifurcation analysis: a
frequency domain approach}, Singapore: World Scientific, 1996.
.
\bibitem{m3} J. L. Moiola, G. R. Chen;
\emph{Frequency domain approach to
computational analysis of bifurcations and limit cycles: a
tutorial}, Int. J. Bifur. Chaos 3 (1993) 843-867.


\end{thebibliography}

\end{document}