\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2006(2006), No. 106, pp. 1--11.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2006 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2006/106\hfil An age-structured model] {An age-structured model of cannibalism} \author[R. Ma\v r\'\i k, L. P\v ribylov\'a \hfil EJDE-2006/106\hfilneg] { Robert Ma\v r\'\i k, Lenka P\v ribylov\'a } % in alphabetical order \address{Robert Ma\v r\'\i k \newline Dept. of Mathematics\\ Mendel University \\ Zem\v ed\v elsk\'a 3\\ 613 00 Brno, Czech Republic} \email{marik@mendelu.cz} \address{Lenka P\v ribylov\'a\newline Dept. of Applied Mathematics\\ Masaryk University \\ Jan\'a\v ckovo n\'am. 2a\\ 662 95 Brno, Czech Republic} \email{pribylova@math.muni.cz} \date{} \thanks{Submitted April 3, 2006. Published September 8, 2006.} \thanks{The first author is supported by Grant 201/04/0580 from the Czech Grant Agency.} \subjclass[2000]{34C23, 37G10} \keywords{Cannibalism; predator-prey; Hopf bifurcation; Lyapunov coefficient} \begin{abstract} We investigate the predator-prey model with cannibalism in the predator population, suggested by Magnusson \cite{M} in 1995. We explore the model by a theory of bifurcations, based mainly on the results of Bautin. Among others, we show that the limit cycle appearing in the model due to the Andronov-Hopf bifurcation may be stable or unstable. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{remark}[theorem]{Remark} \newtheorem{corollary}[theorem]{Corollary} \section{Introduction} Magnusson \cite{M} introduced the predator-prey system with age structure and cannibalism in the predator population in the form \begin{equation} \label{eq:system1} \begin{aligned} \dot X&=AY-\mu_aX+\gamma SXY+VCXZ,\\ \dot Y&=\lambda X-AY-\mu_jY-SXY,\\ \dot Z&=LZ-QZ^2-VZX, \end{aligned} \end{equation} where $X$ is the population of adult predators, $Y$ the population of juvenile predators, $Z$ the population of prey, $T$ is the time and the dot denotes derivative with respect to $T$. In model \eqref{eq:system1} the Lotka--Volterra type of interspecific interaction is considered. The parameters $\mu_a$ and $\mu_j$ describe the natural death rate of the adult and juvenile predators, respectively. The constant $\lambda$ is the birth rate of predators and $A$ is the rate at which juvenile predators mature into adults. The term $VXZ$ describes the rate at which adult predators kill the prey and the constant $C\in(0,1)$ is an efficiency of conversion of sources obtained by killing the prey to the increase of the fitness of population of adult predators. In a similar way, the term $SXY$ is the rate at which adult predators kill juvenile predators and the corresponding increase of fitness of adult predators is proportional to this term by the constant $\gamma\in(0,1)$. Remark that the population of prey is subjected to the logistic growth in the absence of predators. Magnusson in \cite{M} used $Q=0$ and proved that for suitable values of parameter $S$ an Andronov-Hopf bifurcation occurs and the increase of the cannibalism rate can destabilize the system. This result is in the same paper extended to small values of $Q$, i.e. for the logistic growth with large carrying capacity, when the competition in the prey population is not significant. Later Kaewmanee and Tang \cite{KT} reexamined model \eqref{eq:system1} (with general $Q$) and obtained similar results. As pointed out already in \cite{H}, it is usually much easier to show that the Andronov-Hopf bifurcation appears in a particular model than to distinguish whether the periodic solution is stable or whether it is unstable and does not describe a state which may really appear in nature. In both papers \cite{KT}, \cite{M}, authors prove that the Andronov-Hopf bifurcation takes place, but the stability of the limit cycle is not analyzed. All authors provide numerical results showing the stable limit cycle and in this sense some of the results are rough. The aim of this paper is to continue in the study of dynamical system \eqref{eq:system1} and to obtain deeper results concerning its stability, topological properties and types of bifurcations (especially stability and uniqueness of limit cycles using Lyapunov coefficients and theory developed in \cite{B,L}). Among others, we focus our attention to more parameters, not only to the parameter $S$. In order to make the computations manageable we consider $Q=0$ (as in \cite{M}), i.e., we suppose no intraspecific competition in the prey population. The paper is organized as follows. The next chapter recalls some facts concerning the theory of bifurcations. The Magnusson's mathematical model is studied in Section $3$. Section $4$ contains some numeric results and a discussion about possible phase portraits for various values of parameters in the model. Among others, we show that an unstable limit cycle exists for suitable values of parameters. \section{Preliminaries} \subsection{Mathematical models} Many problems of natural science can be solved using mathematical models. Some of deterministic models can be described by dynamical systems of autonomous differential equations. These systems contain studied variables and coefficients given by internal attributes of the problem or by the external conditions. Values of this type (constant according to time) are parameters, since the studied variables and behavior of the whole dynamical system depend on them. General mathematical model can be represented by the following system \begin{equation} \label{eq:par_system} \dot x=F(x,\alpha), \end{equation} where $x \in \mathbb{R}^n$ are studied phase variables, $\alpha \in \mathbb{R}^m$ are parameters and $F$ is sufficiently smooth. \subsection{Bifurcation diagrams and normal forms} Changing parameters $\alpha$, the phase portrait of the system \eqref{eq:par_system} changes. There are two possibilities of the change: the phase portrait stays topologically equivalent to the previous one or not. Generally the space of parameters can be divided into structurally stable domains with topologically equivalent phase portraits. This division together with related phase portraits is called a bifurcation diagram of the system \eqref{eq:par_system}. Structurally unstable boundaries correspond with some particular type of bifurcation that can be found by transforming the system \eqref{eq:par_system} to its normal form. Some typical systems (already in the normal form) describing particular types of bifurcations were studied as model systems for many types of bifurcations. \subsection{Andronov-Hopf bifurcation} Now let us turn our attention to the bifurcation which will be proved in the age-structured model of cannibalism, the Andronov-Hopf bifurcation. First of all, recall some basic facts about this bifurcation. This bifurcation is related to an existence of a pair of purely imaginary eigenvalues. By crossing the imaginary axis and changing the real part of eigenvalues from negative to positive value a stable focus changes to an unstable focus. This change is accompanied by a birth of a cycle. The model system for Andronov-Hopf bifurcation \begin{equation}\label{eq:A-H-model} \begin{aligned} \dot x_1&=\mu x_1-x_2-x_1(x_1^2+x_2^2)\\ \dot x_2&=x_1+\mu x_2-x_2(x_1^2+x_2^2) \end{aligned} \end{equation} possesses a unique stationary point $x_1=0=x_2$ for every real $\mu$. The eigenvalues at this stationary points are $\lambda=\mu\pm i$. Introducing complex-valued variable $z=x_1+ix_2$ the system \eqref{eq:A-H-model} takes the form $\dot z=(\mu+i)z-z|z|^2$ which becomes in polar coordinates to a system \begin{equation} \label{eq:A-H-model-pol} \begin{aligned} \dot\rho &=\rho(\mu-\rho^2)\\ \dot\phi &=1. \end{aligned} \end{equation} From here it follows that the stationary point is a stable focus for $\mu<0$ and an unstable focus for $\mu>0$. If $\mu>0$ then a solution $\rho=\sqrt\mu$ presents a limit cycle around the origin. \begin{figure}[ht] \begin{center} \includegraphics[width=1\textwidth]{fig1} % Postscript created from the jpeg file \end{center} \caption{Hopf-Andronov bifurcation. \label{fig1}} \end{figure} An important problem is whether the limit cycle stable or unstable. One of the concepts which enables us to decide which of these two possibilities really arises is the concept of Lyapunov coefficient $l_1$, shortly introduced in the following paragraphs. \begin{lemma} \label{lem2.1} Consider general one-parameter system \begin{equation} \label{eq:par_system-1} \dot x=f(x,\mu), \quad x\in\mathbb{R}^2,\quad \mu\in\mathbb{R} \end{equation} which possesses a stationary point $x=0$ for every sufficiently small $|\mu|$ and let \begin{equation*} \lambda_{1,2}(\mu)=\psi(\mu)\pm i\omega(\mu), \end{equation*} where $\psi(0)=0$, $\omega(0)=\omega_0>0$, be the eigenvalues of this stationary point. There exists an invertible transformation depending on the parameter $\mu$ which converts system \eqref{eq:par_system-1} into the complex form \begin{equation*} \dot z=\Bigl(\psi(\mu)+i\omega(\mu)\Bigr)z+c_1(\mu)z|z|^2+O(|z|^4) \end{equation*} \end{lemma} \begin{remark}[stability of limit cycle in Andronov-Hopf bifurcation] \label{rem:l1} \rm The number \begin{equation*} l_1(\mu)=\frac{\mathop{\rm Re} c_1(\mu)}{\omega(\mu)} -\psi(\mu)\frac{\mathop{\rm Im} c_1(\mu)}{\omega^2(\mu)} \end{equation*} is called the \textit{first Lyapunov coefficient.} The sign of $l_1(0)$ determines stability or instability of the limit cycle. A stable limit cycle appears near the equilibrium in the case $l_1(0)<0$ for $\mu$ close to zero. This process is referred as ``supercritical bifurcation'', since an equilibrium at stationary point is replaced by small oscillations. On the contrary, if $l_1(0)>0$, then ``subcritical bifurcation'' takes place, since the stable equilibrium enclosed by an \textit{unstable} limit cycle changes to an unstable focus. \end{remark} Bautin \cite{B} introduced the terms ``subcritical'' and ``supercritical'' bifurcation and derived formulas for the first Lyapunov coefficient for systems of two and three equations with analytical coefficients. These formulas are too long to repeat them on this place. For some special systems these formulas take simpler forms than presented in \cite{B}. Another approach which enables to distinguish between the stable and unstable limit cycle is presented in \cite{H}. In this book the first Lyapunov coefficient is replaced by the quantity $\mu_2$ which is the first nonvanishing term of an infinite series, see \cite[Chapter 1]{H} for more details. Finally, recall that the period of the limit cycle approaches $\frac{2\pi}{\omega(0)}$ as the parameter approaches the critical value. \subsection{Shoshitaishvili theorem} \begin{lemma} \label{lem2.2} Let for the system \eqref{eq:par_system} be $\alpha \in \mathbb{R}^1$ and $x=0$ be a stationary point for $\alpha = 0$ with $n_0 \neq 0$ purely imaginary eigenvalues. Then there exists a local invariant manifold $W^c(\alpha)$ in the neighbourhood of zero. The manifold is attractive, if all other eigenvalues have positive real parts. The system \eqref{eq:par_system} can be restricted to $W^c(\alpha)$ for small $|\alpha|$ by a local projection of $W^c(\alpha)$ to $T^c$ (generalized eigenspace corresponding to the union of purely imaginary eigenvalues). On this so-called central manifold, the system \eqref{eq:par_system} can be represented in the new coordinates by a system \begin{equation*} \dot u = \Phi(u, \alpha),\quad u \in \mathbb{R}^{n_0}. \end{equation*} \end{lemma} \begin{theorem}[Shoshitaishvili theorem] \label{thm2.1} The system \eqref{eq:par_system} is locally topologically equivalent to a suspension of the central manifold, which can be represented by a system \begin{equation*} \begin{aligned} \dot u &= \Phi(u, \alpha),\\ \dot y &=-y, \\ \dot z &=z, \end{aligned} \end{equation*} where $u \in \mathbb{R}^{n_0}$, $y \in \mathbb{R}^{n_-}$ and $z \in \mathbb{R}^{n_+}$ ($n_\pm$ is the number of eigenvalues with positive or negative real parts). \end{theorem} For the proof of the above theorem, see \cite{SH}. \section{Mathematical model} Let us return our attention to the system \eqref{eq:system1}. First of all, we non-dimensionalise the system and obtain a smaller amount of parameters. Remark that we are interested especially in the interspecific parameters $V$ and $S$ and from this reason we perform the transformation which preserves these parameters. \begin{lemma} \label{lem3.1} Transformation $X=\frac A\gamma x$, $Y=\frac A{\gamma ^2}y$, $Z=\frac{A}{C\gamma}z$, $T=\frac{\gamma}{A}t$ transforms the dynamical system \eqref{eq:system1} with $Q=0$ into \begin{equation} \label{eq:system2} \begin{aligned} \dot x&=y-mx+Sxy+Vxz,\\ \dot y&=nx-oy-Sxy,\\ \dot z&=pz-Vxz, \end{aligned} \end{equation} where the dot, $\dot{\ }$ represents $\frac{d}{dt}$, $m=\frac{\mu_a\gamma}{A}>0$, $n=\frac{\lambda\gamma^2}{A}>0$, $o=\gamma\left(1+\frac{\mu_j}{A}\right)>0$, $p=\frac{L\gamma}A>0$. For $o\neq 1$ and $m\neq n$, System \eqref{eq:system2} possesses three stationary points $S_0=[x_0, y_0, z_0]$, $S_1=[x_1, y_1, z_1]$ and $S_2=[x_2, y_2, z_2]$, where \begin{gather} \label{eq:st-point} x_0=\frac pV,\quad y_0=\frac{np}{oV+pS},\quad z_0=\frac{(mo-n)V+p(m-n)S}{V(oV+pS)}, \\ \label{eq:st-point2} x_1=\frac{om-n}{S(n-m)},\quad y_1=\frac 1S\cdot\frac{om-n}{o-1},\quad z_1=0, \\ \label{eq:st-point3} x_2= y_2= z_2=0. \end{gather} \end{lemma} The statements of this lemma follow immediately. From a practical point of view, the most interesting stationary point is the point $S_0$. The absence of prey in the stationary point $S_1$ is explained in \cite{KT} and \cite{M} as the existence of an alternative food for the predator. This alternative food is not incorporated in our model. However, the existence of an alternative food is not necessary for real interpretation of the system. Particularly, there are lakes in nature with predator fishes only. In such ecosystems, the juvenile predators are the major food for adults. In fact, the set $z=0$ is an invariant set of the system. In this set system \eqref{eq:system2} becomes \begin{equation} \label{eq:system-z0} \begin{aligned} \dot x &=y-mx+Sxy,\\ \dot y &=nx-oy-Sxy \end{aligned} \end{equation} and describes the states with no prey. However, the main aim is to study the steady states in which all populations are present. Therefore we focus our attention to the stationary point $S_0$. \subsection{Stationary point $S_0$} We will focus our attention to the stationary point $S_0$. This point is in the first octant if \begin{equation} \label{eq:podm-pro-1o} (mo-n)V+p(m-n)S>0. \end{equation} This condition covers three mutually different cases \begin{enumerate} \item $mo-n>0$, $m-n>0$ and $V$, $S$ are arbitrary positive numbers \item $mo-n>0$, $m-n<0$ and $V>\frac{S(n-m)p}{mo-n}$ \item $mo-n<0$, $m-n>0$ and $V<\frac{S(n-m)p}{mo-n}$ \end{enumerate} The transformation \begin{align*} x'&=x-x_0,\\ y'&=y-y_0,\\ z'&=z-z_0 \end{align*} moves the stationary point $S_0$ into the origin in new coordinates. In the new variables $[x',y',z']$ system \eqref{eq:system2} reads as \begin{equation} \label{eq:system-po-tr} \begin{aligned} \dot{x'}&=y'-mx'+S(x'y'+x'y_0+y'x_0)+V(x'z'+x'z_0+x_0z'),\\ \dot{y'}&=nx'-oy'-S(x'y'+x'y_0+x_0y'),\\ \dot{z'}&=pz'-V(x'z'+x'z_0+x_0z'), \end{aligned} \end{equation} where dot is the derivative with respect to $t$. The Jacobi matrix of system \eqref{eq:system-po-tr} evaluated at the origin is \begin{equation} \label{eq:jacobi} \begin{aligned} J(0,0,0) &= \begin{pmatrix} -m+Sy_0+Vz_0&1+Sx_0&Vx_0\\ n-Sx_0&-o-Sx_0&0\\ -Vz_0&0&p-Vx_0 \end{pmatrix}\\ &= \begin{pmatrix} -\frac{y_0}{x_0}&1+Sx_0&Vx_0\\ o\frac{y_0}{x_0}&-n\frac{x_0}{y_0}&0\\ -Vz_0&0&0 \end{pmatrix} \end{aligned} \end{equation} and the characteristic polynomial of this matrix is \begin{equation} \label{eq:char-poly} |J(0,0,0)-\nu I|= \begin{vmatrix} -\frac{y_0}{x_0}-\nu&1+Sx_0&Vx_0\\ o\frac{y_0}{x_0}&-n\frac{x_0}{y_0}-\nu&0\\ -Vz_0&0&-\nu \end{vmatrix} =-(\nu^3+A\nu^2+B\nu+C), \end{equation} where \begin{align*} A&=\frac{y_0}{x_0}+n\frac{x_0}{y_0},\\ B&=n-o\frac{y_0}{x_0}(1+Sx_0)+V^2x_0z_0,\\ C&=V^2n\frac{x_0^2z_0}{y_0}. \end{align*} Consider the most important case, when the point $S_0$ is in the interior of the first octant. In this case clearly $A>0$ and $C>0$. Let us examine the expression $AB-C$. An evaluation shows that \begin{equation} \label{eq:AB-C-0} \begin{aligned} AB-C=&\Big(\frac{y_0}{x_0}+n\frac{x_0}{y_0}\Big) \Big(n-o\frac{y_0}{x_0}(1+Sx_0)+V^2x_0z_0\Big) -V^2n\frac{x_0^2z_0}{y_0}\\ =&\frac{y_0}{x_0}\Bigl(n-o\frac{y_0}{x_0}(1+Sx_0)+V^2x_0z_0\Bigr) +n^2\frac{x_0}{y_0}-no(1+Sx_0)\\ =&\big[1+n\big(\frac{x_0}{y_0}\big)^2\big]\big[ n-o\frac{y_0}{x_0}-Soy_0\big]\frac {y_0}{x_0}+V^2y_0z_0\\ =&\big[1+n\big(\frac{x_0}{y_0}\big)^2\big]\bigl[ Sy_0-Soy_0\bigr]\frac {y_0}{x_0}+V^2y_0z_0\\ =&\big[1+n\big(\frac{x_0}{y_0}\big)^2\big]\frac {y^2_0}{x_0} S(1-o)+V^2y_0z_0\\ =&\frac{1}{x_0}\Bigl[(y_0^2+nx_0^2)S(1-o)+Vpy_0z_0\Bigr]. \end{aligned} \end{equation} Hence if $o<1$, then $AB-C$ is positive. In this case all real parts of eigenvalues are negative and the stationary point $S_0$ is stable. The only possibility which allows the point $S_0$ become to be unstable is for $o>1$. Under this condition the point $S_0$ lies in the first octant if and only if either \begin{equation*} mo>m>n, \end{equation*} or \begin{equation*} mo>n>m \quad\text{and}\quad V>Sp\frac{n-m}{mo-n} . \end{equation*} \subsection{Case $o>1$} We have to write the expression $AB-C$ in terms of parameters of system \eqref{eq:system2} to obtain correct conclusions concerning the influence of parameters on stability. (Remember that all $x_0$, $y_0$ and $z_0$ depend on the parameters and hence \eqref{eq:AB-C-0} cannot be used to obtain correct conclusions. This step seems to be omitted in \cite{KT}.) A direct computation shows that the expression $AB-C$ can be written in the form of the product of a positive factor and the three-degree polynomial $P(\cdot)$ in variable $\frac VS$ as follows: \begin{equation}\label{eq:ABCP} AB-C=\frac{pnS^3}{V(oV+pS)^2} P\Bigl(\frac VS\Bigr), \end{equation} where \begin{equation} \label{eq:Ptau} P(\tau)= (mo-n)\tau^3 +\bigl[p(m-n)+(1-o)(n+o^2)\bigr)]\tau^2 +2op(1-o)\tau+(1-o)p^2\,. \end{equation} According to the Descarte's rule of signs, the polynomial $P(\tau)$ possesses a unique positive zero. An evaluation shows that $P\bigl(p\frac{n-m}{mo-n}\bigr)$ is negative. Really \begin{equation} \label{eq:Pneg} \begin{aligned} P\Bigl(p\frac{n-m}{mo-n}\Bigr)=&p^3\frac{(n-m)^3}{(mo-n)^2} +\Bigl[p(m-n)+(1-o)(n+o^2)\Bigr] \frac{(n-m)^2}{(mo-n)^2}p^2\\ &+ 2op(1-o)p\frac{n-m}{mo-n}+(1-o)p^2\\ =&(1-o)p^2\Bigl[ (n+o^2)\frac{(n-m)^2}{(mo-n)^2}+2o\frac{n-m}{mo-n}+1 \Bigr]\\ =&(1-o)p^2\Bigl[ n\frac{(n-m)^2}{(mo-n)^2}+ \Bigl(o\frac{n-m}{mo-n}+1\Bigr)^2 \Bigr]. \end{aligned} \end{equation} We summarize the above computations in the following theorem. \begin{theorem} \label{th:stab} Let either $ mo>m>n$ or $ mo>n>m$ hold. Denote by $\tau_k$ the unique positive zero of the polynomial \eqref{eq:Ptau}. \begin{enumerate} \item If $\frac VS>\tau_k$, then all real parts of eigenvalues of Jacobi matrix at $S_0$ are negative and $S_0$ is a stable singular point in the first octant. \item If either \begin{equation*} \text{$mo>m>n$\quad and\quad $\frac VS<\tau_k$} \end{equation*} or \begin{equation*} \text{$mo>n>m$\quad and\quad $p\frac{n-m}{mo-n}<\frac VS<\tau_k$,} \end{equation*} then at least one of the eigenvalues of Jacobi matrix at $S_0$ has a positive real part and $S_0$ is an unstable singular point in the first octant. \item If $mo>n>m$ and $\frac VS1$. From equality \eqref{eq:Pneg} it follows that $\tau_k>p\frac{n-m}{mo-n}$. Now Theorem follows from \eqref{eq:ABCP}, \eqref{eq:Pneg} and from the well known Ruth--Hurwitz criterion. \end{proof} \begin{remark} \label{rmk3.1} \rm Note that if $o<1$, then $AB-C$ is always positive and the point $S_0$ is stable. No bifurcation occurs if $o<1$. If $n>mo>m$, then $z_0<0$ for all positive values of $V$ and $S$ and the point $z_0$ lies outside the first octant. If $mo>n>m$, then $\frac VS>p\frac{n-m}{mo-n}$ is a necessary and sufficient condition for $\min(x_0,y_0,z_0)>0$. \end{remark} \section{Andronov-Hopf bifurcation of the system \eqref{eq:system2}} In this section we consider system \eqref{eq:system2} with positive parameters $m,n,p,S,V$ and $o>1$ (according to the previous remark). We make a natural assumption $\min(x_0,y_0,z_0)>0$. \begin{corollary} \label{coro4.1} Let $\tau_k$ be the zero of the polynomial \eqref{eq:Ptau} for arbitrary fixed positive parameters $m,n,o,p,S$. Then $V_k=\tau_kS$ is the critical value of Andronov-Hopf bifurcation of the system \eqref{eq:system2} for the stationary point $S_0$. \end{corollary} \begin{proof} Follows immediately from the Theorem \eqref{th:stab}, since the characteristic polynomial \eqref{eq:char-poly} has two purely imaginary eigenvalues, if $AB-C=0$ (notice that $A>0$ and $C>0$), that is while $P(\tau)=0$. \end{proof} \begin{remark} \label{rmk4.1} \rm Both the supercritical and the subcritical bifurcation types occur in the system \eqref{eq:system2}. The type of Andronov-Hopf bifurcation is determined by the first Lyapunov coefficient $l_1$, as Remark \ref{rem:l1} shows. This coefficient can be calculated numerically for given critical parameters. \end{remark} For example, let $m=7$, $n=3$, $o=3$, $p=5$ and $S=1$. The polynomial \eqref{eq:Ptau} has a unique zero $\tau_k \doteq 2.251$ and $V_k = \tau_k$ is a critical value of Andronov-Hopf bifurcation. Following \cite[Chapter III]{B}, we transfer the stationary point $S_0$ to the origin and find that $l_1 \doteq -1.683\pi$. In this case the supercritical bifurcation takes place and the existence of a stable limit cycle in the neighborhood of an unstable focus for $VV_k$ near $V_k$. Figure 3 shows parts of two trajectories for $V=5.5$. One trajectory (with initial conditions $[x(0)= 2,y(0) =1,z(0)=2]$) is converging to the stable focus $S_0$ on the central manifold and the other is a stable separatrix (with initial conditions $[x(0)= 2,y(0) =2,z(0)=0]$, $z=0$ is an invariant set) of the saddle point in the origin. These trajectories are very close to each other and they are converging to the central manifold, so we can suppose that the unstable limit cycle is somewhere in between on the central manifold (unlike the planar phase space, in the 3d phase space it is very hard to find the unstable limit cycle by drawing trajectories). \begin{figure}[ht] \begin{center} \includegraphics[width=1\textwidth]{fig3} \end{center} \caption{An unstable limit cycle. \label{fig3}} \end{figure} In view of the above computations, we can see that another type of bifurcation must take place here. Changing parameter $m$ from 7 to 4 and counting the critical value $V_k$ we got two qualitatively different phase portraits. You can verify that for $m \doteq 4.254$, the critical value $V_k \doteq 3.792$ and $l_1 = 0$. The figure 4 shows the trajectories near the center (initial conditions $[x(0)= 2,y(0) =2,z(0)=2], \, [x(0)= 2,y(0) =-2,z(0)=2] \, [x(0)= 2,y(0) =-2,z(0)=1], \, [x(0)= 2,y(0) =1,z(0)=2], \, [x(0)= 2,y(0) =2,z(0)=0]$ ). \begin{figure}[ht] \begin{center} \includegraphics[width=1\textwidth]{fig4} \end{center} \caption{Phase portrait of the center. \label{fig4}} \end{figure} This is so called generalized Hopf (Bautin) bifurcation and more limit cycles can arise in the neighbourhood of the stationary point $S_0$. \subsection*{Conclusion} We investigated the system introduced by Magnusson \cite{M} as a model of cannibalism. We derived the conditions under which the Andronov-Hopf bifurcation takes place. These conditions are derived in terms of the parameters of the original model and among others, we provide an equation which describes the critical value for which the bifurcation occurs. Using theory of Bautin we proved that both subcritical and supercritical bifurcations may take place in this model and hence the limits cycle enclosing the stationary point need not to be stable. This phenomenon was not mentioned in any of previous works on this system. \begin{thebibliography}{00} \bibitem{B} N. N. Bautin: {\it Povedenie dinamicheskih sistem vblizi granic oblasti ustoychivosti}, OGIZ GOSTEXIZDAT, Leningrad-Moskva, 1949 \bibitem{H} B. D. Hassard, N. D. Kazarinoff, Y.--H. Wan: {\it Theory and applications of Hopf bifurcation}, Cam. Univ. Press (1981), russian transl. Moskva, Mir (1985). \bibitem{KT} C. Kaewmanee, I. M. Tang: Cannibalism in an age-structured predator-prey system, Ecol. Mod. 167 (2003), 213--220. \bibitem{L} Lyapunov: {\it Povedenie dinamicheskih sistem vblizi granic oblasti ustoychivosti}, OGIZ GOSTEXIZDAT, Leningrad-Moskva, 1949 \bibitem{M} K. G. Magnusson: Destabilizing efect on cannibalism on a structures predator-prey system, Math. Biosci. 155 (1999), 61--75. \bibitem{SH} A. N. Shoshitaishvili; {\it Bifurkacii topologicheskovo tipa vektornovo polya vblizi osoboy tochki}, Trudy seminara im. I.G. Petrovskovo, No. 1 (1975) 279--309. \end{thebibliography} \end{document}