\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 $$\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}$$ 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: $$\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}$$ 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 $$\label{e2.1} \frac{dx(t)}{dt}=Ax(t)+H(x),$$ 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 $$\label{e2.2} \begin{gathered} \frac{dx}{dt}=Ax+Bu,\\ y=-Cx,\\ u=g[y(t-\tau); \tau], \end{gathered}$$ 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 $$\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}.$$ 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 \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} It follows from \eqref{e2.4} and \eqref{e2.5} that $$\label{e2.7} E^2+F^2=S^2+T^2.$$ Then $$\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,$$ where \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} 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 $$\label{e2.10} \tau_0=\frac{1}{\omega_0}[\arccos\frac{ES-FT} {E^2-F^2}+2k\pi](k=0, 1, 2,\dots).$$ 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, $$\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}$$ 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}