\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}