\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2009(2009), No. 76, pp. 1--10.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2009 Texas State University - San Marcos.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2009/76\hfil Hopf bifurcation in simple food chain model] {Hopf bifurcation for simple food chain model with delay} \author[M. Cavani, T. Lara, S. Romero\hfil EJDE-2009/76\hfilneg] {Mario Cavani, Teodoro Lara, Sael Romero} % in alphabetical order \address{Mario Cavani \newline Departamento de Matem\'{a}ticas, N\'{u}cleo de Sucre, Universidad de Oriente, Cuman\'{a} 6101, Venezuela} \email{mcavani@sucre.udo.edu.ve} \address{Teodoro Lara \newline Departamento de F\'{\i}sica y Matem\'{a}ticas, Universidad de los Andes, N\'ucleo Univeristario Rafael Rangel, Trujillo, Venezuela} \email{tlara@ula.ve} \address{Sael Romero \newline Departamento de Matem\'{a}ticas, N\'{u}cleo de Sucre, Universidad de Oriente, Cuman\'{a} 6101, Venezuela} \email{sromero@sucre.udo.edu.ve} \thanks{Submitted May 11, 2009. Published June 16, 2009.} \subjclass[2000]{34D99} \keywords{Simple food chain model; Hopf bifurcation; Holling type I; attractor} \begin{abstract} In this article we consider a chemostat-like model for a simple food chain where there is a well stirred nutrient substance that serves as food for a prey population of microorganisms, which in turn, is the food for a predator population of microorganisms. The nutrient-uptake of each microorganism is of Holling type I (or Lotka-Volterra) form. We show the existence of a global attractor for solutions of this system. Also we show that the positive globally asymptotically stable equilibrium point of the system undergoes a Hopf bifurcation when the dynamics of the microorganisms at the bottom of the chain depends on the history of the prey population by means of a distributed delay that takes an average of the microorganism in the middle of the chain. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \section{Introduction} We consider the food chain model \begin{equation} \label{eq1} \begin{gathered} S'(t)=(S^{0}-S(t))D-\frac{b}{\gamma } S(t)X(t), \\ X'(t)=X(t)(bS(t)-D-\frac{d}{\eta } Y(t)), \\ Y'(t)=Y(t)(dX(t)-D), \end{gathered} \end{equation} where $S(0)=S_{0}\geq 0$, $X(0)=X_{0}\geq 0$, $Y(0)=Y_{0}\geq 0$. These equations are in the form of the chemostat model \cite{smithwaltman95}. The meaning of the variables and parameters is as follows: $S(t)$ denotes the substrate concentration, $X(t)$ is the concentration of a prey microorganism population that grows eating the substrate, and $Y(t)$ the concentration of a predator microorganism population that eats the prey microorganism. The functional responses of the species $X(t)$ and $Y(t)$ are of the so called Holling type I (or Lotka-Volterra) form (see \cite{holling59}), where the parameters $b$ and $d$ denote the per capita growth rate of the prey and the predator respectively; $\gamma $ and $\eta $ represent the growth yield constants of the microorganisms $X(t)$ and $Y(t)$ respectively. The parameter $S^{0}$ is the concentration in the feed bottle and $D$ denotes the input rate from the feed bottle and the washout rate from the growth chamber. This model may be considered as a specialization of some of the systems given in \cite{BWolkowicz86}, we shall perform a wide description of the global dynamic of this model in order to point out the differences between the dynamic of this system and the one corresponding to a system with a distributed delay that we propose in order to model the case when the existence of a significative time lag in the growth of the predator microorganism is considered. The recognition of time delays in the growth response of a population to changes in the environment has led to extensive theoretical and experimental studies, however there has been little emphasis in distributed delays in chemostat models, \cite{Wolkowicz1} is very important reference in this direction. Thus, we are assuming in a more realistic fashion that the growth of the predator is influenced by the amount of prey in the past. More precisely, we suppose as for example in \cite{Cushing} or \cite{Wolkowicz1}, that the predator grows up depending on the weight average over the past by mean of the the function $Z(t)$ given by the following integral \begin{equation} \label{chv1} Z(t):=\int_{-\infty }^{t}dX(\tau )Y(\tau ) e^{-D(t-\tau )}(\alpha e^{-\alpha (t-\tau ) })d\tau ,\quad \alpha >0, \end{equation} therefore, we have the integro-differential system \begin{equation} \label{eq2} \begin{gathered} S'(t)=(S^{0}-S(t))D-\frac{b}{\gamma } S(t)X(t), \\ X'(t)=X(t)(bS(t)-D-\frac{d}{\eta } Y(t)), \\ Y'(t)=\int\nolimits_{-\infty }^{t}dX(\tau)Y(\tau )e^{-D(t-\tau )}(\alpha e^{-\alpha (t-\tau )})d\tau -DY(t), \end{gathered} \end{equation} $S(0)=S_{0}\geq 0$, $X(0)=X_{0}\geq 0$, $Y(0)=Y_{0}(t)=\varphi (t)\geq 0$ $(t\leq 0)$. Clearly these assumptions imply that the influence of the past\ is fading away exponentially and the number $\frac{1}{\alpha }$ could be interpreted as the measure of the influence of the past. So, to smaller $\alpha >0$, the interval in the past in which the values of $X$ are taken, is bigger (\cite{Cushing,Farkasmotion,MacD1,Wolkowicz1}). Also we assume that initial function $\varphi$ is in $ BC_{+}$, the Banach space of the bounded and continuous functions from $(-\infty ,0]$ to $\mathbb{R}_{+}$. To make the model in more treatable way we perform the change of variables: \begin{gather*} \overline{t}= tD,\quad \overline{S} = \frac{S}{S^{0}},\quad \overline{X}= \frac{X}{\gamma S^{0}},\quad \overline{Y}= \frac{Y}{\eta \gamma S^{0}}, \\ \overline{b}=\frac{bS^{0}}{D},\quad \overline{d}= \frac{\gamma dS^{0}}{D}, \quad \overline{\alpha } = \frac{\alpha }{D}. \end{gather*} Omitting the bars, the nondimensional version of models (\ref{eq1}) and (\ref{eq2}) can be rewritten, respectively, as \begin{equation} \label{ndeq1} \begin{gathered} S'(t)=1-S(t)-bS(t)X(t), \\ X'(t)=X(t)(bS(t)-1-dY(t)), \\ Y'(t)=Y(t)(dX(t)-1), \end{gathered} \end{equation} and \begin{equation} \label{ndeq2} \begin{gathered} S'(t)=1-S(t)-bS(t)X(t), \\ X'(t)=X(t)(bS(t)-1-dY(t)), \\ Y'(t)=\int\nolimits_{-\infty }^{t}dX(\tau )Y(\tau )\alpha e^{-(\alpha +1)(t-\tau ) }d\tau -Y(t). \end{gathered} \end{equation} In this article we show the existence of the global attractor for the solutions of the foregoing systems. We also show that the positive globally asymptotically stable equilibrium point of \eqref{ndeq1} loses its stability when we model the food chain by \eqref{ndeq2}. In this case the equilibrium of positive coordinates undergoes a Hopf bifurcation and more realistic periodic solutions gain the stability. \section{A simple food chain without delay} Here we show some properties of system \eqref{ndeq1}. \begin{lemma} \label{lem1} The positive cone, $\mathbb{R}_{+}^{3}$, is positively invariant with respect to \eqref{ndeq1}. \end{lemma} \begin{proof} As we can see, if $S(t^{\ast })=0$ for some $t^{\ast }\geq 0$ then $S(t)\geq 0$ for all $t\geq t^{\ast }$. The positiveness of the functions $X(t)$ and $Y(t)$ are straightforward checked once the corresponding equations are considered. \end{proof} Note that by adding the three equations of \eqref{ndeq1} and defining $W(t)=1-S(t)-X(t)-Y(t)$, we obtain a single equation \begin{equation*} W'(t)=-W(t) \end{equation*} with $W(0)>0$. It is easy to verify that $\lim_{t\to \infty}W(t)=0$ and that the convergence is exponential. This implies that \eqref{ndeq1} has the property of pointwise dissipativity in the sense that there exists a bounded set $B$ to which the solutions eventually enter and remain. Thus we have shown the following \begin{lemma}[Dissipativity] \label{lem2} System \eqref{ndeq1} is pointwise dissipative. Moreover, the attractors of the solutions are located on the manifold \begin{equation} \Sigma =\{ (S,X,Y)\in \mathbb{R}_{+}^{3}:S+X+Y=1\} . \label{sigma} \end{equation} \end{lemma} The pointwise dissipative property implies the existence of a unique global attractor of \eqref{ndeq1} which must lie in the manifold $\Sigma $. \begin{lemma} \label{lem3} If $d\leq 1$, the predator population $Y(t)$ dies out. \end{lemma} \begin{proof} By taking into account the equation for $Y(t)$ in \eqref{ndeq1} and applying comparative arguments the result follows. \end{proof} By virtue of the previous lemma we shall assume for the rest of this article that \begin{equation} d>1. \label{h1} \end{equation} Another important conclusion of the Dissipativity lemma is that the system can be simplified by eliminating one variable. In fact by taking \begin{equation*} S(t)=1-X(t)-Y(t), \end{equation*} we obtain the following system of two ordinary differential equations \begin{equation} \label{twoeq} \begin{gathered} X'(t)=(b-1)X(t)-bX^{2}(t)-(b+d)X(t)Y(t), \\ Y'(t)=Y(t)(dX(t)-1). \end{gathered} \end{equation} Figure \ref{fig1} shows some numerical examples for the above equation with several values of $b$ and $ d$. \begin{figure}[ht] \label{fig1} \begin{center} \includegraphics[width=0.7\textwidth]{fig1} \end{center} \caption{$b< 1$ and both species extinguish} \end{figure} \begin{lemma} \label{lem4} If $b\leq 1$ in \eqref{twoeq}, the prey $X(t)$ and predator populations $ Y(t)$ die out. \end{lemma} The proof runs in the same fashion as in Lemma \ref{lem3}. As a consequence of the previous result we will assume in the sequel that $b>1$. System \eqref{twoeq} has three equilibrium points given by \begin{equation*} E_{0}=(0,0), \quad E_{1}=\big(\frac{b-1}{b},0\big), \quad E_{2}=\big(\frac{1}{d},\frac{d(b-1)-b}{d(b+d)}\big). \end{equation*} The stability properties of these points are summarized as follows. \begin{theorem} \begin{itemize} \item[(i)] If $b<1$, then $E_{0}$ is the unique equilibrium point of \eqref{twoeq} in the positive cone and is globally asymptotically stable. \item[(ii)] If $b=1$, the point $E_{0}$ undergoes a node-saddle bifurcation. And for $b>1$ the equilibrium point $E_{1}$ shows up, and it is globally asymptotically stable for $1\frac{d}{d-1}. \label{desig1} \end{equation} \end{itemize} \end{theorem} \begin{proof} Parts (i) and (ii) follow immediately from a linear analysis of the equilibria solutions $E_{0}$ and $E_{1}$. To show (iii) we use Dulac's Criterion \cite{Farkasmotion}. Let $f_{1}(X,Y)$ and $f_{2}(X,Y)$ be the corresponding functions in the right hand side of \eqref{twoeq} for $X'(t)$ and $Y'(t)$ respectively. In our case we look for a function of the form $h(X,Y)=X^{\alpha }Y^{\delta }$ such that the expression $\frac{\partial hf_{1}}{\partial X}+\frac{\partial hf_{2}}{\partial Y}$ is not zero and does not change its sign while $X>0$ and $Y>0$. In doing so, we see that \begin{equation} \label{dulac} \begin{aligned} &\frac{\partial (hf_{1})(X,Y)}{\partial X}+\frac{\partial (hf_{2})(X,Y)}{ \partial Y} \\ &=[ (\alpha +1)(b-1)-(\delta +1)] X^{\alpha}Y^{\delta } +[ (\delta +1)d-b(\alpha +2)] X^{\alpha +1}Y^{\delta }\\ &\quad -(b+d)(\alpha +1)X^{\alpha }Y^{\delta +1}. \end{aligned} \end{equation} Therefore, while $X>0$ and $Y>0$, \eqref{dulac} will be negative if we can find values of $\alpha $ and $\delta $ such that \begin{gather*} (\alpha +1)-\frac{\delta +1}{b-1} \leq 0, \\ (\alpha +2)-\frac{d}{b}(\delta +1) \geq 0. \end{gather*} But it is an easy task to guarantee the existence of such a values $\alpha ^{\ast }$ and $\delta ^{\ast }$ for which the previous\ inequalities hold, and therefore the function $h(X,Y)=X^{\alpha ^{\ast }}Y^{\delta ^{\ast }}$ satisfies the conditions we were looking for. Hence by applying the Dulac's Criterion to this function we conclude that \eqref{twoeq} has no periodic orbits and the Poincar\'{e}-Bendixson theory implies that equilibrium point $E_{2}$ is globally asymptotically stable. \end{proof} \section{A simple food chain with delay} Now we considered the delayed model given by \eqref{ndeq2}. If we take $Z(t)$ in (\ref{chv1}) as a change of variable, then following system shows up. \begin{equation} \label{eq3} \begin{gathered} S'(t)=1-S(t)-bS(t)X(t), \\ X'(t)=X(t)(bS(t)-1-dY(t)), \\ Y'(t)=Z(t)-Y(t), \\ Z'(t)=\alpha dX(t)Y(t)-(\alpha +1)Z(t), \end{gathered} \end{equation} with $S(0)=S_{0}\geq 0$, $X(0)=X_{0}\geq 0$, $Y(0)=Y_{0}\geq 0$, $Z(0)=\varphi (0)\geq 0$. The relations between the solutions of this system and those of \eqref{ndeq2} are as in the corresponding description given in \cite{CLS}. The properties of positiveness and pointwise dissipativeness hold as in the non-delay model and consequently similar results can be stated and proved. \begin{lemma} The positive cone, $\mathbb{R}_{+}^{4}$, is positively invariant with respect to \eqref{eq3}. \end{lemma} Now we set $U(t)=1-S(t)-X(t)-Y(t)-\frac{Z(t)}{\alpha }$, and obtain single equation, \begin{equation*} U'(t)=-U(t) \end{equation*} with $U(0)>0$. From here, $\lim_{t\to \infty }U(t)=0$ and the convergence is exponential, and again as before \eqref{eq3} has the property of pointwise dissipativity. \begin{lemma}[Dissipativity] \label{lem7} The system \eqref{eq3} is pointwise dissipative. Moreover the attractors of the system are located on the manifold \begin{equation} \Lambda =\{ (S,X,Y,Z)\in \mathbb{R}_{+}^{4}:S+X+Y+\frac{Z}{\alpha } =1\} . \label{lambda} \end{equation} \end{lemma} The pointwise dissipative property implies the existence of a unique global attractor of \eqref{eq3} which must lie in the manifold $\Lambda $. As before the system can be simplified by the elimination of one variable. In this case we take \begin{equation*} S(t)=1-X(t)-Y(t)-\frac{Z(t)}{\alpha } \end{equation*} and obtain a system of three ordinary differential equations, \begin{equation} \label{3deq} \begin{gathered} X'(t)=(b-1)X(t)-bX^{2}(t) -(b+d)X(t)Y(t)-\frac{b}{\alpha }X(t)Z(t), \\ Y'(t)=Z(t)-Y(t), \\ Z'(t)=\alpha dX(t)Y(t)-(\alpha +1)Z(t). \end{gathered} \end{equation} See illustration in Figure \ref{fig2}. \begin{figure}[ht] \label{fig2} \begin{center} \includegraphics[width=0.7\textwidth]{fig1} \end{center} \caption{$b< 1$ and all species extinguish} \end{figure} This system has three equilibrium points: \[ P_{0} =(0,0,0),\quad P_{1}=(\frac{b-1}{b},0,0), \quad P_{2}=(X^{\ast },Y^{\ast },Y^{\ast }), \] where \[ X^{\ast } =\frac{\alpha +1}{\alpha d},\quad Y^{\ast }=\frac{\alpha d(b-1)-b(\alpha +1)}{d(b(\alpha +1)+\alpha d)}. \] The expression for $P_{2}$ makes sense only when $\alpha d(b-1)-b(\alpha +1)>0$, inequality that is equivalent to \begin{equation}\label{Eq} X^{\ast }<\frac{b-1}{b}\,. \end{equation} For $b>1$, the right hand side of \ref{Eq} implies $X^{\ast }<1$. The stability properties of these equilibrium points are summarized as follows. \begin{theorem} \label{thm8} \begin{itemize} \item[(i)] If $b<1$, then $P_{0}$ is the unique equilibrium point of \eqref{3deq} in the positive cone and it is globally asymptotically stable. \item[(ii)] If $b=1$, $P_{0}$ undergoes a node-saddle bifurcation. For $b>1$, $P_{1}$ appears and it is globally asymptotically stable for $b(\alpha +1)-\alpha d(b-1)\geq 0$, or equivalently \begin{equation*} 1\frac{(\alpha +1)b}{\alpha (b-1)} \label{ineqE2} \end{equation} the point $P_{2}$ appears and has positive coordinates. \end{itemize} \end{theorem} The proof of the above theorem follows from a linear analysis of the equilibria solutions. Stability properties of the equilibrium point $E_{2}$ are given in the following result. \begin{theorem} \label{thm9} There exists a value $\ d^{\ast }$ satisfying \eqref{ineqE2} and if $d0, \quad \lim_{d\to +\infty }\Psi (d)=-\infty , \\ \lim_{d\to -\infty }\Psi '(d)=-\infty ,\quad \Psi '(0)>0,\quad \lim_{d\to +\infty }\Psi '(d)=-\infty \,. \end{gather*} Therefore, there exists a unique value $d_{1}$ such that $\Psi '(d)>0 $ for $0d_{1}$. This means that $\Psi (d)$ increases for $0\leq dd_{1}$. But the we can guarantee the existence of a unique value $d^{\ast }>d_{1}$, such that $\Psi (d)>0$ for $0d^{\ast }$. Moreover, since $c_{2}>2\alpha ^{2}(\alpha +1) $, \begin{equation*} d^{\ast }>d_{1}>\frac{(\alpha +1)b}{\alpha (b-1)}>1, \end{equation*} so Routh-Hurwitz implies that $E_{2}$ is locally asymptotically stable for $dd^{\ast }$. For $d=d^{\ast }$ the equilibrium undergoes a Hopf bifurcation. In fact, \begin{equation*} \lambda ^{3}+a_{2}(d^{\ast })\lambda ^{2}+a_{1}(d^{\ast })\lambda +a_{0}(d^{\ast })=(\lambda ^{2}+a_{1}(d^{\ast }))(\lambda +a_{2}(d^{\ast })). \end{equation*} To check the Hopf bifurcation, we set $\lambda _{1}(d^{\ast })$ as the root of (\ref{poly}) that assume the value $i\omega $, $\omega ^{2}=a_{1}(d^{\ast })$, at $d^{\ast }$ and by \begin{equation*} F(\lambda ,d)=\lambda ^{3}+a_{2}(d)\lambda ^{2}+a_{1}(d)\lambda +a_{0}(d) \end{equation*} the characteristic polynomial (\ref{poly}) as a function of $d$. Thus the derivative of the implicit function $\lambda _{1}$ at $d^{\ast } $ is \[ \lambda _{1}'(d^{\ast }) =-\frac{F_{d}'(i\omega ,d^{\ast})}{F_{\lambda }'(i\omega ,d^{\ast })} =-\frac{a_{2}'(d^{\ast })(i\omega )^{2}+a_{1}'(d^{\ast })(i\omega )+a_{0}'(d^{\ast })}{3(i\omega )^{2}+2a_{2}(d^{\ast })(i\omega )+a_{1}'(d^{\ast })}, \] and \begin{align*} ( \mathop{\rm Re}(\lambda _{1}(d^{\ast }))_{d}' &=\mathop{\rm Re}((\lambda _{1})_{d}'(d^{\ast })) \\ &=-\frac{(a_{1}(d^{\ast })a_{2}(d^{\ast })-a_{0}(d^{\ast }))_{d}'}{ a_{1}(d^{\ast })(1+a_{2}^{2}(d^{\ast }))} \\ &=-\frac{a_{2}(d^{\ast })b\left[ X_{d}^{\ast \prime }(\alpha +2+dY^{\ast })+X^{\ast }(Y^{\ast }+dY_{d}^{\ast \prime })\right] }{a_{1}(d^{\ast })(1+a_{2}^{2}(d^{\ast }))}, \end{align*} after some boring calculations we can see that the quantity between brackets is negative, therefore $(\mathop{\rm Re}(\lambda _{1}(d^{\ast }))_{d}'>0 $, and so the transversality condition required for the Hopf bifurcation holds. To verify the properties of stability of the periodic orbit we need to translate \eqref{3deq} and locate its origin at the equilibrium point $(X^{\ast },Y^{\ast },Y^{\ast })$. In this case the new system is \begin{equation} \label{eqtrans} \begin{gathered} x' =-bX^{\ast }x-(b+d)X^{\ast }y-\frac{b}{\alpha }X^{\ast }z-bx^{2}-(b+d)xy-\frac{b}{\alpha }xz, \\ y' =-y+z, \\ z' =\alpha dY^{\ast }x+\alpha dX^{\ast }y-(\alpha +1)z+\alpha dxy. \end{gathered} \end{equation} We perform the following change of variables to transform the above system in a normal form \begin{equation*} \begin{bmatrix} x \\ y \\ z \end{bmatrix} =T\begin{bmatrix} u \\ v \\ w \end{bmatrix} , \end{equation*} where $T$ is the $3\times 3$ matrix whose first and second columns are the real and imaginary part of the eigenvector associated with the eigenvalue $\lambda _{1}(d^{\ast })$ and the third column is the eigenvector associated to the eigenvalue $\lambda _{0}(d^{\ast })=-a_{2}(d^{\ast})$. Indeed \begin{equation*} T=\begin{bmatrix} A & B & C \\ \frac{1}{\omega } & 1 & 1 \\ \frac{1}{\omega } & 0 & D \end{bmatrix} \end{equation*} where \begin{gather*} A =-\frac{bX^{\ast }(\omega ^{2}+(\alpha d+b(\alpha +1))X^{\ast })}{ \omega ^{2}+bX^{\ast }},\quad B=\frac{1}{\alpha }+\frac{\omega (\omega ^{2}+(\alpha d+b(\alpha +1))X^{\ast })}{\omega ^{2}+bX^{\ast }}, \\ C =-\frac{b+d}{b}+(\frac{1}{\alpha }-\frac{\alpha +2+bX^{\ast }}{ bX^{\ast }})(\alpha +1+bX^{\ast }),\quad D=-(\alpha +1+bX^{\ast }). \end{gather*} Therefore, the system (\ref{eqtrans}) takes the form \begin{equation} \label{eq4} \begin{gathered} u' =\omega v+G_{1}(u,v,w), \\ v' =-\omega u+G_{2}(u,v,w), \\ w' =-a_{2}(d^{\ast })w+G_{3}(u,v,w), \end{gathered} \end{equation} with \begin{align*} \begin{bmatrix} G_{1}(u,v,w) \\ G_{2}(u,v,w) \\ G_{3}(u,v,w) \end{bmatrix} &=T^{-1}\begin{bmatrix} F_{1}(Au+Bv+Cw,\frac{u}{\omega }+v+w,\frac{u}{\omega }+Dw) \\ F_{2}(Au+Bv+Cw,\frac{u}{\omega }+v+w,\frac{u}{\omega }+Dw) \\ F_{3}(Au+Bv+Cw,\frac{u}{\omega }+v+w,\frac{u}{\omega }+Dw) \end{bmatrix} \\ &=\begin{bmatrix} H_{1}u^{2}+H_{2}v^{2}+H_{3}w^{2}+H_{4}uv+H_{5}vw+H_{6}uw \\ 0 \\ \alpha dCw^{2}+L_{1}vw+L_{2}uw+ \alpha dBv^{2}+L_{3}uv+\alpha dA\frac{u^{2}}{\omega } \end{bmatrix} \end{align*} and \begin{gather*} L_{1} =\alpha d(B+C), \quad L_{2}=\alpha d(A+\frac{C}{\omega }),\quad L_{3}=\alpha d(A+\frac{B}{\omega }), \\ H_{1} =-A(bA+\frac{1}{\omega }(\frac{b}{\alpha }+b+d)),\quad H_{2}=-B(b+d+bB), \\ H_{3} =-C(b+d+bC+\frac{b}{\alpha }D),\quad H_{4}=-A(b+d+2bB)-\frac{B}{\omega }(b+d+\frac{b}{\alpha }), \\ H_{5} =-(B(b+d+\frac{b}{\alpha }D)+C(b+d+2bB)), \\ H_{6} =-A(b+d+b(2C+\frac{D}{\alpha })) -\frac{C}{\omega }(b+d+\frac{b}{\alpha }). \end{gather*} Now we determine approximately the $d=d^{\ast }$-section of the center manifold $M$ which is tangent to the $(u,v)$- plane at the origin. This is $w=h(u,v),$ $h(0,0)=h_{u}'(0,0)=h_{v}'(0,0)=0$, and $h$ sufficiently smooth. Then \begin{equation} w=h(u,v)=\frac{1}{2}(h_{11}u^{2}+2h_{12}uv+h_{22}v^{2})+o(\mid (u,v)\mid ^{2}). \label{centr} \end{equation} Restricting the system to the center manifold, if $(u(t),v(t),w(t))$ is a solution of \eqref{eq4} near the origin with a value on $M$, then it stays locally in $M$, i.e., \begin{equation*} w(t)\equiv h(u(t),v(t)). \end{equation*} But then \begin{equation} w'(t)-h_{u}'(u(t),v(t))u'(t)-h_{v}'(u(t),v(t))v'(t)\equiv 0, \label{eq5} \end{equation} so by using \eqref{eq4}, omitting terms of order at least three and equating the coefficients to zero, \begin{gather*} h_{11} = 2\frac{-2\alpha dB\omega ^{3}-L_{3}a_{2}(d^{\ast })\omega ^{2}+2\alpha dA\omega ^{2}+\alpha dAa_{2}^{2}(d^{\ast })}{\omega a_{2}^{3}(d^{\ast })}, \\ h_{12} = \frac{1}{a_{2}^{2}(d^{\ast })}(-2\alpha dA+2\alpha dB\omega +a_{2}(d^{\ast })L_{3}), \\ h_{22} = 2\frac{2\omega \alpha dA-2\alpha dB\omega ^{2}-L_{3}a_{2}(d^{\ast })\omega +a_{2}^{2}(d^{\ast })\alpha dB}{ a_{2}^{3}(d^{\ast })}. \end{gather*} To restrict \eqref{eq4} to the $d=d^{\ast }$-section of the center manifold $M$, we introduce new coordinates, $y_{1}=u$, $y_{2}=v$, $y_{3}=w-h(u,v)$ and \ref{eq4}) becomes \begin{equation} \label{eqcentr} \begin{aligned} y_{1}' &= \omega y_{2}+\frac{1}{2}h_{22}H_{5}y_{2}^{3}+( h_{12}H_{5}+\frac{1}{2}h_{22}H_{6})y_{1}y_{2}^{2}+ ( \frac{1}{2}H_{5}h_{11}+H_{6}h_{12})y_{1}^{2}y_{2} \\ &\quad +\frac{1}{2}h_{11}H_{6}y_{1}^{3} +H_{1}y_{1}^{2}+H_{2}y_{2}^{2}+H_{4}y_{1}y_{2}+O(\mid y\mid ^{4}) \\ y_{2}' &=-\omega y_{1}. \end{aligned} \end{equation} Applying the Bautin's formula \cite[Lemma 7.2.7]{Farkasmotion}, we can see that \begin{equation} \frac{4\omega }{3\pi }V_{y_{1}y_{1}y_{1}}'''(0,0) =(3h_{11}h_{12}H_{5}+h_{22})H_{6}+\frac{2}{\omega }(H_{1}+H_{2})H_{4} \label{bautin} \end{equation} \end{proof} \begin{thebibliography}{0} \bibitem{BWolkowicz86} G. J. Butler, G. Wolkowicz; \textit{Predator-mediated competition in the chemostat}, J. Math. Biol., \textbf{24} (1986), 167-191. \bibitem{CLS} M. Cavani, M. Lizana, H. L. Smith; \textit{Stable Periodic Orbits for a Predator-Prey Model with Delay}, J. of Math. Anal. and Appl., \textbf{249} (2000), 324-339. \bibitem{Cushing} J. M. Cushing; \textit{Integrodifferential Equations and Delay Models in Population Dynamics }\textbf{20}, Springer-Verlag, Heidelberg, 1977. \bibitem{Farkasmotion} M. Farkas; Periodic Motion, Springer-Verlag, New York, 1994. \bibitem{holling59} C. S. Holling; \textit{The components of predation as revealed by a study of small-mammal predation of the European pine sawfly}, Canadian Entomologist, \textbf{91} (1959), 293-320. \bibitem{MacD1} N. MacDonald; \textit{Time Lags in Biological Models}, Lecture Notes in Biomathematics \textbf{27}, Springer-Verlag, Heidelberg,1978. \bibitem{smithwaltman95} H. L. Smith, P. Waltman; \textit{The Theory of the Chemostat}, Cambridge University Press, Cambridge, UK, 1995. \bibitem{Wolkowicz1} G. Wolkowicz, H. Xia, S. Ruan; \textit{Competition in the Chemostat: A distributed delay model and its global asymptotic behavior}, SIAM J. Appl. Math. \textbf{57}(1997), 1281-1310. \end{thebibliography} \end{document}