text:
\begin{document}
\maketitle

\begin{abstract}
This paper, as a complement to the works by Hsu et al [SIAM. J. Appl. Math. 55 (1995)] and Huang et al [J. Differential Equations 257 (2014)], aims to examine the Hopf bifurcation and global
dynamics of a predator-prey model of Leslie type with generalized Holling type III functional response for the two cases: (A) when it has a unique non-degenerate positive equilibrium; (B)
when it has three distinct positive equilibria.
For each case, the type and stability of each equilibrium, Hopf bifurcation at each weak focus, and the number and distribution of limit cycles in the first quadrant are studied.
It is shown that every equilibrium is not a center.
For the case (A), $i$ limit cycle(s) can appear near the unique positive equilibrium, $i=1, \cdots, 4$.
For $i=3$ or $4$, the model has two stable limit cycles, which gives a positive answer to the open problem proposed by Coleman [Differential equations model,1983]: finding at least two ecologically stable cycles in any natural (or laboratory) predator-prey system.
For the case (B), one positive equilibrium is a saddle and the others are both anti-saddle.
If one of the two anti-saddles is weak focus and the other is not, then the order of the weak focus is at most $3$.
If both anti-saddles are weak foci, then they are unstable weak foci of order one.
Moreover, one limit cycle can bifurcate from each of them simultaneously.
Numerical simulations demonstrate  that there is also a big stable limit cycle enclosing these two limit cycles.
Our results indicate that the maximum number of limit cycles in the model of this kind is at least $4$, which improves the preceding results that this number is at least $2$.
\end{abstract}
\maketitle
{\bf Keywords}: two ecologically stable cycles, Hopf bifurcation, predator-prey system, generalized Holling type III functional response\\
\rule[1.5mm]{\textwidth}{0.1pt}

\section{Introduction and statement of the main results}

In this work we consider the predator-prey system of Leslie type with generalized Holling type III or sigmoidal functional response \cite{Huang2014},
\begin{equation}\label{model}
\left\{
\begin{aligned}
&\dot{x}=r x \left(1-\dfrac{x}{K}\right)-\dfrac{m x^2 y}{a x^2+b x+1}, \\
&\dot{y}=s y \left(1-\dfrac{y}{h x}\right)
\end{aligned}
\right.
\end{equation}
with $(x, y)\in \mathcal{A}=\{(x, y)\left|x>0, y\geq 0 \right.\}$ and $r,\ K,\ m,\ a,\ s,\ h>0$, $b>-2\sqrt{a}$.

In system \eqref{model} $x(t)$ and $y(t)$ represent the population densities of the prey and predator at time $t>0$, respectively; the effect of the predation is given by the function
\begin{equation*}
p(x)=\frac{m x^2}{a x^2+b x+1},
\end{equation*}
and is called  the generalized Holling type III or sigmoidal functional response \cite{Bazykin1998}. The parameters  $r$ and $s$ are the intrinsic growth rates or biotic potential of the prey and
predator, respectively,  $K$ is the prey environment carrying capacity, and   $h$ is a measure of the food quality of the prey for conversion into
predator births.

Here the predator growth equation is of Leslie form originated by Leslie \cite{Leslie1948}, but the conventional environmental carrying capacity $h x$ is proportional to prey abundance $x$ \cite{May1973}.

The system \eqref{model} has been investigated by Hsu et al \cite{Hsu1995} and Huang et al \cite{Huang2014}. It has been conjectured that for predator-prey systems with a unique positive
equilibrium, local and global stability are equivalent \cite{Saez1999}.
When $b>0$, Hsu et al \cite{Hsu1995} prove that it is true under certain conditions.
When $b>-2\sqrt{a}$ (so that $a x^2+b x+1>0$ and hence $p(x)>0$ for all $x>0$), Huang et al \cite{Huang2014} show that it is not true by analysing the nonlinear dynamics
of system \eqref{model} when it has at least one degenerate positive equilibrium.

The predator-prey systems have been widely studied, such as investigating the
stability and bifurcations, proving the global stability of the unique positive equilibrium, and studying the uniqueness or nonexistence of limit cycles, see in \cite{May1973,Bazykin1998,Hsu1995} and the references therein.
However, it's not an easy task to study the number and distribution of limit cycles in a given predation model. This problem is related to Hilbert's 16th Problem \cite{Gaiko}, and it is a question that has remained unsolved for the predation model \cite{Eduardo2011}.

For some two-dimensional predator-prey systems, many authors proved that there can exist at least two limit cycles, see  \cite{Kuang1988,Saez1999,Ruan2001,Wrzosek1990,Zhu2002}, etc.
In \cite{Aguirre2009}, Aguirre et al showed that there can exist three limit cycles including only one stable cycle.
In \cite{Lu2001,Lloyd1996} the existence of some Kolmogorov type systems with at least two stable limit cycles surrounding the singularity in the positive quadrant is given. We note that the systems are not predator-prey systems. So far, whether there are more than $3$ limit cycles in predator-prey systems is still unknown.
The study of the number and distribution of limit cycles in a given predation model may be the most difficult part. This is principally because these systems are generally not polynomial differential equations, and the expressions of some positive equilibria of such systems are so complicated and even not formulated explicitly that one can not study them.
In this paper, we will propose some available methods for studying the existence, stability, number and distribution of limit cycles in a given predator-prey systems.

This paper aims to study the system \eqref{model} when it has no degenerate positive equilibrium. If $E(x_{*}, y_{*})$ is a positive equilibrium of system \eqref{model}, then $y_{*}=hx_{*}$ and $x_{*}$ is a root of the equation
\begin{equation}\label{dai1.2}
a x^3+\left(\frac{K m h}{r}+b-K a\right) x^2+(1-K b) x-K=0
\end{equation}
in the interval $(0, K)$.
We denote the determinant and trace of the Jacobian matrix of system \eqref{model} at $E$ by $\textrm{det}\left(J(E)\right)$ and $\textrm{tr}\left(J(E)\right)$, respectively.
If $\textrm{det}\left(J(E)\right)\neq0$, then  $E$ is called an elementary equilibrium,  otherwise  it is a degenerate equilibrium.
Specially, $E$ is called a hyperbolic saddle  if $\textrm{det}\left(J(E)\right)<0$ and called \emph{center or focus type} if $\textrm{det}\left(J(E)\right)>0$ and $\textrm{tr}\left(J(E)\right)=0$, respectively.

The number of positive  equilibria of
system \eqref{model} is determined by the number of roots of Eq. \eqref{dai1.2} in the interval $(0, K)$.
Note that  Eq. \eqref{dai1.2} can have one, two or
three positive roots in the interval $(0, K)$. Correspondingly, system \eqref{model} can have one, two, or three positive equilibria.
For system \eqref{model}, it can be inferred from Hsu and Huang in \cite{Hsu1995} that if $x_{*}$ is a multiple root of Eq. \eqref{dai1.2} in the interval $(0, K)$, then $E(x_{*}, y_{*})$ must be a degenerate equilibrium of system \eqref{model}.
Let
\begin{equation*}
\Delta=-4A^3+\left[27 a^2 K+9 a (1-K b) \left(\frac{K m h}{r}+b-K a\right)-
2 \left(\frac{K m h}{r}+b-K a \right)^3\right]^2,
\end{equation*}
where $A=\left(K m h/r+b-K a\right)^2 +3a(Kb-1)$.
Then the case $\Delta=0$ corresponds to the case that system has a degenerate positive equilibrium, which has been studied by Huang et al \cite{Huang2014}.
Hence we only need to consider the case $\Delta\neq 0$, i.e., system \eqref{model} has either a unique non-degenerate positive equilibrium or three distinct positive equilibria.

First, we give the main results about the system with a unique non-degenerate positive equilibrium, which are stated by the following two theorems.
\begin{theo}\label{T3.2}
Suppose $\Delta>0$, then system \eqref{model} has a unique positive equilibrium, which is an elementary and anti-saddle equilibrium.
If the unique positive equilibrium is unstable, then there exists at least one stable limit cycle in the first quadrant.
\end{theo}
\begin{theo}\label{T4.1}
If  $\Delta>0$ and the unique positive equilibrium is center or focus type, then
\begin{itemize}
\item[(1)] it is not a center;
\item[(2)] it is a weak focus of order at most $4$ and there exists a unique class of  parameter values  such that its order is $4$;
\item[(3)] there exist some parameter values  such that system \eqref{model} has $i$ small limit cycles around it, for each $i=1, \cdots, 4$.
\end{itemize}
\end{theo}

It's worth to point out that the unique class of parameter values in the statement (2) of Theorem \ref{T4.1} satisfies $b>0$, see the proof of this theorem in Section $3$. Therefore, for the case $b\leq 0$, if the unique non-degenerate positive equilibrium of system \eqref{model} is center or focus type, then it is a weak focus of order at most $3$, implying at most $3$ limit cycles can bifurcate from it.

The phenomenon that one limit cycle can appear around the unique positive equilibrium was also proved by Hsu et al \cite{Hsu1995}  and  Huang et al \cite{Huang2014}.
The phenomenon that two limit cycles can appear near the unique positive equilibrium
was also observed by Huang et al \cite{Huang2014} through subcritical Hopf bifurcation and numerical simulations.
However, the phenomena that three or four limit cycles can appear near the unique positive equilibrium, observed in the present paper, have not been found by other authors.
Furthermore, the coexistence of four limit cycles including two stable cycles or three limit cycles including two stable cycles, which is not yet reported in the predation model, gives a positive answer to the open problem proposed by Coleman \cite{Coleman1983}, suggesting: ``Find a predator-prey or other interacting system in nature, or construct one in the laboratory, with at least two ecologically stable cycles''.

Second, we will consider the system with three different positive equilibria. The main results are given by the following two theorems.
\begin{theo}\label{T4.3}
Suppose  $\Delta<0$ and $-2\sqrt{a}<b<K a-K m h/r$, then system \eqref{model} has  three different positive equilibria, of which one is a saddle and the others are anti-saddles. If one of the two positive anti-saddles is center or focus type and the other is not, then
\begin{itemize}
\item[(1)] the center or focus type equilibrium is not a center. It is a weak focus of order at most $3$;
\item[(2)] two limit cycles can bifurcate from the weak focus.
\end{itemize}
\end{theo}

\begin{rema}\label{rem3}
When system has three distinct positive equilibria, two of them are anti-saddles.
The phenomenon that one limit cycle can bifurcate from one of the two positive anti-saddles was also observed by Huang et al \cite{Huang2014}.
However, the phenomenon that two limit cycles can bifurcate from one of the two positive anti-saddles, observed in the present paper, has
not been found by other authors. Unfortunately, we can't prove whether there exist three limit cycles surrounding one of the two positive anti-saddles.
\end{rema}
\begin{figure}\label{Fig1.1}
\begin{center}
{\includegraphics[height=7cm,width=9cm]{Fig1.1-2.eps}}
\end{center}
{\small {\bf Fig. 1.1.} Coexistence of one big stable limit cycle enclosing two unstable limit cycles in system \eqref{model}. Here, we fix $r=1$, $K=10$, $m=0.54$, $a=1.25$, $b=-1.65$, $h=1$ and $s=0.4$.}
\end{figure}

\begin{theo}\label{T4.4}
Suppose $\Delta<0$ and $-2\sqrt{a}<b<K a-K m h/r$, i.e., system \eqref{model} has one  positive saddle and two positive anti-saddles. If the two positive anti-saddles are center or focus type, then
\begin{itemize}
\item[(1)] they are unstable weak foci of order one;
\item[(2)] one limit cycle can bifurcate from each of them simultaneously.
\end{itemize}
\end{theo}

When system has three distinct positive equilibria, of which two are anti-saddles, one limit cycle can bifurcate from each of the two positive anti-saddles simultaneously (see Fig. 1.1). This has not been observed  by Hsu et al \cite{Hsu1995} and
Huang et al. \cite{Huang2014}. Numerical simulations show that there is also a big stable limit cycle enclosing these two limit cycles (see Fig. 1.1), which is a new phenomenon observed in the present paper. Unfortunately, we can't prove it.

The rest of this paper is organized as follows. Section $2$ is devoted to some  preliminary results including model reduction, analysis of equilibria of an equivalent polynomial differential system of \eqref{model} and  computation of Lyapunov constants at positive equilibrium. In Section $3$,  we give the proof of our main results.
The paper ends with a discussion.

\section{Some preliminary results}

The preliminary results provided in this section are helpful for the proof of our results.

\subsection{Model reduction}

In order to calculate the Liapunov constants at the nonhyperbolic focus and
study the Hopf bifurcation in the parameter space in a simpler way, it is necessary to reduce system \eqref{model} to a polynomial differential system.

Without loss of generality, assume that $E(x_{*}, y_{*})$ is an arbitrary positive equilibrium of system \eqref{model}.
Let us consider the change of variables
\begin{equation}\label{change1}
u=\dfrac{x}{x_{*}},\,\,\,\, v=\dfrac{y}{y_{*}},\,\,\,\, \tau=r t,
\end{equation}
and the parameter transformation given by $\zeta: \mathbb{R}_+^6\times \mathbb{R}\rightarrow \mathbb{R}_+^6\times \mathbb{R}$,
\begin{equation*}
\zeta (r, K, m, a, s, h, b)=\left(r, \frac{K}{x_{*}}, \frac{m x_{*}y_{*}}{r}, a x_{*}^2, \frac{s}{r}, \frac{h x_{*}}{y_{*}}, b x_{*}\right):=(r, K, m, a, s, h, b)
\end{equation*}
with Jacobian $\textrm{det}D \zeta (r, K, m, a, s, h, b)=\frac{x_{*}^4}{r^2}>0$, which implies that $\zeta$ is invertible.
For simplicity, let us rename the new variables $u\rightarrow x$, $v\rightarrow y$,  $\tau\rightarrow t$. Then in the new variables, system \eqref{model} is reduced to the following system
\begin{equation}\label{norm1}
\left\{
\begin{aligned}
&\dot{x}=x\left(1-\dfrac{x}{K}\right)-\dfrac{m x^2 y}{a x^2+b x+1}, \\
&\dot{y}=s y\left(1-\dfrac{y}{h x}\right).
\end{aligned}
\right.
\end{equation}
Noting that system \eqref{norm1} has an equilibrium at $E^{*}(1, 1)$, we have $m=(a+b+1) (1-1/K)$ and  $h=1$. Furthermore, $K>x_{*}$ becomes $K>1$ and $a>0$, $b>-2\sqrt{a}$,
$s>0$ remain the same, respectively, under the projection $\zeta$.
That is to say, system \eqref{norm1} can be written as
\begin{equation}\label{norm2}
\left\{
\begin{aligned}
&\dot{x}=x\left(1-\dfrac{x}{K}\right)-
\dfrac{(a+b+1)(K-1)x^2y}
{K(ax^2+bx+1)}, \\
&\dot{y}=sy\left(1-\dfrac{y}{x}\right)
\end{aligned}
\right.
\end{equation}
with $(x, y)\in \mathcal{A}$ and the new vector of parameters $\eta=(s, K, a, b)\in \mathcal{D}_0$ is given by the natural projection, where
\begin{equation*}
\mathcal{D}_0=\left\{\eta\in \mathbb{R}_+^3\times \mathbb{R}\left|s>0,\, K>1, \, a>0,\,  b>-2\sqrt{a}\right.\right\}.
\end{equation*}
Since the transformation \eqref{change1} is a linear sign-reserving transformation,  system \eqref{norm2} and system \eqref{model} have the same qualitative property.
By scaling
\begin{equation*}
d\tau=\dfrac{dt}{K x (a x^2+b x+1)},
\end{equation*}
(we will still use $t$ to denote $\tau$ for ease of notation) system \eqref{norm2} can be transformed into  the following quintic polynomial differential system
\begin{equation}\label{norm}
\left\{
\begin{aligned}
&\dot{x}=(ax^2+bx+1)\left(K-x\right)x^2-(K-1)(a+b+1)x^3y, \\
&\dot{y}=K s y\left(x-y\right)(ax^2+bx+1),
\end{aligned}
\right.
\end{equation}
with $(x, y)\in \bar{\mathcal{A}}=\{(x, y)\left|x, y\geq 0 \right.\}$ and $\eta\in \mathcal{D}_0$.
Notice that system \eqref{norm} has the same topological structure as system \eqref{norm2} in $\mathcal{A}$ because $K x (ax^2+bx+1)>0$ for all $x>0$.
In other words, system \eqref{norm} is equivalent to system \eqref{model} in $\mathcal{A}$. Hence we only need to study the polynomial system \eqref{norm} in the region $\bar{\mathcal{A}}$ with $\eta\in \mathcal{D}_0$.

\subsection{Analysis of equilibria  of  system \eqref{norm}}

Obviously, the $x$-axis, $y$-axis and the interior of $\bar{\mathcal{A}}$ are all invariant under system \eqref{norm}.
It's  standard  to  show  that  all solutions of  \eqref{norm} with positive initial values are positive and bounded, and  will eventually tend into the region $\Omega=\{(x(t),  y(t)): 0 \leq x(t)\leq K,\, 0\leq y(t)\leq K\}$. Therefore, $\Omega$ is a positive invariant set of system \eqref{norm} and the limit cycle of  \eqref{norm}, if it exists, must be inside $\Omega$.

Notice that system \eqref{norm} always has a boundary equilibrium $E_0=(K, 0)$ for all parameters, which is always a hyperbolic saddle and divides the positive $x$-axis into two parts which are two stable manifolds of $E_0$. Furthermore, there exists a unique unstable manifold of $E_0$ in the interior of  $\bar{\mathcal{A}}$.

The origin $O(0, 0)$ is an isolated critical point of higher order of system \eqref{norm}.
By introducing the polar coordinates $x=r\cos\theta$, $y=r\sin\theta$, we get the characteristic equation of system \eqref{norm} as follows.
\begin{equation}\label{dai3.1}
\begin{aligned}
G(\theta)=K \sin\theta\cos\theta[(s-1)\cos\theta-s\sin\theta]=0.
\end{aligned}
\end{equation}
Noting that $G(\theta^{*})=0$ is a necessary condition for $\theta=\theta^{*}$ to be a characteristic direction \cite{ZhangZF},  there are at most $3$ possible directions $\theta=0,  \pi/2, \theta_3$, where
$\theta_3=\arctan(1-1/s)$, along which an orbit of system \eqref{norm} may approach the origin in the first quadrant.
Let  $S^+_{\delta} (O)=\{(r, \theta):0<r<\delta,\
0\leq\theta\leq\pi/2\}$, where $0<\delta\ll 1$.
By Theorems 3.4, 3.7 and 3.8 in \cite{ZhangZF}, we get the following results.
\begin{lemm}\label{T3.1}  For system \eqref{norm}, the origin $O$ is an isolated critical point.
\begin{itemize}
\item[(1)] If  $0<s<1$, then
\begin{itemize}
\item[(a)] the positive $x$-axis to the left of $E_0$ is a unique orbit of system \eqref{norm} tending to $O$ along $\theta=0$ as  $t\rightarrow-\infty$;
\item[(b)] the positive $y$-axis is a unique orbit of system \eqref{norm} tending to $O$ along $\theta=\pi/2$ as $t\rightarrow+\infty$.
\end{itemize}
\item[(2)] If $s=1$, then
\begin{itemize}
\item[(a)] there is an infinite number of orbits of system \eqref{norm} in $S^+_{\delta} (O)$ tending to $O$ along $\theta=0$ as $t\rightarrow-\infty$;
\item[(b)]  the positive $y$-axis is a unique orbit of system \eqref{norm} tending to $O$ along $\theta=\pi/2$ as $t\rightarrow+\infty$.
\end{itemize}
\item[(3)] If $s>1$, then
\begin{itemize}
\item[(a)]  there is an infinite number of orbits of system \eqref{norm} in $S^+_{\delta} (O)$ tending to $O$ along  $\theta=0$ as $t\rightarrow-\infty$;
\item[(b)]  the positive $y$-axis is a unique orbit of system \eqref{norm} tending to $O$ along $\theta=\pi/2$ as
  $t\rightarrow+\infty$;
\item[(c)] there is a unique orbit of system \eqref{norm} tending to $O$ along $\theta=\theta_3$ as $t\rightarrow-\infty$. And this orbit is a separatrix that divides $S^+_{\delta} (O)$
    into two parts: one is a hyperbolic sector, the other is a repulsing parabolic sector.
\end{itemize}
\end{itemize}

The phase portrait of system \eqref{norm} near $O$ is shown in Fig. 2.1 for all cases.
\end{lemm}
\begin{figure}\begin{center}
{\includegraphics[height=5cm,width=13cm]{Fig3.1.eps}}
\end{center}
{\small {\bf Fig. 2.1.} The sketch graph of the vector field on the boundary of $\Omega$ and topological structure of the orbits of system \eqref{norm} near the origin. }
\end{figure}

Now we are going to study the positive (i.e., interior) equilibria of system \eqref{norm}.
Assume that $\tilde{E}(\tilde{x}, \tilde{y})$ is a positive equilibrium of system \eqref{norm}, then $\tilde{y}=\tilde{x}$ and $\tilde{x}$ is a positive root of the equation
\begin{equation}\label{dai2.13}
\begin{aligned}(x-1)(a x^2+(K b+K-1)x+K)=0
\end{aligned}
\end{equation}
on the interval $(0, K)$.
The Jacobian matrix  of system \eqref{norm} at $E^*(1, 1)$ is
\begin{equation*}
\textrm{J}(E^*)=\left(\begin{array}{cc}
K a-K-2 a-b & -(K-1) (a+b+1)\\
K s (a+b+1)  & -K s (a+b+1)
\end{array}\right).
\end{equation*}
Then the characteristic polynomial of $\textrm{J}(E^*)$ is
\begin{equation*}
\Theta(\lambda)=\lambda^2-\textrm{tr}\left(J(E^{*})\right)\lambda+\textrm{det}\left(J(E^{*})\right),
\end{equation*}
where
\begin{equation}\label{TR}
\begin{aligned}
&\textrm{tr}\left(J(E^{*})\right)=
K a-K-2 a-b-K s (a+b+1),\\
&\textrm{det}\left(J(E^{*})\right)=
K s (a+b+1) (K b+2 K+a-1).
\end{aligned}
\end{equation}
Hence, $E^*$ is a saddle if $\textrm{det}\left(J(E^{*})\right)<0$, is a degenerate equilibrium if $\textrm{det}(J(E^{*})) =0$, and is an elementary equilibrium if
$\textrm{det}\left(J(E^{*})\right)>0$, respectively. Specially, $E^*$ is \emph{center or focus type} if $\textrm{det}\left(J(E^{*})\right)>0$ and $\textrm{tr}\left(J(E^{*})\right)=0$.

Now we assume that $\textrm{det}\left(J(E^{*})\right)>0$, i.e., $K b+2 K+a-1>0$. Then the stability of $E^*$ is determined by the sign of $\textrm{tr}\left(J(E^{*})\right)$.
This gives that $E^*$ is locally asymptotically stable (or unstable) if $\textrm{tr}\left(J(E^{*})\right)<$ (or $>$) $0$.
Specially, $E^{*}$ is locally asymptotically stable if $K a-K-2 a-b\leq 0$.

Furthermore, a straightforward calculation shows that
\begin{equation}\label{dai2.10}
\left(\textrm{tr}\left(J(E^{*})\right)\right)^2-4 \textrm{det}\left(J(E^{*})\right)=\mu_2 {s}^{2}+ \mu_1 s+ \mu_0 \triangleq \psi(s),
\end{equation}
where $\mu_2={K}^{2} (a+b+1)^{2}>0$, $\mu_1=-2 K ( a+b+1)
 ( K a+2 K b+3 K-b-2 )$, $\mu_0=( K a-K-2 a-b )^{2}\geq 0$.
Note that $\tilde{\Delta}=\mu_1^2-4 \mu_2 \mu_0=16\,{K}^{2}(a+b+1)^{3} (K-1) (K b+2 K+a-1)>0$. From \eqref{dai2.10}, if $K a+2 K b+3 K-b-2\leq 0$, then $\psi(s)>0$ and hence $E^{*}$ is a node.

Let  $K a+2 K b+3 K-b-2>0$.
If $K a-K-2 a-b=0$, then $\psi(s)$ has a unique positive zero at $s_1$, where $s_1=-\mu_1/\mu_2$. Furthermore, $\psi(s)\geq 0$  for $s\in[s_1, +\infty)$ and
$\psi(s)<0$ for $s\in(0, s_1)$, respectively.
Notice that $E^{*}$ is locally asymptotically stable for this case. Hence, if $s>s_1$, then $E^*$ is a locally asymptotically stable node. Otherwise, it is a locally asymptotically stable focus.
Let $K a-K-2 a-b\neq 0$, then $\psi(s)$ has two positive zeros at $s_2$ and $s_3$ with $0<s_2<s_3$, where
\begin{align*}
&s_2=\frac{K a+2 K b+3 K-b-2-2 \sqrt{( K-1)( a+b+1)(K b+2 K+a-1)}}{ K (a+b+1)},\\
&s_3=\frac{K a+2 K b+3 K-b-2+2 \sqrt{( K-1)( a+b+1)(K b+2 K+a-1)}}{ K (a+b+1)}.
\end{align*}
Furthermore, $\psi(s)\geq 0$ for $s\in(0, s_2]\cup[s_3, +\infty)$ and $\psi(s)<0$ for $s\in(s_2, s_3)$, respectively. Thus $E^*$ is a node for $s\in(0, s_2)\cup(s_3, +\infty)$ and is a focus or center for  $s\in(s_2, s_3)$.

We note that, if $K a-K-2 a-b>0$, then $K a+2 K b+3 K-b-2=2 (K b+2 K+a-1)+(K a-K-2 a-b)>0$.
For this case, $\textrm{tr}\left(J(E^{*})\right)$ as a function of $s$ has a positive zero at $s^{*}$,  where $s^{*}=(K a-K-2a-b)/(K(a+b+1))$. Since $\psi(s^{*})=-4 (K a-K-2 a-b) (K b+2
K+a-1)<0$, we have $s^{*}\in(s_2, s_3)$.

Summarizing the above  discussions, we arrive at the following results.
\begin{lemm}\label{L2.1} For the positive equilibrium $E^{*}(1, 1)$, the following statements hold.
\begin{itemize}
\item[(a)] Suppose $K b+2 K+a-1<0$, then $E^{*}(1, 1)$ is a hyperbolic saddle;
\item[(b)] Suppose $K b+2 K+a-1=0$, then $E^{*}(1, 1)$ is a degenerate equilibrium;
\item[(c)] Suppose $K b+2 K+a-1>0$, then $E^{*}(1, 1)$ is a node, if either $K a+2 K b+3 K-b-2\leq 0$ or $K a+2 K b+3 K-b-2>0$ and $s\in(0, s_2)\cup(s_3, +\infty)$.  And it is a focus or center, if $s\in(s_2, s_3)$. More precisely,
\begin{itemize}
\item[(c1)] if $K a-K-2a-b<0$,  then $E^*$ is a locally asymptotically stable node, if either $K a+2 K b+3 K-b-2\leq 0$ or $K a+2 K b+3 K-b-2>0$ and $s\in(0, s_2)\cup(s_3, +\infty)$.  And it is a locally asymptotically stable focus, if  $K a+2 K b+3 K-b-2>0$ and  $s\in(s_2, s_3)$.
\item[(c2)] if $K a-K-2a-b=0$, then $E^{*}$ is a locally asymptotically stable node (or focus), if $s\in(s_1, +\infty)$ (or $s\in(0, s_1)$).
\item[(c3)]  Let  $K a-K-2a-b>0$.
\begin{itemize}
\item[(i)] If $s^{*}<s<s_3$ (or $s>s_3$), then $E^{*}(1, 1)$ is a locally asymptotically stable focus (or node);
\item[(ii)]  If $s_2<s<s^{*}$ (or $s<s_2$), then $E^{*}(1, 1)$ is an unstable focus (or node);
\item[(iii)] If $s=s^{*}$, then $E^{*}(1, 1)$ is a weak focus or center.
\end{itemize}
\end{itemize}
\end{itemize}
\end{lemm}
The number of positive  equilibria of
system \eqref{norm} is determined by the number of roots of Eq. \eqref{dai2.13} on the interval $(0, K)$.
Obviously,  Eq. \eqref{dai2.13} can have one, two or
three positive roots on the interval $(0, K)$. Correspondingly, system \eqref{norm} can have one, two, or three
positive equilibria. For system \eqref{norm}, it can be inferred from \cite{Hsu1995} that if $\tilde{x}$ is a multiple positive root of Eq. \eqref{dai2.13}, then $\tilde{E}(\tilde{x},
\tilde{y})$ must be a degenerate positive equilibrium of system \eqref{norm}.
Applying the results of Lemma \ref{L2.1}, after tedious analysis, we can get the following results which are equivalent to Lemma 2.1 of \cite{Huang2014}.
\begin{lemm}\label{L2.2}
Let $\bar{\Delta}=(K b+K-1)^2-4aK$. Then the following statements hold.
\begin{itemize}
\item[(a)] Suppose  $b>1/K-1-2\sqrt{a/K}$, then system \eqref{norm} has a unique positive equilibrium $E^{*}(1, 1)$, which is an elementary and anti-saddle equilibrium;
\item[(b)] Suppose  $b=1/K-1-2\sqrt{a/K}$.
\begin{itemize}
\item[(b1)] If $a=K$, i.e., $b=1/K-3$, then system \eqref{norm} has a unique positive equilibrium $E^{*}(1, 1)$, which is a degenerate equilibrium;
\item[(b2)] If $a\neq K$, then system \eqref{norm} has two different positive equilibria: a degenerate equilibrium $E_2^{*}(\sqrt{K/a}, \sqrt{K/a})$ and an elementary anti-saddle
    equilibrium $E^{*}(1, 1)$;
\end{itemize}
\item[(c)] Suppose  $b<1/K-1-2\sqrt{a/K}$.
\begin{itemize}
\item[(c1)] If $b=(1-a)/K-2$, then system \eqref{norm} has two different positive equilibria: a degenerate equilibrium $E^{*}(1, 1)$ and an elementary anti-saddle equilibrium $E_2^{*}(K/a,
    K/a)$;
\item[(c2)] If $b\neq (1-a)/K-2$, then system \eqref{norm} has three distinct positive equilibria $E^{*}(1, 1)$, $E_2^{*}((1-Kb-K-\sqrt{\bar{\Delta}})/(2a), (1-Kb-K-\sqrt{\bar{\Delta}})/(2a))$ and
    $E_3^{*}((1-K b-K+\sqrt{\bar{\Delta}})/(2 a), (1-K b-K+\sqrt{\bar{\Delta}})/(2 a))$, which are all elementary equilibria. Without loss of generality, assume that $1$ is the minimum root of  Eq.
    \eqref{dai2.13}. Then both $E^{*}$ and $E_3^{*}$ are anti-saddle, and $E_2^{*}$ is a saddle.
\end{itemize}
\end{itemize}
\end{lemm}
The cases (b) and (c1) of Lemma \ref{L2.2} have been investigated by Huang et al \cite{Huang2014}. In the rest of this paper  we only need to focus on the cases (a) and (c2) of Lemma \ref{L2.2}.
Noting that Huang et al \cite{Huang2014} considered all possible cases when system has at least one degenerate positive equilibrium, we only need to consider the case that $E^{*}(1, 1)$
is non-degenerate.
In order words, the type and stability of $E^{*}(1, 1)$ have been determined for all the cases of Lemma \ref{L2.1} except the subcase (iii). Hence, we only need to  determine if $E^{*}(1, 1)$ is a center or a weak focus under the corresponding conditions.

\begin{rema}\label{rem2.2}
The qualitative behavior near the origin shows that the origin is always a repeller in the interior of the first quadrant, which corrects the results in \cite{Eduardo2014} that the origin can be a global attractor.
It's pointed out that we analyze the isolated critical point $O$ of higher order with a different approach. By Lemma \ref{L2.1}, we know that the origin $O$ cannot be an $\omega$ limit set of the unstable manifold of the hyperbolic saddle $E_0$, i.e., there is no heteroclinic orbit. This also corrects the wrong results in \cite{Eduardo2014} that there exists a heteroclinic curve joining the hyperbolic saddle $E_0$ and the origin $O$.
\end{rema}

\subsection{Computation of Lyapunov constants}

In this subsection, we are going to study Hopf bifurcation of system \eqref{norm}.
Without loss of generality,  assume that $E^{*}(1, 1)$ is center or focus type, which implies that the condition (iii) in Lemma \ref{L2.1} (c3) holds, i.e.,
\begin{equation}\label{cond2.3}
\begin{aligned}
K b+2 K+a-1>0,\ K a-K-2 a-b>0,\  s=s^{*}.
\end{aligned}
\end{equation}

To determine if $E^{*}(1, 1)$ is a center or a weak focus, one needs to calculate the focal values of system \eqref{norm} at this equilibrium.
In this paper, we will use the Lyapunov
constants instead of the focal values to solve these problems.
The equivalence between the Lyapunov constants and the focal values can be seen in \cite{Liu2001} and Chapter $1$ in \cite{Liu2008}.

Denote the $k$-th Lyapunov constants of system \eqref{norm} at  $E^{*}(1, 1)$ by $V_{2k+1}$, $k=0, 1, \cdots$. Then $V_1=\textrm{Tr}\left(J(E^{*}(1, 1))\right)$. Under the hypothesis $V_1=0$,
i.e., $s=s^{*}$, we will compute higher Lyapunov constants of system \eqref{norm} at $E^{*}(1, 1)$. To this end, we consider a
conjugation $\psi: \mathbb{R}^2\rightarrow \mathbb{R}^2$, such that
\begin{equation*}
\psi (u, v)=(u, v)+(1, 1)=(x, y).
\end{equation*}
For brevity, we denote $d\triangleq Kb+2K+a-1>0,\ e\triangleq K a-K-2a-b>0$.
Applying the algorithm in \cite{Sang1} and using the software Maple for symbolic calculas, we have
\begin{equation}\label{dai2.19}
\begin{aligned}
V_3=-\dfrac{(K-1)^2 \varphi_1}{4d},
V_5=\dfrac{(K-1)^3 \varphi_2}{48d^3e},
V_7=-\dfrac{(K-1)^4 \varphi_3}{9216d^5e^2},
V_9=\dfrac{(K-1)^5 \varphi_4}{1105920d^7e^3},
\end{aligned}
\end{equation}
where the quantity $V_k$ is reduced w.r.t. the Gr\"{o}bner basis of $\{V_j: j<k\}$ and the expressions of $\varphi_i$'s, $i=1, 2, 3, 4$, we omit here for their cumbersome expressions are available upon request.
This requires a large-scale computation. For example, for $i=1, 2, 3, 4$, the number of  terms of $\varphi_i$ is, respectively, $29$, $340$, $1216$ and $2947$. From \eqref{cond2.3}, we shall consider the Lyapunov constants in \eqref{dai2.19} with $\mu=(K, a, b)\in \Lambda$, where
\begin{equation}\label{basecond}
\Lambda=\left\{\mu\in\mathbb{R}^3\left|K>1, a>0, \textup{max}\bigg\{\frac{1-a}{K}-2, -2\sqrt{a}\bigg\}<b<K a-K- 2 a\right.\right\}.
\end{equation}

\begin{rema}\label{R3.2}
From Lemmas \ref{L2.1} and \ref{L2.2}, if the conditions in the case (a) of Lemma \ref{L2.2} and the case (ii) of Lemma \ref{L2.1} (c2) hold, i.e.,
\begin{equation*}
\begin{aligned}
\textup{max}\bigg\{\frac{1}{K}-1-2\sqrt{\frac{a}{K}},\ -2\sqrt{a}\bigg\}<b<Ka-K-2a, \ s<s^{*},
\end{aligned}
\end{equation*}
then $E^{*}(1, 1)$ is the unique positive equilibrium of system \eqref{norm} and is unstable. Hence system \eqref{norm} has at least one stable limit cycle in $\Omega$ from Theorem \ref{T3.2}.
\end{rema}

\section{Proofs}

Notice that system \eqref{norm2} is equivalent to system \eqref{model} in $\mathcal{A}$. Next we only need to consider the polynomial system \eqref{norm}, instead of system \eqref{model}.
\subsection{Proof of Theorem \ref{T3.2}}
If $\Delta>0$, then the equation \eqref{dai1.2} has a unique positive root in the interval $(0, K)$.
Correspondingly, system \eqref{model} has a unique positive equilibrium.
Similar to the analysis of Lemma 2.1 in \cite{Huang2014}, the first assertion is derived immediately.

Next we will prove the second assertion. From Lemma \ref{T3.1}, every solution $(x(t),  y(t))$ of system \eqref{norm} with positive initial values will eventually be away from the origin (see Fig. 2.1). Hence, the boundary
of the region $\Omega$ can be used as the outer boundary of a Poincar\'{e}-Bendixson annular region. By applying the Poincar\'{e}-Bendixson Theorem, the conclusion of this theorem is clearly established.

\subsection{Proof of Theorem \ref{T4.1}}

From the analysis in Section $2$, the conditions of this theorem imply that $E^{*}(1, 1)$ is the unique positive equilibrium of system \eqref{norm} and is center or focus type.
Noting that $b>1/K-1-2\sqrt{a/K}$ implies that $b>(1-a)/K-2$. From Lemmas \ref{L2.1} and \ref{L2.2}, if the conditions in Theorem \ref{T4.1} hold, then $s=s^{*}$ and $\mu\in\Sigma_1$, where
\begin{equation*}
\Sigma_1=\left\{\mu\in\mathbb{R}^3\left|K>1, a>0, \textup{max}\left\{\frac{1}{K}-1-2\sqrt{\frac{a}{K}},\
-2\sqrt{a}\right\}<b<K a-K-2 a\right.\right\}.
\end{equation*}

To prove this theorem, we first prove that $\varphi_1, \varphi_2, \varphi_3$ have a unique common real root in $\Sigma_1$.
We consider these polynomials in the ring $\textbf{R}[K, a, b]$.

Firstly, computing the mutual resultants of $\varphi_{1}, \varphi_{2}, \varphi_{3}$ with respect to $b$, we get
\begin{equation}\label{dai4.2}
\left\{
\begin{aligned}
&\varphi_{12}:=\textrm{Res}(\varphi_{1}, \varphi_{2}, b)
=32K^{4}a(K-1)^{13}(a-1)^{18}(a-K)^2g_{1b}g_{2}^2g_{31}, \\
&\varphi_{13}:=\textrm{Res}(\varphi_{1}, \varphi_{3}, b)=512K^{5}a(K-1)^{21}(a-1)^{28}(a-K)^3g_{1b}g_{2}^2g_{32}, \\
&\varphi_{23}:=\textrm{Res}(\varphi_{2}, \varphi_{3}, b)=-480K^{10}a(K-1)^{63}(a-1)^{64}(a-K)^{10}g_{1b}g_{2}^2g_{33},
\end{aligned}
\right.
\end{equation}
where $g_{1b}=K^3 a-4 K^2 a^2-K^3+10 K a^2+6 K^2-11 K a-4 a^2-3 K+6 a$, $g_{2}=2{K}^{2}+K-1-(3 K-1)a$, and
$g_{31}$, $g_{32}$, $g_{33}$ are three complicated polynomials in $K$, $a$ of degree $19$, $39$,  $103$, respectively. And more, $\varphi_{12}, \varphi_{13}, \varphi_{23}$ are well factored
over the rational field.
Using the resultant elimination theory (Theorem \ref{L4.3} in Appendix A),  we have $\textbf{V}(\varphi_1, \varphi_2, \varphi_3)=\textbf{V}(\varphi_1, \varphi_2, \varphi_3, \varphi_{12},
\varphi_{13}, \varphi_{23})$.
To find the common real roots of $\{\varphi_1, \varphi_2, \varphi_3, \varphi_{12}, \varphi_{13}, \varphi_{23}\}$, from \eqref{dai4.2}, we have the following five cases to discuss.

\textbf{Case} (a): $a=1$. Substituting it to $\varphi_1$, we get
$\varphi_{1}=-(b+4) (b+2)^2(K^2-4 K+1-K b)$.
Hence $\varphi_1=0$ implies that $b=-2$, $-4$ or $(K^2-4 K+1)/K$.
The cases $b=-2$ and  $b=-4$ can be excluded since $b>-2\sqrt{a}$.
The case $b=(K^2-4 K+1)/K$ follows that $K a-K-2 a-b=-(K-1)^2/K<0$, which contradicts to
$K a-K-2 a-b>0$. Therefore, for this case, we have $\textbf{V}(\varphi_1)=\emptyset$.

\textbf{Case} (b): $a=K$. With this substitution, we have
\begin{equation}\label{case(b)}
\varphi_{1}=(K b+3K-1)\phi_1, \,\,
\varphi_{2}=(K b+3K-1)^2(K^2-3 K-b)\phi_2,
\end{equation}
where $\phi_{1}$ and $\phi_{2}$ are polynomials in $K$, $b$ of degree $4$ and $9$, respectively.
If $b=1/K-3$, then $K b+2 K+a-1=0$, which contradicts to \eqref{cond2.3}.
If $b=K^2-3K$, then $K a-K-2 a-b=0$, which contradicts to  \eqref{cond2.3}.
Thus, it follows from \eqref{case(b)} that $\textbf{V}(\varphi_1, \varphi_2)=\textbf{V}(\phi_1, \phi_2)$. By Sturm's Theorem, we conclude that the resultant
$\textrm{Res}(\phi_{1}, \phi_{2}, b)$ is a nonzero constant, implying that $\textbf{V}(\phi_1, \phi_2)=\emptyset$. Therefore, for this case, we have $\textbf{V}(\varphi_1,
\varphi_2)=\emptyset$.

\textbf{Case} (c): $g_2=0$, i.e., $a=(2K^2+K-1)/(3K-1)$. Then we have
\begin{equation}\label{case(c)}
\varphi_{1}=
(3 K b+8 K-b-4)\phi_{4}/(3K-1)^3, \,\,
\varphi_{2}=(3 K b+8 K-b-4)\phi_{5}/(3K-1)^8,
\end{equation}
where $\phi_{4}$, $\phi_{5}$ are  polynomials in $K$, $b$ of degree $7$, $21$, respectively.
If $3 K b+8 K-b-4=0$, then $b=-4(2K-1)/(3K-1)$. This follows $K b+2 K+a-1=0$, which contradicts to \eqref{cond2.3}. Thus, it follows from \eqref{case(c)} that $\textbf{V}(\varphi_1,
\varphi_2)=\textbf{V}(\phi_4, \phi_5)$. Similar to the case (b), we conclude that $\textbf{V}(\phi_4, \phi_5)=\emptyset$. Hence, for this case, we have  $\textbf{V}(\varphi_1,
\varphi_2)=\emptyset$.

\textbf{Case} (d): $g_{31}=g_{32}=g_{33}=0$. Calculating the  mutual resultants of $g_{31}, g_{32}, g_{33}$  with respect to $a$, we get
\begin{equation}\label{dai4.3}
\begin{aligned}
&\textrm{Res}(g_{31}, g_{32}, a)=
c_1 K^{12}(2K-1)^2(K-1)^{147}(19 K^2-6 K+2) h_1 h_2 h_3 h_4, \\
&\textrm{Res}(g_{31}, g_{33}, a)=
c_2 K^{12}(2K-1)^2(K-1)^{372}(19 K^2-6 K+2) h_1 h_2 h_3 h_5,
\end{aligned}
\end{equation}
where $c_1$, $c_2$ are two constants, and $h_1, h_2, h_3$, $h_4$ and $h_5$ are polynomials in $K$ of degree $3$, $8$, $42$, $274$ and $830$, respectively.  Obviously, we have $19 K^2-6 K+2>0$
since $K>1$. In addition, by Sturm's Theorem, we conclude that $h_1>0$ and $h_2>0$ for $(1, +\infty)$. Furthermore, $h_4$ and $h_5$ have no common zeros because the resultant $\textrm{Res}(h_4,
h_5, K)$ is a nonzero constant.
Thus, from \eqref{dai4.3} and Theorem \ref{L4.3}  in Appendix A, in this case we have $\textbf{V}(\varphi_1, \varphi_2, \varphi_3)=\textbf{V}(\varphi_1, \varphi_2, \varphi_3, g_{31}, g_{32},
g_{33}, h_3)=\textbf{V}(h_3, g_{31}, \varphi_1)\cap\textbf{V}(g_{32}, g_{33}, \varphi_2, \varphi_3)$.
Next, we first find out the common real roots of $\{h_3, g_{31}, \varphi_1\}$ in $\Sigma_1$, and then verify whether they are really the common real roots of $\{\varphi_2, \varphi_3\}$. Even
though $\{h_3, g_{31}, \varphi_1\}$ is an irreducible ascending set, we can not use the $\textsl{RealRootIsolate}$ command in Maple to obtain the isolated real roots in $\Sigma_1$ directly for
their extremely complicated expressions. Here, we will use the real root isolation algorithm of multivariate polynomial systems, proposed by Lu et al \cite{Lu2007},   to solve this problem.
According to the sign of $b$, we have the following three subcases to discuss.

\textbf{Subcase} (d1): $b>0$. Without specification, taking the accuracy of command \emph{realroot} in Maple always being  $1/10^{20}$.
By command $\textsl{realroot}$, we get that $h_3(K)$ has two real roots in $(1, +\infty)$ as follows.
\begin{align*}
&K_{1}\in \left[{\dfrac {354982973071688257601}{295147905179352825856}}, {\dfrac {
1419931892286753030405}{1180591620717411303424}}\right]\triangleq \left[\underline{K}_{1}, \overline{K}_{1}\right],\\
&K_{2}\in \left[{\dfrac {244229161108537248689527983}{1180591620717411303424}},{\dfrac
{15264322569283578043095499}{73786976294838206464}}\right]\triangleq \left[\underline{K}_{2}, \overline{K}_{2}\right].
\end{align*}
For the real root interval $[\underline{K}_{1}, \overline{K}_{1}]$ of $h_3(K)$,
we get the following two positive real root isolation intervals of $\{h_3(K), g_{31}(K, a)\}$.
\begin{align*}
&a_{1}\in\left[{\dfrac {75253211057447250009}{295147905179352825856}},{\dfrac {
150506422114898607437}{590295810358705651712}}\right]\triangleq \left[\underline{a}_{1}, \overline{a}_{1}\right],\\
&a_{2}\in\left[{\dfrac {494428213075491687541}{295147905179352825856}},{\dfrac {
1977712852306468880883}{1180591620717411303424}}\right]\triangleq \left[\underline{a}_{2}, \overline{a}_{2}\right].
\end{align*}
From Definition 2.4 of \cite{Lu2007}, we denote the maximal and minimal polynomials
of $\varphi_{1}$ on $G_{1}=[\underline{K}_{1}, \overline{K}_{1}]\times[\underline{a}_{1}, \overline{a}_{1}]$
by $\overline{\varphi}_1(b)$ and $\underline{\varphi}_1(b)$, respectively.
By command \emph{realroot}, we know that both $\overline{\varphi}_{1}(b)$ and $\underline{\varphi}_{1}(b)$ have no positive real root. Hence $\varphi_{1}$ has no positive real  root for the
real root isolation interval $G_{1}$. Similarly, $\varphi_{1}$ has no positive real root for the real root isolation interval $G_{2}=[\underline{K}_{1},
\overline{K}_{1}]\times[\underline{a}_{2}, \overline{a}_{2}]$.

Similarly, for the real root interval $[\underline{K}_{2}, \overline{K}_{2}]$ of $h_3(K)$, we can get three positive real root isolation intervals of $\{h_3(K), g_{31}(K, a)\}$ as below.
\begin{align*}
&a_{3}\in\left[{\dfrac {302357397559963930735}{590295810358705651712}},{\dfrac {604714795119927861471}{1180591620717411303424}}
\right]\triangleq \left[\underline{a}_{3}, \overline{a}_{3}\right],\\
&a_{4}\in\left[{\dfrac {111474188614417001592773}{590295810358705651712}}
,{\dfrac {222948377228834003185547}{1180591620717411303424}}
\right]\triangleq \left[\underline{a}_{4}, \overline{a}_{4}\right],\\
&a_{5}\in\left[{\dfrac {4494291717948679361239971}{1180591620717411303424}}
,{\dfrac {1123572929487169840309993}{295147905179352825856}}\right]\triangleq \left[\underline{a}_{5}, \overline{a}_{5}\right].
\end{align*}
For the real root isolation interval $G_{3}=[\underline{K}_{2}, \overline{K}_{2}]\times[\underline{a}_{3}, \overline{a}_{3}]$ of $\{h_3, g_{31}\}$, similar to the approach for $G_{1}$, we get
that $\varphi_{1}$ has a positive real root at $b_{1}$, where
\begin{equation*}
b_{1}\in\left[{\dfrac {125098419344119784764828183}{1180591620717411303424}}
,{\dfrac {62549209672059892382517527}{590295810358705651712}}
\right]\triangleq \left[\underline{b}_{1}, \overline{b}_{1}\right].
\end{equation*}
By Theorem 2.3 of \cite{Lu2007}, we have $Ka-K-2a-b\leq \overline{K}_{2}\overline{a}_{3}-\underline{K}_{2}-2\underline{a}_{3}-\underline{b}_{1}
<0$, which implies that $(K_{2}, a_3, b_1)\not\in\Sigma_1$.

For the real root isolation interval $G_{4}=[\underline{K}_{2}, \overline{K}_{2}]\times[\underline{a}_{4}, \overline{a}_{4}]$ of $\{h_3(K), g_{31}(K, a)\}$, $\varphi_{1}$ has three positive
real roots at $b_{i}$ with $b_{i}\in[\underline{b}_{i}, \overline{b}_{i}]$, $i=2, 3, 4$, where
\begin{align*}
&\left[\underline{b}_{2}, \overline{b}_{2}\right]=\left[{\dfrac {192301334985162977663}{2361183241434822606848}}
,{\dfrac {769205339940651910653}{9444732965739290427392}}
\right],\\
&\left[\underline{b}_{3}, \overline{b}_{3}\right]=\left[{\dfrac {12252658027754888276693}{1180591620717411303424}}
,{\dfrac {6126329013877444138347}{590295810358705651712}}
\right],\\
&\left[\underline{b}_{4}, \overline{b}_{4}\right]=\left[{\dfrac {46120025694613890773675092515}{1180591620717411303424}}
,{\dfrac {46120025694613890773675299959}{1180591620717411303424}}
\right].
\end{align*}
By Theorem 2.3 of \cite{Lu2007}, we have
$\varphi_2(K_{2}, a_4, b_3)\geq {\varphi}^+_{2}(\underline{K}_{2}, \underline{a}_{4}, \underline{b}_{3})+{\varphi}^-_{2}(\overline{K}_{2}, \overline{a}_{4}, \overline{b}_{3})\approx 1.83\times
10^{46}>0$, implying that $(K_{2}, a_4, b_3)$ is not a common real root of $\{\varphi_1, \varphi_2, \varphi_3\}$. In addition,
$(K_{2}, a_4, b_4)\not\in \Sigma_1$ since $K a-K-2 a-b\leq \overline{K}_{2}\overline{a}_{4}-\underline{K}_{2}-2\underline{a}_{4}-\underline{b}_{4}<0$.

It's not difficult to verify that the real root $(K_2, a_4, b_2)\in\Sigma_1$. Similarly, after tedious calculation, we can verify that $(K_2, a_4, b_2)$ is a common real root of $\{h_3, g_{31},
\varphi_2\}$ and also a common real root of $\{h_3, g_{31}, \varphi_3\}$. Thus, we conclude that $(K_2, a_4, b_2)$ is a common real root of $\{\varphi_1, \varphi_2, \varphi_3\}$ in $\Sigma_1$.

For the real root isolation interval $G_{5}=[\underline{K}_{2}, \overline{K}_{2}]\times[\underline{a}_{5}, \overline{a}_{5}]]$ of $\{h_3(K), g_{31}(K, a)\}$, $\varphi_{1}$ has two positive real
roots at $b_{i}$ with $b_{i}\in[\underline{b}_{i}, \overline{b}_{i}]$, $i=5, 6$, where
\begin{align*}
&\left[\underline{b}_{5}, \overline{b}_{5}\right]=\left[{\dfrac {71394370384755146991531}{1180591620717411303424}}
,{\dfrac {17848592596188786747883}{295147905179352825856}}
\right],\\
&\left[\underline{b}_{6}, \overline{b}_{6}\right]=\left[{\dfrac {464853909302125803012485143759}{590295810358705651712}}
,{\dfrac {929707818604251606024970505815}{1180591620717411303424}}
\right].
\end{align*}
The real root $(K_{2}, a_{5}, b_{5})$ is not a common real root of $\{\varphi_1, \varphi_2, \varphi_3\}$ since
$\varphi_3(K_{2}, a_{5}, b_{5})\geq {\varphi}^+_{3}(\underline{K}_{2}, \underline{a}_{5}, \underline{b}_{5})+{\varphi}^-_{3}(\overline{K}_{2}, \overline{a}_{5}, \overline{b}_{5})\approx
2.085\times 10^{91}>0.$
On the other hand, $(K_{2}, a_{5}, b_{6})\not\in \Sigma_1$ because
$Ka-K-2a-b\leq \overline{K}_{2}\overline{a}_{5}-\underline{K}_{2}-2\underline{a}_{5}-\underline{b}_{6}\approx -1.92\times 10^{5}
<0$.

Therefore, for the subcase (d1), we know that $p^{*}=(K_2, a_4, b_2)$ is the unique common real root of $\{\varphi_1, \varphi_2, \varphi_3\}$ in $\Sigma_1$.

\textbf{Subcase} (d2): $b<0$. For this subcase, we will find the common real roots of $\{h_3, g_{31}, \varphi_1\}$ in the cone $(1, +\infty)\times(0, +\infty)\times(-\infty, 0)$  by taking the
transformation $b\rightarrow -b$.
We denote
$\hat{\varphi}_1(K, a, b)=\varphi_1(K, a, -b)$.
Since the variable $b$ exactly appears in $\varphi_1(K, a, b)$ for the triangular set $\{h_3, g_{31}, \hat{\varphi}_1\}$, we only need to find the common real roots of  $\{h_3, g_{31},
\hat{\varphi}_1\}$ in the cone $(1, +\infty)\times(0, +\infty)\times(0, +\infty)$.
Repeating the same arguments as the subcase (d1), we conclude that $\{h_3, g_{31}, \hat{\varphi}_1\}$ have fourteen classes of real roots in the cone $(1, +\infty)\times(0, +\infty)\times(0,
+\infty)$: $(K_1, a_1, b_i)$, $(K_1, a_2, b_j)$, $(K_2, a_3, b_k)$, $(K_2, a_4, b_{18})$, $(K_2, a_5, b_{19})$, $(K_2, a_5, b_{20})$, $i=7, \cdots, 10, j=11, \cdots, 14, k=15, \cdots, 17$,
where $K_1$, $K_2$, $a_3$,  $a_4$ and $a_5$ are the same as in the subcase (d1) and
\begin{align*}
&b_{7}\in\left[{\dfrac {342692738901203887719}{590295810358705651712}}
,{\dfrac {685385477803542424951}{1180591620717411303424}}
\right]\triangleq \left[\underline{b}_{7}, \overline{b}_{7}\right],\\
&b_{8}\in\left[{\dfrac {1338789934782848805353}{1180591620717411303424}}
,{\dfrac {334697483698195198185}{295147905179352825856}}
\right]\triangleq \left[\underline{b}_{8}, \overline{b}_{8}\right],\\
&b_{9}\in\left[{\dfrac {194243729812573167855}{147573952589676412928}}
,{\dfrac {776974919255073511755}{590295810358705651712}}
\right]\triangleq \left[\underline{b}_{9}, \overline{b}_{9}\right],\\
&b_{10}\in\left[{\dfrac {3967870668570439988007}{1180591620717411303424}}
,{\dfrac {3967870668571263500067}{1180591620717411303424}}
\right]\triangleq \left[\underline{b}_{10}, \overline{b}_{10}\right],\\
&b_{11}\in\left[{\dfrac {2889636985159625502057}{1180591620717411303424}}
,{\dfrac {2889637013691542396369}{1180591620717411303424}}
\right]\triangleq \left[\underline{b}_{11}, \overline{b}_{11}\right],\\
&b_{12}\in\left[{\dfrac {386107933062469966493}{147573952589676412928}}
,{\dfrac {1544431759404950266357}{590295810358705651712}}
\right]\triangleq \left[\underline{b}_{12}, \overline{b}_{12}\right],\\
&b_{13}\in\left[{\dfrac {3400585485349590843693}{1180591620717411303424}},{\dfrac {1700292756101773392567}{590295810358705651712}}
\right]\triangleq \left[\underline{b}_{13}, \overline{b}_{13}\right],\\
&b_{14}\in\left[{\dfrac {3105247046679173264103}{590295810358705651712}},{\dfrac {6210494094466506191137}{1180591620717411303424}}\right]\triangleq \left[\underline{b}_{14},
\overline{b}_{14}\right],\\
&b_{15}\in\left[{\dfrac {821227566792260176395}{590295810358705651712}},{\dfrac {1642455133584520352861}{1180591620717411303424}}\right]\triangleq \left[\underline{b}_{15},
\overline{b}_{15}\right],\\
&b_{16}\in\left[{\dfrac {497030357386605526709}{295147905179352825856}},{\dfrac {1988121429546422106927}{1180591620717411303424}}\right]\triangleq \left[\underline{b}_{16},
\overline{b}_{16}\right],\\
&b_{17}\in\left[{\dfrac {4031280671528406533751}{590295810358705651712}}
,{\dfrac {8062561343056813067555}{1180591620717411303424}}\right]\triangleq \left[\underline{b}_{17}, \overline{b}_{17}\right],\\
&b_{18}\in\left[{\dfrac {19438401444953006468837}{1180591620717411303424}}
,{\dfrac {9719200722476503234419}{590295810358705651712}}\right]\triangleq \left[\underline{b}_{18}, \overline{b}_{18}\right],\\
&b_{19}\in\left[{\dfrac {973745532080307675583}{4722366482869645213696}},{\dfrac {15214773938754807431}{73786976294838206464}}\right]\triangleq \left[\underline{b}_{19},
\overline{b}_{19}\right],\\
&b_{20}\in\left[\dfrac {78104754415268652314327}{1180591620717411303424}
,\dfrac {9763094301908581539291}{147573952589676412928}\right]\triangleq \left[\underline{b}_{20}, \overline{b}_{20}\right].
\end{align*}
It's pointed out that $(K_i, a_j, b_k)$ is a common real root of  $\{h_3, g_{31}, \hat{\varphi}_1\}$ if and only if $(K_i, a_j, -b_k)$ is a common real root of  $\{h_3, g_{31}, \varphi_1\}$.
It's easy to verify that the following eight classes of
real roots of $\{h_3, g_{31}, \varphi_1\}$ contradict with $b>-2\sqrt{a}$.
\begin{equation*}
(K_1, a_1, -b_{i}), (K_1, a_2, -b_{j}), (K_2, a_3, -b_{k}),\ i=8, 9, 10,\ j=12, 13, 14,\ k=16, 17.
\end{equation*}
By Theorem 2.3 of \cite{Lu2007}, we can verify that the real roots $(K_1, a_1, -b_{7})$, $(K_1, a_2, -b_{11})$ and $(K_2, a_3, -b_{15})$  contradict to $K a-K-2 a-b>0$. In addition, the real
roots $(K_2, a_4, -b_{18})$ and $(K_2, a_5, -b_{20})$ contradict to $K b+2 K+a-1>0$.
Furthermore, $(K_{2}, a_{5}, b_{19})$ is not a common real root of $\{\varphi_1, \varphi_2, \varphi_3\}$ since
$\hat{\varphi}_2(K_{2}, a_{5}, b_{19})\leq \hat{\varphi}^+_{2}(\overline{K}_{2}, \overline{a}_{5}, \overline{b}_{19})+\hat{\varphi}^-_{2}(\underline{K}_{2}, \underline{a}_{5},
\underline{b}_{19})\approx -2.07\times 10^{50}<0,$
where $\hat{\varphi}_2(K, a, b)=\varphi_2(K, a, -b)$.

Therefore, for the subcase (d2), we have  $\textbf{V}(\varphi_1, \varphi_2, \varphi_3)=\emptyset$.

\textbf{Subcase} (d3): $b=0$.  With this substitution,  we have
$\varphi_{1}=-4a \psi_1$ and $\varphi_{2}=-4a \psi_2$,
where $\psi_1$ and $\psi_2$ are polynomials in $K$, $a$ of degree $3$ and $10$, respectively.
Thus, $\textbf{V}(\varphi_1, \varphi_2)=\textbf{V}(\psi_1, \psi_2)$.
The resultant of $\psi_1$  and $\psi_2$ with respect to $a$ is
\begin{equation}\label{H2}
\textrm{Res}(\psi_1, \psi_2, a)=4096(3 K-1)^{2}(K-1)^{13}H_{1}H_{2},
\end{equation}
where $H_{1}$ is a polynomial in $K$ of degree $5$ and $H_{2}=K^2-6 K+3$. By Sturm's Theorem, we know that $\textrm{Res}(\psi_1, \psi_2, a)=0$ implies $H_{2}=0$. This follows $K=3+\sqrt{6}$.
With this substitution, we get that $\psi_1=(8+3\sqrt{6})(a+3)(3+2\sqrt{6}-5a)/5=0$ implies $a=(3+2\sqrt{6})/5$.
Hence, $K a-K-2 a-b=0$, which contradicts to \eqref{cond2.3}. Therefore, for this case, we have $\textbf{V}(\varphi_1, \varphi_2)=\emptyset$.

To summarize, $\{\varphi_1, \varphi_2, \varphi_3\}$ have a unique common real root
$p^{*}=(K_2, a_4, b_2)\in\Sigma_1$ for all the cases above.

\textbf{Case} (e): $g_{1b}=0$.
From above, if $\{\varphi_1, \varphi_2, \varphi_3\}$ have other common real roots in $\Sigma_1$ except the real root $p^{*}$, then there must be $g_{1b}=0$.

Computing the mutual resultants of $\varphi_{1}, \varphi_{2}, \varphi_{3}$ with respect to $a$, we get
\begin{equation}\label{dai4.4}
\left\{
\begin{aligned}
&\textrm{Res}(\varphi_{1}, \varphi_{2}, a)=-8b(b+2)^{18}(K-1)^{13}(K b+3 K-1)^2 g_{1a}\tilde{g}^2_{2}\tilde{g}_{31}, \\
&\textrm{Res}(\varphi_{1}, \varphi_{3}, a)=-64b(b+2)^{28}(K-1)^{21}(K b+3 K-1)^2 g_{1a}\tilde{g}^2_{2}\tilde{g}_{32}, \\
&\textrm{Res}(\varphi_{2}, \varphi_{3}, a)=-120b(b+2)^{64}(K-1)^{63}
(K b+3K-1)^8 g_{1a}\tilde{g}^2_{2}\tilde{g}_4^2\tilde{g}_{33},
\end{aligned}
\right.
\end{equation}
where $g_{1a}={K}^{3}b+4\,{K}^{3}-8{K}^{2}b-4K{b}^{2}-24{K}^{2}-7Kb+2{b}^{2}+12K+6b$, $\tilde{g}_{2}=3 K b+8 K-b-4$,
$\tilde{g}_{4}=K^2-3 K-b$, and $\tilde{g}_{31}$, $\tilde{g}_{32}$, $\tilde{g}_{33}$ are three polynomials in $K$, $b$ of degree $18$, $37$, $102$, respectively.
Similar to the proof of eliminating  the variable $b$ above, we can prove that there must be $g_{1a}=0$, provided that $\{\varphi_1, \varphi_2, \varphi_3\}$ have other common real roots in
$\Sigma_1$  except $p^{*}$.

In what follows, we will prove, by contradiction, that $\{\varphi_1, \varphi_2, \varphi_3\}$ have no common real roots in $\Sigma_1$ such that  $g_{1b}=g_{1a}=0$.
Otherwise, assume that there exist parameter values in $\Sigma_1$ such that $\varphi_1= \varphi_2=\varphi_3=g_{1b}=g_{1a}=0$ hold.

Firstly, applying  pseudo division and by command \emph{prem} in Maple, we get
\begin{equation*}
(2-4K)^3\varphi_1=q_1g_{1a}+r_1,
\end{equation*}
where $q_1$ and $r_1$ are polynomials in $K$, $a$, $b$ of degree $7$ and  $11$, respectively. It follows from $\varphi_1=g_{1a}=0$ that $r_1=0$, which yields $\textbf{V}(\varphi_1,
g_{1a})=\textbf{V}(g_{1a}, r_1)$. Noting that $r_1=H_3b-H_4=0$,
where both $H_3$ and $H_4$ are polynomials in $K$, $a$ of degree $10$.
If $H_3\neq 0$, then $b=H_4/H_3$ and hence $Ka-K-2a-b=-g_{1b}H_5/H_3=0$, where  $H_5$ is a polynomial in $K$, $a$ of degree $8$. This contradicts to \eqref{cond2.3}. So there must be $H_3=0$,
then $r_1=H_4=0$, which yields $\textbf{V}(g_{1a}, r_1)=\textbf{V}(g_{1a}, H_3, H_4)$.
Applying  pseudo division and by command \emph{prem} in Maple, we obtain
\begin{equation}\label{dai4.5}
(K-2)^2H_4=q_2g_{1b}+4(K-1)^2H_2H_6,
\end{equation}
where $H_6$ and $q_2$ are polynomials in $K$, $a$ of degree $8$ and $6$, respectively, and $H_2$ is given in Eq. \eqref{H2}.
If $K=2$, then $g_{1b}=-2(4a-5)$ and $g_{1a}=-2(3b+10)(b+2)$.
It follows from $g_{1b}=g_{1a}=0$ that $a=5/4$ and $b=-10/3$ or $-2$.
These two cases can be excluded since $b>-2\sqrt{a}$ and $Ka-K-2a-b>0$. Thus, we only need to consider the case $K\neq 2$.  From \eqref{dai4.5} and $H_4=g_{1b}=0$, we have $H_2H_6=0$. This
follows $H_2=0$ or $H_6=0$.
The case $H_2=0$ can be excluded. In fact, if $H_2=0$, i.e., $K=3+\sqrt{6}$, then we have
\begin{equation*}
g_{1b}=2(7\sqrt{6}+17) a (2\sqrt{6}-5a+3)/5, \,\, g_{1a}=-2(2\sqrt{6}+5) b(\sqrt{6}+b+3).
\end{equation*}
From the subcase (3) above, the case $b=0$ can be excluded.
Then $g_{1b}=g_{1a}=0$ implies $a=(3+2\sqrt{6})/5$ and $b=-3-\sqrt{6}$. This follows $K b+2 K+a-1=-(47+18\sqrt{6})/5<0$, which contradicts to \eqref{cond2.3}. Therefore, there must be $H_6=0$
and hence $\textbf{V}(g_{1a}, g_{1b}, H_4)=\textbf{V}(g_{1a}, g_{1b}, H_6)$. Notice that $H_6=-H_7a+H_8$, where $H_7$, $H_8$ are polynomials in $K$ of degree $7$, $8$, respectively.
The case $H_7=0$, implies $H_8=0$, can be excluded since the resultant $\textrm{Res}(H_7, H_8, K)$ is a nonzero constant. Thus, we have $H_7\neq0$ and $a=H_8/H_7$. Substituting it to
$g_{1b}=0$, we get
\begin{equation}\label{dai4.6}
g_{1b}=
2K(2K-1)(K-1)^2(K-2)^2H_9/H_7^2=0,
\end{equation}
where $H_9$ is a polynomial in $K$ of degree $9$.
The analysis above tells us that $K\neq 2$, then from \eqref{dai4.6}, there must be $H_9=0$ and hence $\textbf{V}(g_{1b}, H_6)=\textbf{V}(H_6, H_9)$.

To sum up,  we have $\textbf{V}(\varphi_1, g_{1a}, g_{1b})=\textbf{V}(g_{1a}, g_{1b}, r_1)=\textbf{V}(g_{1a}, g_{1b}, H_3, H_4)=\textbf{V}(g_{1a}, g_{1b}, H_3, H_6)=\textbf{V}(g_{1a}, H_3, H_6,
H_9)=\textbf{V}(g_{1a}, H_6, H_9)\cap\textbf{V}(H_3)$.

Now we will show that $\textbf{V}(g_{1a}, H_6, H_9)=\emptyset$. For convenience, denote  by  $S_a$ a semi-algebraic system whose polynomial equations, non-negative polynomial inequalities,
positive polynomial inequalities and polynomial inequations are given by $F_a:=[H_9, H_6, g_{1a}]$, $N:=[\ ]$, $P:=[K-1, a, K a-K-2a-b, K b+2 K+a-1]$, and $H:=[\ ]$, respectively, where $[\ ]$
represents the null set.
It's obvious that the regular chain $\{H_9, H_6, g_{1a}\}$ is squarefree.
Using the \textit{RealRootIsolate} program in Maple and taking the accuracy being $1/10^{20}$, we find that
the solution set of the semi-algebraic system $S_a$ is null, implying that $\textbf{V}(g_{1a}, H_6, H_9)=\emptyset$.  Hence, $\textbf{V}(\varphi_1, \varphi_2, \varphi_3, g_{1a},
g_{1b})=\textbf{V}(\varphi_2, \varphi_3, H_3)\cap\textbf{V}(g_{1a}, H_6, H_9)=\textbf{V}(\varphi_2, \varphi_3, H_3)\cap\emptyset=\emptyset$.
This leads to a contradiction and proves our assertion.

To sum up, we know that $V_{2i+1}=0$, $i=1, 2, 3$, if and only if
\begin{equation}\label{dai4.7}
(K, a, b)=(K_{2}, a_{4}, b_{2}).
\end{equation}
By Theorem 2.3 of \cite{Lu2007},  we have
\begin{equation*}
\varphi_4(K_{2}, a_{4}, b_{2})\leq {\varphi}^+_{4}(\overline{K}_{2}, \overline{a}_{4}, \overline{b}_{2})+{\varphi}^-_{4}(\underline{K}_{2}, \underline{a}_{4}, \underline{b}_{2})\approx -{
1.725035818\times 10^{103}}<0,
\end{equation*}
which yields $V_9(K_2, a_4, b_2)<0$.
This implies that the unique positive equilibrium $E^{*}(1, 1)$ of system \eqref{norm} is not a center but a weak focus of order at most $4$, provided that it is center or focus type.
Furthermore, it is a weak focus of order $4$ if and only if $s=s^{*}$ and \eqref{dai4.7} hold.
By Theorem 2.3.2 of \cite{Han2013}, at most $4$ limit cycles can bifurcate from $E^{*}(1, 1)$. This proves the assertions (1) and (2).

Next we will prove the assertion (3).
Similar to determining the sign of $\varphi_4(K_{2}, a_{4}, b_{2})$, we can check by Maple that the Jacobian determinant of $\varphi_1$, $\varphi_2$,  $\varphi_3$ with respect to $K$, $a$, $b$
at $(K_2, a_4, b_2)$ is negative, which means that the Jacobian matrix of  $V_1$, $V_3$, $V_5$,  $V_7$ with respect to $s$, $K$, $a$, $b$ has its full rank $4$, i.e.,
\begin{eqnarray}\label{dai4.8}
\textrm{rank}\left[\frac{\partial (V_1, V_3, V_5, V_7)}{\partial (s, K, a, b)}\right]_{(\hat{s}, K_2, a_4, b_2)}=4,
\end{eqnarray}
where $\hat{s}=s^{*}(K_2, a_4, b_2)$.
By Theorem 2.3.2 of \cite{Han2013}, there exist some parameter values such that system \eqref{norm} has $4$ small limit cycles around  $E^{*}(1, 1)$
and two of them are stable (see Fig. 3.1. (a)).
Obviously, we can choose appropriate parameter perturbations such that system \eqref{norm} has $i$ small limit cycles around  $E^{*}(1, 1)$ for each $i=1, 2, 3$. This completes the proof.

\subsection{Proof of Theorem \ref{T4.3}}

If $\Delta<0$ and $-2\sqrt{a}<b<K a-K m h/r$, then the equation \eqref{dai1.2} has three different positive roots in the interval $(0, K)$.
Correspondingly, system \eqref{model} has three different positive equilibria, i.e., system \eqref{norm} has three distinct positive equilibria.
By Lemma \ref{L2.2}, we know that  one of the three positive equilibria is a saddle and the others are anti-saddles.
And we have $\mu\in \Sigma_2$, where
\begin{equation}\label{cond4.9}
\Sigma_2=\left\{\mu\in\mathbb{R}^3\left|K>1,\, a>0,\, -2\sqrt{a}<
b<\dfrac{1}{K}-1-2\sqrt{\dfrac{a}{K}},\, b\neq \frac{1}{K}-2-\frac{a}{K}\right.\right\}.
\end{equation}
Without loss of generality, assume that $E^{*}(1, 1)$ is one of the two positive anti-saddles.
If $E^{*}(1, 1)$ is center or focus type, then $s=s^{*}$ and $\mu\in \Lambda$.
Moreover, further assume that the value $1$ is the minimal root of  Eq. \eqref{dai2.13} in $(0, K)$, whereas the proof of the case that $1$ is the maximum  root of  Eq. \eqref{dai2.13} is
similar.
That is, we have $(1-Kb-K-\sqrt{\bar{\Delta}})/(2a)>1$, which implies $\mu\in \Sigma_3$, where
\begin{equation}\label{cond2.30}
\Sigma_3=\left\{\mu\in\mathbb{R}^3\left|K b+2 K+a-1>0, \, K b+K+2 a-1<0\right.\right\}.
\end{equation}

To prove the assertion (1), it suffices to prove the first three Lyapunov constants  in \eqref{dai2.19}  have no common real root, i.e., $\{\varphi_1, \varphi_2, \varphi_3\}$  have no common
real root,  such that $\mu\in \Lambda\cap\Sigma_2\cap\Sigma_3$.
Using the same arguments as the proof of Theorem \ref{T4.1}, we only need to check whether the parameter values $(K_2, a_4, b_2)$ locate in $\Lambda\cap\Sigma_2\cap\Sigma_3$.
Since $b_2>0$, we have $(K_2, a_4, b_2)\not\in\Sigma_2$.
This implies that the center or focus type equilibrium $E^{*}(1, 1)$ is not a center but a weak focus of order at most $3$. By Theorem 2.3.2 of \cite{Han2013}, at most $3$ limit cycles can
bifurcate from it.

To prove the assertion (2), we first find some parameter values such that $E^{*}(1, 1)$ is an order two weak focus.
To this end, we set $K=100$, $a=60$  and hence $\varphi_1$ is a polynomial of $b$. By command $\textsl{realroot}$ with accuracy $1/10^{20}$ in Maple and Theorem 2.3 of \cite{Lu2007}, we can
find $\varphi_1(100, 60, b)$ has a unique real root $b^{\ast}$ such that $(100, 60, b^{\ast})\in \Lambda\cap\Sigma_2\cap\Sigma_3$, where
\begin{align*}
b^{\ast}\in\left[-{\frac {1510607606947999582739}{590295810358705651712}}, -{\frac {
3021215213895999165477}{1180591620717411303424}}
\right]\triangleq \left[\underline{b}^{\ast}, \overline{b^{\ast}}\right].
\end{align*}
Furthermore, by Theorem 2.3 of \cite{Lu2007}, we have
\begin{equation*}
\varphi_2(100, 60, b^{\ast})\leq {\varphi}^+_{2}(100, 60, \overline{b^{\ast}})+{\varphi}^-_{2}(100, 60, \underline{b}^{\ast})\approx -{1.546063838\times 10^{19}}<0,
\end{equation*}
which yields $V_5(100, 60, b^{\ast})<0$.
We first perturb $b$ small near $b^{\ast}$ such that $V_3V_5<0$ and adjust $s$ such that $V_1=0$ holds.  A limit cycle bifurcates. For the second limit cycle, perturb $s$ so that $V_1$ is of the opposite sign of $V_3$. Therefore, two limit cycles can bifurcate from $E^{*}(1, 1)$ (see Fig. 3.2). This completes the proof.

\subsection{Proof of Theorem \ref{T4.4}}

By Theorem \ref{T4.3}, the conditions in Theorem \ref{T4.4} imply that system \eqref{norm} has three distinct positive equilibria,  one of them is a saddle and the others are anti-saddles.
To simplify the calculations, we introduce the new parameters as follows.
Under  the hypothesis that system \eqref{norm} has three distinct positive equilibria,  we assume that  Eq. \eqref{dai2.13} has three positive zeros $1, \alpha$ and $\beta$ in the interval $(0,
K)$  with $1<\alpha<\beta<K$.
Then $E^{*}_2(\alpha, \alpha)$ is a saddle, $E^{*}(1, 1)$ and $E^{*}_3(\beta, \beta)$ are both anti-saddle, respectively.
From Eq. \eqref{dai2.13} and Vieta's formulas for quadratic polynomial, we have
\begin{equation*}
\alpha\beta=\dfrac{K}{a}, \,\,
\alpha+\beta=-\dfrac{Kb+K-1}{a},
\end{equation*}
i.e.,
\begin{equation}\label{dai4.11}
a=\dfrac{K}{\alpha\beta}, \,\,
b=\dfrac{1}{K}-\dfrac{1}{\alpha}-\dfrac{1}{\beta}-1.
\end{equation}
Then, from Lemma \ref{L2.2}, the conditions that system \eqref{norm} has three distinct positive equilibria are equivalent to
\begin{equation}\label{cond2.31}
1<\alpha<\beta<K,\,\,  s>0,\,\, \alpha \beta -K \alpha-K \beta -K\alpha \beta+2 K \sqrt{K\alpha\beta}>0.
\end{equation}
With a time scaling transformation $d\tau=(K-\alpha) ( K-\beta) d t/(K^2 {\alpha}^{2}{\beta}^{2})$, system \eqref{norm} which has been replaced by Eq. \eqref{dai4.11} can be reduced to the following equivalent differential system (we will still use $t$ to denote $\tau$ for ease of
notation).
\begin{equation}\label{dai4.12}
\left\{
\begin{aligned}
\dot{x}=&x^2( K-x )\left[{K}
^{2}{x}^{2}-(K\alpha\,\beta+K\alpha+K\beta-\alpha\,
\beta)x+K\alpha\,\beta\right]\\
&-(K-1)(K-\alpha)(K-\beta) {x}^{3}y
, \\[2pt]
\dot{y}=& K s y (x-y)\left[{K}
^{2}{x}^{2}-(K\alpha\,\beta+K\alpha+K\beta-\alpha\,
\beta)x+K\alpha\,\beta\right],
\end{aligned}
\right.
\end{equation}
where the parameters $K$, $\alpha$, $\beta$ and $s$ satisfy \eqref{cond2.31}.
Next we only need to consider system \eqref{dai4.12}, instead of system \eqref{norm}.

For convenience, we denote the $n$-th Lyapunov constant of system \eqref{dai4.12} at  $E^{*}(1, 1)$ and $E_3^{*}(\beta, \beta)$  by $V_{2n+1}^{(1)}$ and $V_{2n+1}^{(3)}$, respectively, $n=0, 1,
\cdots$. Then  $V_{1}^{(1)}=\textrm{Tr}(J(E^{*}))$ and $V_{1}^{(3)}=\textrm{Tr}(J(E^{*}_3))$, where $J(E^{*})$ is the Jacobian matrix of system \eqref{dai4.12} at $E^{*}$.
From system \eqref{dai4.12}, it's not difficult to obtain
\begin{equation}\label{dai4.13}
\left\{
\begin{aligned}
V_{1}^{(1)}=&-{K}^{2}\alpha\beta+{K}^{3}+K\alpha\beta-2
{K}^{2}+K\alpha+K\beta-\alpha\beta-s K (K-\alpha) (K-\beta),\\
V_{1}^{(3)}=&(K^3\beta-2K^2{\beta}^{2}+K\alpha{\beta}^{2}-{K}^{2}\alpha+K\alpha\beta+
K{\beta}^2-\alpha{\beta}^{2}){\beta}^{2}-\\
&s K{\beta}^3(K-1)(K-\alpha).
\end{aligned}
\right.
\end{equation}
If both $E^{*}(1, 1)$ and $E_3^{*}(\beta, \beta)$ are center or focus type, then $V_{1}^{(1)}=V_{1}^{(3)}$, which implies
\begin{equation}\label{dai4.14}
\left\{
\begin{aligned}
&s={\dfrac {{K}^{3}\beta-2{K}^{2}{\beta}^{2}+{K}^{3}-2{K}^{2}\beta+3
K{\beta}^{2}-2{K}^{2}+3 K\beta-2\,{\beta}^{2}}{K \left( (\beta+1) {K}^{2}
-4\,\beta\,K+{\beta}^{2}+\beta \right) }}\triangleq s_0, \\
&\alpha={\dfrac {K \beta ( 3 {K}^{2}-2 K\beta-2 K+\beta ) }{
(K \beta+K-\beta) ( {K}^{2}-\beta ) }}\triangleq \alpha_0.
\end{aligned}
\right.
\end{equation}
Denote $l(K)\triangleq (\beta+1) {K}^{2}
-4\,\beta\,K+{\beta}^{2}+\beta$, where $l(\cdot)$ is a polynomial function. Since $l'(K)=2(\beta+1)K-4\beta
=2[\beta(K-1)+(K-\beta)]>0$ for $K>\beta>1$, we have $l(K)>l(1)=(\beta-1)^2>0$. Thus, the denominators of the right-hand side of Eq. \eqref{dai4.14} are both
positive.

We first prove the assertion (1). It suffices to prove $V_{3}^{(1)}>0$ and $V_{3}^{(3)}>0$.

Substituting \eqref{dai4.14}  to the system \eqref{dai4.12} and using the software Maple for symbolic calculate, we get
\begin{equation*}
V_{3}^{(1)}=\dfrac {2 {K}^{6}\beta(K-\beta)^{2}(\beta-1)^{2}F_1}{
({K}^{2}-\beta) ^{4}(K\beta+K-\beta)^{3}},\ \
V_{3}^{(3)}=\dfrac {2{K}^{6}{\beta}^{4}(K-1)^{2}(\beta-1)^{2}G_1}{({K}^{2}-\beta)^{4}
(K\beta+K-\beta)^{3}},
\end{equation*}
where $F_1$ and $G_1$ are both polynomials in $K$ and $\beta$ of degree $11$.  From \eqref{cond2.31}, it's easy to see that the denominators of $V_{3}^{(1)}$ and $V_{3}^{(3)}$ are both positive
and the signs of $V_{3}^{(1)}$ and $V_{3}^{(3)}$ are the same as that of $F_1$ and $G_1$, respectively. Hence we only need to prove $F_1>0$ and $G_1>0$.

Denote $\epsilon\triangleq (K-\beta)/(\beta-1)>0$, then $\beta=(\epsilon+K)/(\epsilon+1)$. Substituting it to $F_1$, we obtain
$F_1=(K-1)^{6}\overline{F}_1/(\epsilon+1)^{6}$, where $\overline{F}_1$ is collected by $\epsilon$ as follows.
\begin{eqnarray*}
\begin{aligned}
\overline{F}_1=&(24{K}^{3}-4K+4
){\epsilon}^{6}+(24{K}^{4}+45{K}^{3}-21{K}^{2}+7
K+5) {\epsilon}^{5}+(6{K}^{5}+75{K}^{4}-\\
&
4{K}^{3}-11{K}^{2}+12K+2){\epsilon}^{4}+(33{K}^{5}+42{K}^{4}-18{K}^{3}+3{K}^{2}+6K) {\epsilon}^{3}+\\
&(29{K}^{5}+7{K}^{4}-9{K}^{3}+6{K}^{2}) {\epsilon}^{2}+(11{K}^{5}-4{K}^{4}+2{K}^{3})\epsilon+{K}^{5}.
\end{aligned}
\end{eqnarray*}
It's easy to see that all the coefficients of the power of $\epsilon$ in $\overline{F}_1$ are positive since $K>1$. Noting that $\epsilon>0$, we have $\overline{F}_1>0$, which yields $F_1>0$.

Using the same arguments as the proof of $F_1>0$, we can prove $G_1>0$. This completes the proof of the assertion (1).

Next we will prove the assertion (2). For any given values of $K$ and $\beta$, let $(s, \alpha)=(s_0+\varepsilon, \alpha_0)$, where $s_0$, $\alpha_0$ are given by Eq. \eqref{dai4.14}, and
$\varepsilon$ is a perturbation parameter. From \eqref{dai4.13}, we have
\begin{equation}\label{dai4.15}
L_1=-K ( K-\beta)(K-\alpha_0)\varepsilon,\ \
T_1=-K{\beta}^{3} ( K-1)( K-\alpha_0)\varepsilon.
\end{equation}
By continuity, we can choose $\varepsilon$ small enough such that $V_{3}^{(1)}>0$ and $V_{3}^{(3)}>0$.
From \eqref{dai4.15},  it follows from $K>\beta>\alpha_0>1$ that $V_{1}^{(1)}<0$ and $V_{1}^{(3)}<0$ hold as long as $\varepsilon>0$.
That is to say, there exist parameter values $(s, \alpha)$ near $(s_0, \alpha_0)$ such that $L_1V_{3}^{(1)}<0$ and $T_1V_{3}^{(3)}<0$. This implies that one limit cycle can bifurcate from
$E^{*}(1, 1)$ and $E_3^{*}(\beta, \beta)$, respectively. This ends the proof.

\section{Discussion}

In this paper, we studied the Hopf bifurcation and global dynamics of system \eqref{model} with generalized Holling functional response of type III when it has no degenerate positive equilibrium.
We not only rigorously proved the existence of some phenomena that have been observed in \cite{Hsu1995,Huang2014}, but also observed some new and complicated dynamical behaviors,
such as the coexistence of four limit cycles including two stable cycles, three limit cycles including two stable cycles, one  big stable limit cycle enclosing two unstable limit cycles, or three hyperbolic positive equilibria and two limit cycles surrounding one of them.
Therefore, our results can be a complement to the works by Hsu et al \cite{Hsu1995} and  Huang et al \cite{Huang2014} for this model, and show that the maximum number of limit cycles in the
model of this kind is at least $4$, which improves the preceding results that this number is at least $2$.
Our results also indicate that the nonlinear dynamics of such biological and epidemiological models not only depend on more bifurcation parameters but also are very sensitive to parameter perturbations, which are important for the control of biological species or infectious diseases.
Furthermore, the coexistence of bistable states (one stable limit cycle and one stable equilibrium, or two stable limit cycles) and tristable states (two stable limit cycles and one stable equilibrium, or two stable equilibria  and one stable limit cycle) show that the model of this kind is also highly sensitive to initial values.

In this work, the coexistence of two stable limit cycles gives a positive answer to one of the almost impossible projects proposed in \cite{Coleman1983} which is given in the first section.
This shows the ecological relevance of the coexistence of multiple limit cycles in
predator-prey systems and the importance of our results, which should serve for the outlined
problem to be actually feasible in a biological lab with appropriate little creatures \cite{Coleman1983}. It will be interesting to see if two stable limit cycles occurs in realistic predator-prey systems.

This paper provided some available methods for studying the existence, stability, number and distribution of limit cycles in a given predator-prey systems.
It is worth mentioning that these methods can be applied to other models to study the Hopf bifurcation at positive equilibria with complicated expressions or no explicit expressions. However, there are two interesting problems that remain open. One is to prove the existence of the big stable limit cycle for some given parameter values, see Fig 1.1. If we can rule out the possibility  that there are  two homoclinic loop connecting with the saddle $E_2^*$, called eight-loop in general, then Poincar\'{e}-Bendixson Theorem implies the existence of the big stable limit cycle. We conjecture that this
problem may be solved by constructing an appropriate inner boundary of an annular region. The other is whether $3$ limit cycles can bifurcate from an arbitrary anti-saddle positive equilibrium if system has three distinct positive equilibria, see Remark \ref{rem3}. This problem may be very challenging since both the conditions and the expressions of $V_3$ and $V_5$ are extremely complicated.
Also, is there a natural (or laboratory) predator-prey system with at least three ecologically stable cycles?

config: arxiv


TITLE: Is it possible to solve $ax^2+hxy+by^2+c=0$ in integers?
QUESTION [0 upvotes]: Last time I got stuck in this problem which I have posted earlier. Today I have come accross to this new situation.
How to solve the diophantine equation $ax^2+hxy+by^2+c=0$ in integers ? Given all of $a,b,c,h\in \mathbb Z$.
The motivation was to solve $30x^2+21y^2-57xy+729=0$ in integers which through MAPLE i got as $(x,y)=(-22,-17), (-10,-11), (10,11), (22,17)$.

REPLY [1 votes]: In general, given an initial solution $x_0,y_0 = m,n$ to,
$$a x^2 + b x y + c y^2 + d = 0\tag1$$
if the discriminant $D= b^2-4ac\,$ is a non-square $D>0$, then an infinite more can be found as,
$$x = mu^2 - 2(b m + 2c n)u v + D m v^2$$
$$y = n u^2 + 2(2a m + b n)u v + D n v^2$$
where $u,v$ solve the Pell equation,
$$u^2 - D v^2 = \pm1\tag2$$
However, if your $D$ is a square, or $D<0$, then $(1)$ will only have a finite number of integer solutions. (In your case, it is a square $D=27^2$.)


TITLE: given a infinite dimensional vector space over Complex field, can you find a list of linear independent list of any length?
QUESTION [1 upvotes]: I'm trying to prove the following question with a different path,
Suppose there exists a linear map on V whose null space and range are both finite-dimensional, prove that V is finite-dimensional.
Now, instead of trying to build a basis from null space and rang T, I want to say that since V is infinite-dimensional, I can construct a list of linearity independent vectors of any length k>0. Specifically, $k>\dim $ range$ T =n$, now suppose $v_1,...,v_k$ are such linearity independent vectors, then $Tv_1,\dots, Tv_k$ are linearly independent in the range, but this is a contradiction since we can't have a length of independent list bigger than the dimension of the space.

REPLY [1 votes]: I found the problem interesting and I think I found a proof.
Firstly, I have been told that this theorem only holds for vector spaces over the same field. So the question becomes the following.
Let $V$ and $W$ be vector spaces defined over the same field $F$. Let $T : V\to W$ be a linear map. Suppose that range$T$ and null$T$ are finite dimensional. Prove that $V$ is finite dimensional.
Proof: Given that range$T$ is finite dimensional, there exists a basis of range$T$. Let $Tv_1,. . .,Tv_n$ be the basis of range$T$. Assuming that the function $T$ is not a partial function. We see that for every $v\in V$, $Tv\in$ range$T$. Because $Tv_1,. . .,Tv_n$ is a basis, there exists scalars $a_1,. . .,a_n$ such that $Tv=a_1Tv_1+. . .+a_nTv_n$. Subtracting $Tv$ from both sides and using the linearity of $T$ we get $T(a_1v_1+. . .+a_nv_n-v)=0$.
Given that null$T$ is finite dimensional. Let $w_1,. . .,w_m$ be a basis of null$T$. Then, $a_1v_a+. . .+a_nv_n-v=c_1w_1+. . .+c_mw_m$ $\implies$ $a_1v_a+. . .+a_nv_n-c_1w_1+. . .+c_mw_m=v$.
The list $v_1,. . .,v_n,w_1,. . .,w_m$ spans $V$ because $v$ was an arbitrary vector. Hence, $V$ is finite dimensional.
Note that this result is also consistent with the rank-nullity theorem. Note that you cannot use the rank nullity theorem to prove this as the rank nullity theorem already assumes that $V$ is finite dimensional.
I did post this proof as another question. Although I think the proof is correct. If someone finds a mistake, I will update you on it.


text:
\section{Statistical verification and approximation of safe sets}

In order to validate our deep neural network controllers considering safety we compared them to the benchmark explicit MPC in terms of constraint satisfaction along state trajectories $\xi = \{x_0,\dots,x_{k_{\text{end}}}\}$.
We ran closed-loop simulations from $n_{\text{r}}=20000$ random initial points $\hat{x}_{i}$ for each of the both examples with the optimal and approximated method and obtained the datasets $\mathcal{D}_{\text{dnn}} = \{\xi^{\text{dnn}}_i,\dots,\xi^{\text{dnn}}_{n_{\text{r}}}\}$ and
$\mathcal{D}_{\text{exp}} = \{\xi^{\text{exp}}_i,\dots,\xi^{\text{exp}}_{n_{\text{r}}}\}$.
These sets of state-trajectories were classified by the logical requirement
\begin{align}
\gamma(\xi) = \square_{\interval{0}{k_{\text{end}}}}(|C_{xu}x_k| \leq c_{xu}), \label{eq:log_requirement}
\end{align}
where $C_{xu} = [C_x^T, C_u^T]^T$ and $c_{xu} = [c_x, c_u]^T$ define the constraints including states.
This means that a closed-loop simulation run is assessed as valid only if no constraint violations occur.
In the following, sets which only contain valid initial values/state trajectories will be marked in the exponent with a '$+$' and those containing invalid ones with a '$-$'.
Sets denoted by $\hat{\mathcal{D}}_{(\cdot)}^{(\cdot)} = \{\hat{x}_{i}^{(\cdot)},\dots,\hat{x}_{i_{\text{end}}}^{(\cdot}\}$ contain the initial values $\hat{x}_i^{(\cdot)}\}$ of trajectories $\xi_i^{(\cdot}\}$ contained in sets denoted by $\mathcal{D}_{(\cdot)}^{(\cdot)} = \{\xi_{i}^{(\cdot)},\dots,\xi_{i_{\text{end}}}^{(\cdot)}\}$.
For all pairs '$+$'/'$-$' following equations hold:
\begin{subequations}
	\begin{align}
		\hat{\mathcal{D}} = \hat{\mathcal{D}}^+ \cup \hat{\mathcal{D}}^-, \\
		\hat{\mathcal{D}}^+ \cap \hat{\mathcal{D}}^- = \emptyset.
	\end{align}
\end{subequations}

Our estimations of the volume will always be point-wise measure with respect to $|\hat{\mathcal{D}}_{\text{exp}}^+|$, which is the cardinality of all initial values of the valid trajectories generated by the explicit MPC solution.
For finding explicit expressions for safe sets for our approximated solution, we split the set $\hat{\mathcal{D}}_{\text{dnn}}$ into $\hat{\mathcal{D}}_{\text{s}}$ for computing the safe sets and $\hat{\mathcal{D}}_{\text{t}}$ for testing them.

\subsection{Implicit description}

In case of an optimal approximation $\mathcal{L}_7$/$\mathcal{L}_{10}$ and  $\mathcal{N}_{6,6}$/$\mathcal{N}_{10,6}$ would be mathematically equivalent.
This would mean all initial values $\hat{x}_i$ which fulfill \eqref{eq:log_requirement} are part of the safe set.
Since $\mathcal{N}_{6,6}$/$\mathcal{N}_{10,6}$ are only approximations, we compute the percentage of points within the safe set which can be covered by our approximation.
This set contains the initial values $\hat{x}_{i}$ of trajectories which did not lead to constraint violations throughout the closed-loop run.
To compare the performance of the approximated controller we define accordingly the subset $\hat{\mathcal{D}}_{\text{dnn}}^+$.
This allows us to obtain a first naive measure $m_{\text{dir}}$ of our safe set by computing the ratio of the cardinality of the two sets:
\begin{align}
m_{\text{dir}} = \frac{|\hat{\mathcal{D}}_{\text{dnn}}^+|}{|\hat{\mathcal{D}}_{\text{exp}}^+|}. \label{eq:m_dir}
\end{align}
For both examples the ratio is $m_{\text{dir}} = \SI{97.7}{\percent}$ as pointed out in the first row of Tab.~\ref{tab:safe_set_and_approximation}.
This shows that the performance of the approximated controllers is close to the optimal solution.