\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{amssymb}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 255, pp. 1--20.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2016/255\hfil Density-dependent predator-prey systems]
{Stability analysis and Hopf bifurcation
 of  density-dependent predator-prey systems
 with Beddington-DeAngelis functional response}

\author[X. Jiang, Z. K. She, Z. Feng \hfil EJDE-2016/255\hfilneg]
{Xin Jiang, Zhikun She, Zhaosheng Feng}

\address{Xin Jiang \newline
 School of Mathematics and Systems Science, Beihang University,
Beijing 100191, China}
\email{jiangxin1991@126.com}

\address{Zhikun She (corresponding author) \newline
School of Mathematics and Systems Science,
Beihang University, Beijing 100191, China}
\email{zhikun.she@buaa.edu.cn}

\address{Zhaosheng Feng \newline
School of Mathematical and Statistical Sciences,
University of Texas-Rio Grande Valley,
Edinburg, TX 78539, USA}
\email{zhaosheng.feng@utrgv.edu}

\thanks{Submitted March 6, 2016. Published September 21, 2016.}
\subjclass[2010]{34D20, 34E05, 37G15}
\keywords{Density-dependent; local and global stability; Hopf bifurcation;
\hfill\break\indent  monotonicity; geometrical restriction}

\begin{abstract}
 In this article, we study a density-dependent predator-prey
 system with the Beddington-DeAngelis functional response for stability
 and Hopf bifurcation under certain parametric conditions.
 We start with the condition of the existence of the unique
 positive equilibrium, and provide two sufficient conditions
 for its local stability by the Lyapunov function method and the Routh-Hurwitz
 criterion, respectively. Then, we establish sufficient conditions for
 the global stability of the positive equilibrium by proving the
 non-existence of closed orbits in the first quadrant $\mathbb{R}^2_{+}$.
 Afterwards, we analyze the Hopf bifurcation geometrically by
 exploring the monotonic property of the trace of the Jacobean matrix with
 respect to $r$ and analytically verifying that there is a unique $r^{*}$
 such that the trace is equal to $0$.
 We also introduce an auxiliary map by restricting all the five parameters
 to a special one-dimensional geometrical structure and analyze the Hopf
 bifurcation with respect to all these five parameters. Finally, some
 numerical simulations are illustrated which are in agreement with our
 analytical results.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{example}[theorem]{Example}
\allowdisplaybreaks

\section{Introduction}

The dynamic relationship between predators and their preys is
one of the dominant themes in both mathematical biology and
theoretical ecology. Better understanding exact trends of
population dynamics can contribute to the environmental protection
and resource utilization. In the past decades, considerable attention
has been dedicated to various predator-prey
models \cite{Method1, G1,du1,du2,Method2,L15,G2, Method3,G3,G4,G5},
of which the following predator-prey system with the
Beddington-DeAngelis functional response,
originally proposed by Beddington \cite{b} and DeAngelis
\cite{d}, independently, has been extensively studied by applied mathematicians
and biologists theoretically and experimentally:
\begin{equation}
\begin{gathered}
\frac{dx(t)}{dt}=x(t) \Big(c-bx(t)-\frac{sy(t)}{m_1+m_2x(t)+m_3y(t)}\Big),
\\
\frac{dy(t)}{dt}=y(t) \Big(-d+\frac{fx(t)}{m_1+m_2x(t)+m_3y(t)} \Big),
\end{gathered}\label{1.1}
\end{equation}
where $c,b,s,m_1,m_2,m_3,d$ and $f$ are positive constants,
$x(t)$ and $y(t)$ represent the population
density of the prey and the predator at time $t$ respectively, $c$ is the
intrinsic growth rate of the prey, $d$ is the death rate of
the predator, and $b$ stands for the mutual interference between preys.
The predator consumes the prey with the functional response of
Beddington-DeAngelis type
$\frac{sx(t)y(t)}{m_1+m_2x(t)+m_3y(t)}$ and contributes to its
growth with the rate $\frac{fx(t)y(t)}{m_1+m_2x(t)+m_3y(t)}$. This
is similar to the well-known Holling type II model with an extra
term $m_3y(t)$ in the denominator.

System \eqref{1.1} is a population model which has
received much attention in biology and ecology.
Cantrell and Cosner \cite{cc} presented qualitative
analysis of system \eqref{1.1} on permanence, which implies the
existence of a locally asymptotically stable positive equilibrium
or periodic orbits. Hwang \cite{first} considered
the local and global asymptotic stability of the positive
equilibrium $(\hat{x},\hat{y})$ by the divergence criterion.
That is, for system \eqref{1.1}, under the conditions
$(f-dm_2)\frac{c}{b}>dm_1$ and $tr(J(\hat{x},\hat{y}))\leq0$,
$(\hat{x},\hat{y})$ is globally asymptotically stable.
Recently, Hwang \cite{second} presented the condition
$(f-dm_2)\frac{c}{b}>dm_1$
and $tr(J(\hat{x},\hat{y}))>0$ to ensure the
uniqueness of the limit cycle of system \eqref{1.1}.
For more details about biological background of system \eqref{1.1}, we refer
the reader to \cite{Bio1,cc,Bio2,d,first}.


However, abundant evidence suggests that predators do interfere
with each other's activities so as to result in competition efforts
and that the predator may be of density dependence because of the
environmental factors \cite{bs,bss}.
Kratina et al have demonstrated a fact
that the predator density dependence is significant at both
high predator densities and low predator densities \cite{pmab}.
Hence, the following more realistic predator and prey density-dependent
model with the Beddington-DeAngelis functional response
was proposed \cite{L14,L11}:
\begin{equation}
\begin{gathered}
\frac{dx(t)}{dt}=x(t) \Big(c-bx(t)-\frac{sy(t)}{m_1+m_2x(t)+m_3y(t)} \Big),\\
\frac{dy(t)}{dt}=y(t) \Big(-d-ry(t)+\frac{fx(t)}{m_1+m_2x(t)+m_3y(t)} \Big),
\end{gathered}\label{1.2}
\end{equation}
where $r$ represents the rate of predator density dependence. Compared
with system \eqref{1.1}, system \eqref{1.2} contains not
only $bx^2(t)$ which stands for intraspecific action of prey
species, but also $ry^2(t)$ which stands for intraspecific action
of predator species.


Li and She \cite{L14} considered dynamics of system
\eqref{1.2} by showing $(f-dm_2)\frac{c}{b}>dm_1$ as
the sufficient and necessary condition for the permanence and
existence of the unique positive equilibrium, and
$(f-dm_2)\frac{c}{b}\leq dm_1$
as the sufficient and necessary condition for the
global asymptotical stability of boundary solution.
Furthermore, the dynamics of the stage-structured model
of system \eqref{1.2} was studied in \cite{L13} and
the non-autonomous case of system \eqref{1.2} was tackled in \cite{L1502}.

From the point of view of ecological managers, it may be desirable
to have a unique positive equilibrium which is globally
asymptotically stable.
Although considerable attention has been undertaken on system \eqref{1.2},
it seems that sufficient and necessary conditions for local stability and
even global stability of the
positive equilibrium have not been comprehensively presented yet. In addition,
the existence of (stable or unstable) limit cycles and
Hopf bifurcation with respect to the parameters are rarely discussed.


In this article, we will discuss the local and global
stability of the positive equilibrium and analyze the Hopf bifurcation
of system \eqref{1.2} in the first quadrant. For simplicity, by the re-scaling
$ t \to c t$,
$x \to \frac{b}{c}x$,
$y \to \frac{bm_3}{c m_2}y$,
$s = \frac{s}{c m_3}$,
$a = \frac{bm_1}{c m_2}$,
$b = \frac{f}{c m_2}$,
$d = \frac{m_2d}{f}$
and $r = \frac{c rm_2^2}{bfm_3}$,
system \eqref{1.2} is nondimensionalized to
\begin{equation}
\begin{gathered}
\frac{dx(t)}{dt}=P(x(t),y(t)):= x(t)(1-x(t))-\frac{sx(t)y(t)}{x(t)+y(t)+a},\\
\frac{dy(t)}{dt}=Q(x(t),y(t)):= by(t)
 \Big(-d-ry(t)+\frac{x(t)}{x(t)+y(t)+a}\Big).
\end{gathered}
\label{1.3}
\end{equation}

We start with the existence of the unique
positive equilibrium of system \eqref{1.3}
and show that this unique positive equilibrium cannot be a saddle under the
given conditions. To ensure the local stability of this positive equilibrium,
we first provide a sufficient condition by constructing a Lyapunov function,
and then give a concrete condition depending only on the parameters by the
Routh-Hurwitz criterion.
By virtue of two classical criteria, we discuss the global asymptotic stability
of the positive equilibrium and present several nonequivalent sufficient conditions.
That is, under the permanence condition, we firstly present a condition by Dulac's
criterion for the global attractiveness of the positive equilibrium.
Secondly, by the divergency criterion, we provide a sufficient condition for the
stability of all possibly existing closed orbits, under which we prove that
the positive equilibrium is locally and globally asymptotic stable because
of the non-existence of stable closed orbits.
Thirdly, under a stronger permanence condition for providing concrete bounds
of all possible closed orbits, we apply Grammer's
rule and Green's theorem to present the divergency integral,
and obtain another sufficient condition on global asymptotical stability
of the positive equilibrium.


Afterwards, we explore the Hopf bifurcation with respect to the parameter $r$.
We first geometrically explore the monotonic property of the trace of the
Jacobian matrix with respect
to $r$, and then analytically verify that there is a unique $r^{*}$ such that
the trace is equal to $0$. Based on these arguments, we analyze the Hopf
bifurcation with respect to the parameter $r$.
Moreover, in order to generalize our conclusion on the Hopf bifurcation,
 we consider the Hopf bifurcation with respect to all five parameters by
introducing an auxiliary map and restrict
the five parameters to a special one-dimensional geometrical structure.
Some numerical simulations are performed to
illustrate our analytical results.

Note that Hwang \cite{first} verified that for system \eqref{1.1}, the local
stability and global stability of the positive equilibrium coincide. However,
from the Hopf bifurcation analysis, we find that for system \eqref{1.3},
the coincidence between local stability and global stability does not hold
due to the existence of the parameter $r$.
Consequently, the analysis results show that the rate $r$ of predator density
dependence has a significant
effort on the global dynamics of system \eqref{1.3}.

The rest of this paper is organized as follows.
In Section 2, we consider the local stability of
the positive equilibrium.
In Section 3, we establish several sufficient conditions for the global stability of
the positive equilibrium by two classical criteria.
In Section 4, we analyze Hopf bifurcations with
respect to the parameter $r$ and all five parameters, respectively.
Section 5 presents some numerical simulations and
Section 6 is a brief conclusion.

\section{Local stability of positive equilibrium}
\label{sec:1}

We know from \cite{L14} that system \eqref{1.3} has a unique positive equilibrium
if and only if the condition
\begin{equation}
 ad < 1-d \label{H1}
\end{equation}
holds, which is the sufficient and necessary condition for permanence of the system.
This unique positive equilibrium cannot be a saddle point.
According to the definition of permanence \cite{L14} and the Poinc\'{a}re-Bendixson
theorem \cite{DC}, for any trajectory starting in
$\mathbb{R}^2_{+}:=\{(x,y)\in\mathbb{R}^2:x>0,y>0\}$,
there are three cases for its $\omega$-limit set in $\mathbb{R}^2_{+}$:
the positive equilibrium which is not a saddle,
a closed orbit or a saddle together with possible homoclinic orbits. For the third
case, in addition to this saddle, there exists at least another positive equilibrium
inside the region enclosed by the homoclinic orbit simultaneously.

Since the unique positive equilibrium is not a saddle,
we now explore the conditions to guarantee the local
asymptotic stability of this positive equilibrium.
Let $(x^{*},y^{*})$ be the unique positive equilibrium,
which is short for the expression
$(x^{*}(a,b,d,r,s),y^{*}(a,b,d,r,s))$ and satisfies
\begin{equation}
\begin{gathered}
1-x^{\ast}-\frac{sy^{\ast}}{x^{\ast}+y^{\ast}+a}=0,\\
-d-ry^{\ast}+\frac{x^{\ast}}{x^{\ast}+y^{\ast}+a}=0.
\end{gathered} \label{2.1}
\end{equation}
Clearly, $0<x^*<1$ and $0<y^*<\frac{1-d}{r}$.
By \cite[Theorem 1.2]{zhangzhifen}, $(x^{*},y^{*})$
smoothly depends on the parameters $a,s,r$ and $d$.
Let $x(t)\to x^{*}+x(t)$ and $y(t)\to y^{*}+y(t)$.
Then system \eqref{1.3} is equivalent to
\begin{equation}
\begin{gathered}
\frac{dx}{dt}=F^{*}_{x}x^{*}x+F^{*}_{y}x^{*}y+g^1(x,y),\\
\frac{dy}{dt}=G^{*}_{x}y^{*}x+G^{*}_{y}y^{*}y+g^2(x,y),
\end{gathered}\label{2.2}
\end{equation}
where
\begin{equation}
\begin{gathered}
F^{*}_{x}=\frac{sy^{*}}{(x^{*}+y^{*}+a)^2}-1,\quad
 F^{*}_{y}=-\frac{s(x^{*}+a)}{(x^{*}+y^{*}+a)^2},\\
G^{*}_{x}=\frac{b(y^{*}+a)}{(x^{*}+y^{*}+a)^2},\quad
 G^{*}_{y}=-\frac{bx^{*}}{(x^{*}+y^{*}+a)^2}-br,
\end{gathered}\label{2.3}
\end{equation}
and both $g^1(x,y)$ and $g^2(x,y)$ are of $o(x,y)$, given by:
\begin{gather*}
g^1(x,y)=-F^{*}_{x}x^{*}x-F^{*}_{y}x^{*}y
+\Big(x-x^2-\frac{sxy}{x+y+a}\Big),\\
g^2(x,y)=-G^{*}_{x}y^{*}x-G^{*}_{y}y^{*}y
+\Big(-bdy-bry^2+\frac{bxy}{x+y+a} \Big).
\end{gather*}

Applying the Lyapunov stability theorem \cite{wi}, we can find a condition to
ensure the local asymptotic stability of $(0,0)$ for system \eqref{2.2} as follows:

\begin{theorem} \label{Theorem2.1}
If condition \eqref{H1} holds and $F^{*}_{x}<0$, then the positive equilibrium
$(x^{*},y^{*})$ of system \eqref{1.3} is locally asymptotically stable.
\end{theorem}

\begin{proof}
Let $V(x,y)=\frac{1}{2}x^2-\frac{x^{*}F^{*}_{y}}{2y^{*}G^{*}_{x}}y^2$.
For system \eqref{2.2}, we have
\begin{align*}
\frac{dV}{dt}(x,y)
&:= \frac{\partial V}{\partial x}\frac{dx}{dt}
+\frac{\partial V}{\partial y}\frac{dy}{dt}\\
&= x^{*}F^{*}_{x}x^2-\frac{x^{*}F^{*}_{y}G^{*}_{y}}{G^{*}_{x}}y^2
+o(x^2,xy,y^2).
\end{align*}
Obviously, $V(x,y)\geq 0$ and $V(x,y)=0$ if and only if $(x,y)=(0,0)$.
Under the condition $F^{*}_{x}<0$, there exists a neighborhood $N$
of the origin such that $\frac{dV}{dt}(x,y)\leq 0$ holds in $N$ and
$\frac{dV}{dt}(x,y)=0$ if and only if $x(t)=y(t)=0$.
This implies that $V(x,y)$ is a Lyapunov function for system \eqref{2.2},
and thus $(0,0)$ is locally asymptotical stable. That is,
the positive equilibrium
$(x^{*},y^{*})$ of system \eqref{1.3} is locally asymptotically stable.
\end{proof}

\begin{remark} \label{rmk2.2} \rm
Since $|x^{*}F^{*}_{y}+y^{*}G^{*}_{x}| < \min \{-2y^{*}G^{*}_{y}, -2x^{*}F^{*}_{x} \}$
implies that $F^{*}_{x}<0$, we here obtain a weaker local asymptotic stability
condition than the one in \cite{L14}.
\end{remark}

Note that the condition of local stability in Theorem \ref{Theorem2.2} depends
on $x^{*}$ and $y^{*}$, and the positive equilibrium $(x^{*},y^{*})$
usually needs to be solved numerically at first.
It is evident that the condition $F^{*}_{x}<0$
can be directly derived by $s\leq 2a$. So, in the following,
we attempt to seek a weaker condition than $s\leq 2a$ by
applying the Routh-Hurwitz criterion for the linerization of
system \eqref{2.2} instead of constructing a Lyapunov function.

Clearly, the characteristic equation of the linerization of system
\eqref{2.2} is
\begin{equation}
\lambda^2+a_1\lambda+a_2=0, \label{2.4}
\end{equation}
where
\begin{equation}
\begin{gathered}
a_1= x^{*}+bry^{*}+\frac{(b-s)x^{*}y^{*}}{(x^{*}+y^{*}+a)^2}, \\
a_2= \big[r+\frac{x^{*}-rsy^{*}}{(x^{*}+y^{*}+a)^2}+\frac{as}{(x^{*}
+y^{*}+a)^{3}}\big]bx^{*}y^{*}.
\end{gathered} \label{2.5}
\end{equation}

Since $a_1 \leq 0$ is equivalent to
 $x^{*}+bry^{*}+\frac{(b-s)x^{*}y^{*}}{(x^{*}+y^{*}+a)^2} \leq 0$,
it follows from the first equation in \eqref{2.1} that
\[
a_1 \leq 0 \Leftrightarrow
\frac{s}{1-x^{*}} \leq \frac{(s-b)x^{*}(x^{*}+s-1)}{s(x^{*}+bry^{*})(x^{*}+a)}.
\]
Because of the relation $s>\frac{sy^{*}}{x^{*}+y^{*}+a}=1-x^{*}$,
$a_1 \leq 0$ implies that $s>b$ and $s>1+a$.
Thus, if $s\leq b$ or $s \leq 1+a$, then $a_1>0$,
i.e. the trace of the Jacobian matrix at $(x^*,y^*)$ is negative.

From \eqref{2.1}, we have
$x^{*}-rsy^{*}=1+sd-\frac{s(x^{*}+y^{*})}{x^{*}+y^{*}+a}>1+sd-s$.
According to
\[
a_2 > 0
\Leftrightarrow \big[r+\frac{x^{*}-rsy^{*}}{(x^{*}+y^{*}+a)^2}
+\frac{as}{(x^{*}+y^{*}+a)^{3}}\big]bx^{*}y^{*}>0,
\]
then $r+\frac{1+sd-s}{(x^{*}+y^{*}+a)^2}+\frac{as}{(x^{*}+y^{*}+a)^{3}}>0$
implies $a_2 > 0$. Since
 $r+\frac{1+sd-s}{(x^{*}+y^{*}+a)^2}+\frac{as}{(x^{*}+y^{*}+a)^{3}}>0$
can be derived by $s\leq1+sd+ra^2$,
we see that, if $s\leq 1+sd+ra^2$, then $a_2>0$.
Especially, when $s\leq \frac{1}{1-d}$, then $a_2>0$, i.e. the
determinant of the Jacobian matrix at $(x^*,y^*)$ is positive.


Let $S_1=\max \{2a,b,1+a\}$ and $S_2=\max \{2a,\frac{1+ra^2}{1-d}\}$.
By the Routh-Hurwitz criterion, it is straightforward to reach a condition,
depending only on parameters,
for the local stability of $(x^*,y^*)$ as follows:

\begin{theorem} \label{Theorem2.2}
 If condition \eqref{H1} and
 \begin{equation}
 s\leq \min\{S_1,S_2\} \label{H2}
\end{equation}
 hold, then the unique positive equilibrium $(x^{*},y^{*})$ of
 system \eqref{1.3} is locally asymptotically stable.
 \end{theorem}

\begin{remark} \label{rmk2.4} \rm
In \cite{L11}, it shows that for system \eqref{1.3},
the positive equilibrium is locally asymptotically stable if
$d<1\wedge ad<(1-d)\left(1-s-\frac{d}{r} \right)$ or
$s+\frac{d}{r}<1\wedge ad<(1-d)\left(1-s-\frac{d}{r} \right)$.
Clearly, this condition can also generate condition \eqref{H2}, but
\eqref{H2} looks succinct.
\end{remark}

\begin{remark} \label{rmk2.5} \rm
Notice that the unique positive equilibrium is not a saddle under
condition \eqref{H1} and we have derived that when $s\leq \frac{1}{1-d}$
(and even $s\leq S_2$), then $a_2>0$, which will be used for analyzing Hopf
bifurcation in Section 4. Here it is still of our interest whether $a_2 > 0$
 always holds or not while we analyze Hopf bifurcation.
\end{remark}


\section{Global stability of the positive equilibrium}
\label{sec:2}

In the previous section, we have presented conditions for local
asymptotic stability of the positive equilibrium. In this section,
we try to establish sufficient conditions for its global asymptotic stability.
For this, by \cite[Theorem 4.3]{L14}, we need to provide conditions
to ensure that there is no closed orbit except for the positive
equilibrium $(x^{*},y^{*})$ in the first quadrant $\mathbb{R}^2_{+}$.
We will use Dulac's criterion and the divergency criterion
for analyzing the global attractiveness of $(x^{*},y^{*})$, respectively.

\begin{theorem} \label{Theorem3.1}
For system \eqref{1.3}, if condition \eqref{H1} and
\begin{equation}
 \min \{ bd, 2a-b \} \geq 1 \label{H3}
\end{equation}
 hold, then the positive equilibrium is
globally asymptotically stable.
\end{theorem}

\begin{proof}
For system \eqref{1.3}, by choosing $B(x,y)=x+y+a$, there holds
\begin{align*}
&\frac{\partial BP}{\partial x}+\frac{\partial BQ}{\partial y}\\
&=(x+y+a)(1-bd-2x-2bry) +(1+b)x-(x^2+bry^2+sy+bdy)
<0
\end{align*}
in the simply connected region $\mathbb{R}^2_{+}$
due to conditions \eqref{H1} and \eqref{H3}.
Thus, by Dulac's criterion \cite {DC}, there is no periodic solution
in $\mathbb{R}^2_{+}$, and this implies that the positive equilibrium
$(x^{*},y^{*})$ is globally asymptotically stable.
\end{proof}

\begin{lemma}[divergency criterion \cite{DC}] \label{DCLemma}
Assume that $L$ is a closed orbit with period $T$.
If the condition
\begin{equation*}
 \oint_0^{T}\operatorname{div}(x,y)dt<0\quad (>0)
\end{equation*}
holds, then $L$ is a single stable (unstable) limit cycle.
\end{lemma}

Based on Lemma~\ref{DCLemma}, we obtain the following result.

\begin{theorem} \label{Theorem3.2}
For system \eqref{1.3}, if condition \eqref{H1} and
\begin{equation}
 s \leq \max \{ b+4abr, b+4a \} \label{H4}
\end{equation}
hold, then the local and global asymptotic stability of the positive
equilibrium $(x^{*},y^{*})$ coincide.
\end{theorem}

\begin{proof}
For system \eqref{1.3}, the Jacobian matrix is
\begin{equation*}
J(x(t),y(t))
=\begin{pmatrix}
J_1 & J_2 \\
J_3 & J_4
\end{pmatrix},
\end{equation*}
where
\begin{gather*}
J_1=1-2x(t)-\frac{sy(t)}{x(t)+y(t)+a}+\frac{sx(t)y(t)}{(x(t)+y(t)+a)^2},\quad
J_2=-\frac{sx(t)(x(t)+a)}{(x(t)+y(t)+a)^2}, \\
J_3=\frac{by(t)(y(t)+a)}{(x(t)+y(t)+a)^2}, \quad
J_4=b[-d-2ry+\frac{x(t)}{x(t)+y(t)+a}-\frac{x(t)y(t)}{(x(t)+y(t)+a)^2}].
\end{gather*}
Assume that $l(t)=(x(t),y(t))$ is an arbitrary but fixed nontrivial periodic
orbit of system \eqref{1.3} with period $T>0$, then there holds
\begin{gather*}
\oint_0^{T}\frac{x'(t)}{x(t)}dt
=\oint_0^{T} \big[1-x(t)-\frac{sy(t)}{x(t)+y(t)+a} \big]dt=0,\\
\oint_0^{T}\frac{y'(t)}{y(t)}dt
=\oint_0^{T}b \big[-d-ry(t)+\frac{x(t)}{x(t)+y(t)+a}\big]dt=0.
\end{gather*}

If $s\leq b$, then
\[
\oint_0^{T}trJ(x(t),y(t))dt=\oint_0^{T}[-x(t)-bry(t)
+\frac{(s-b)x(t)y(t)}{(x(t)+y(t)+a)^2}]dt<0.
\]
If $s> b$, then
\[
\oint_0^{T}trJ(x(t),y(t))dt\leq\oint_0^{T}[-x(t)-\frac{(s-b)y(t)}{4a}
+\frac{(s-b)y(t)}{4a}]dt<0
\]
because of the condition $s\leq b+4abr$ and
$\oint_0^{T}trJ(x(t),y(t))dt\leq\oint_0^{T}[-x(t)+\frac{(s-b)x(t)}{4a}
-bry(t)]dt<0$
by the condition $s \leq b+4a$.


Consequently, by the divergency criterion, the closed orbit $l(t)$ is stable,
which yields a contradiction with the local asymptotic stability of the
positive equilibrium.
So, if $(x^{*},y^{*})$ is locally asymptotically stable, system \eqref{1.3}
has no nontrivial periodic orbit in $\mathbb{R}^2_{+}$.
This indicates that the positive equilibrium
must be also globally asymptotically stable.
\end{proof}


Then, from Theorems \ref{Theorem2.2} and \ref{Theorem3.2},
we can directly obtain the following corollary.

\begin{corollary} \label{coro3.4}
For system \eqref{1.3}, if conditions \eqref{H1},
\eqref{H2} and \eqref{H4} hold,
then the unique positive equilibrium is globally asymptotically stable.
\end{corollary}


\begin{remark} \label{rmk3.5} \rm
Provided that conditions \eqref{H1} and \eqref{H4} hold, we can additionally
obtain that system \eqref{1.3} has the unique (stable) limit cycle in the
first quadrant if the positive equilibrium is unstable, which will be
described as Corollary~\ref{Corollary4.10}.
\end{remark}


Note that in the proof of Theorem~\ref{Theorem3.2}, we simply apply
the mean inequality to assure the negativeness of the divergency integral.
In fact, the divergency integral can be further expanded
by applying Grammer's rule and Green's theorem.
However, this requires boundary values of periodic orbits of system \eqref{1.3}.
For this, similar to \cite[Theorem 2.2]{L15}, by defining a set
$$
\Gamma := \{(x,y) \in \mathbb{R}^2_{+} : \underline{x} \leq x \leq \overline{x},
\underline{y} \leq y \leq \overline{y} \},
$$
where $\underline{x}=1-s,~\overline{x}=1$,
\[
\underline{y}=\frac{1}{2}[-\frac{d+r(a+1-s)}{r}
+\sqrt{[\frac{d+r(a+1-s)}{r}]^2+4\frac{(1-d)(1-s)-da}{r}}]
\]
and $\overline{y}=\frac{1-d(a+1)}{d+r(a+1)}$,
we have the following stronger permanent condition for providing concrete
boundary values.

\begin{theorem} \label{Theorem3.3}
Suppose that system \eqref{1.3} satisfies the condition
\begin{equation}
 ad<(1-d)(1-s)\quad \text{and}\quad s<1. \label{H5}
\end{equation}
Then for any solution $(x(t),y(t))$
of system \eqref{1.3} with the positive initial condition (i.e.
$x(0)>0$ and $y(0)>0$), there is a $T_0>0$ such that for all $t>T_0$ ,
$(x(t),y(t)) \in \Gamma$ holds.
\end{theorem}

Obviously, since $ad<(1-d)(1-s)$ implies that condition \eqref{H1} holds,
and $s<1$ implies that condition \eqref{H2} holds, by combining
Theorem~\ref{Theorem2.2} and Corollary~\ref{coro3.4},
we directly have the following corollary.

\begin{corollary} \label{coro3.7}
If condition \eqref{H5} holds,
then the unique positive equilibrium is locally asymptotically stable.
If conditions \eqref{H4} and \eqref{H5} hold,
then the unique positive equilibrium is globally asymptotically stable.
\end{corollary}


From Theorem~\ref{Theorem3.2} and Corollary~\ref{coro3.7}, we know that
if conditions \eqref{H5} holds and $s\leq b$,
the unique positive equilibrium of system \eqref{1.3} is globally
asymptotically stable.
Further, we can derive from Theorem~\ref{Theorem3.3} that if condition
\eqref{H5} holds, all possible periodic orbits must lie in $\Gamma$.
Thus, instead of \eqref{H1}, we attempt to employ condition \eqref{H5}
for the case of $s>b$ to derive the divergency integral and obtain
the following theorem.

\begin{theorem} \label{Theorem3.4}
For system \eqref{1.3}, if the condition
\begin{equation} \label{H6}
\begin{gathered}
(s-b)(1-x^{\ast})-sb=0,\\
r<\min \{ \frac{as}{(\overline{x}+\overline{y}+a)^2(x^{*}+y^{*}+a)},
\frac{s\underline{x}}{b(\overline{x}+\overline{y}+a)^2} \},
\\
\text{or} \hfill
\\
(s-b)(1-x^{\ast})-sb>0,\\
r<\min \big\{ \frac{as}{(\overline{x}+\overline{y}+a)^2(x^{*}+y^{*}+a)},
\frac{s\underline{x}}{b(\overline{x}+\overline{y}+a)^2} \big\},\\
\underline{x}(x^{*}+\underline{x}+\underline{y}+a-1)>b\overline{y}(1-d-ry^{*}),
\end{gathered}
\end{equation}
holds, then the unique positive equilibrium
$(x^{*},y^{*})$ is globally asymptotically stable, provided that
condition \eqref{H5} holds.
\end{theorem}

\begin{proof}
Assume that $l(t)=(x(t),y(t))$ is an arbitrary but fixed nontrivial
periodic orbit of system $\eqref{1.3}$ with period $T>0$.
Similar to the proof of Theorem~\ref{Theorem3.2}, it suffices to show that
under conditions \eqref{H5} and \eqref{H6},
$\oint_0^{T}trJ(x(t),y(t))dt<0$.

Since $(x(t),y(t))$ is an orbit of system \eqref{1.3} and
$$
\oint_0^{T}trJ(x(t),y(t))dt=\oint_0^{T}
\Big[-x(t)-bry(t)+\frac{(s-b)x(t)y(t)}{(x(t)+y(t)+a)^2}\Big]dt,
$$
we have
\begin{align*}
 &\oint_0^{T}trJ(x(t),y(t))dt\\
&=\oint_0^{T}\Big[-x(t)-bry(t)+(s-b)d\frac{y(t)}{x(t)+y(t)+a}\Big]dt\\
&\quad +\oint_0^{T}\Big[(s-b)\frac{y(t)}{x(t)+y(t)+a}
 (\frac{x(t)}{x(t)+y(t)+a}-d)\Big]dt\\
&=\oint_0^{T}\Big[-x(t)-bry(t)+(s-b)\frac{d}{s}(1-x(t)-\frac{x'(t)}{x(t)})\Big]dt\\
&\quad +\oint_0^{T}\Big[(s-b)\frac{y(t)}{x(t)+y(t)+a}(ry(t)
 +\frac{y'(t)}{by(t)})\Big]dt\\
&=\oint_0^{T}\Big[-x(t)-bry(t)+(s-b)\frac{d}{s}(1-x(t))\Big]dt \\
&\quad + \oint_0^{T}\Big[(s-b)\frac{ry(t)^2}{x(t)+y(t)+a}+\frac{s-b}{b}
 \frac{y'(t)}{x(t)+y(t)+a}\Big]dt\\
&=\oint_0^{T}\Big[-x(t)-bry(t)+(s-b)\frac{1}{s}(1-x(t))(d+ry(t))\Big]dt\\
&\quad +\oint_0^{T}\Big[-(s-b)\frac{ry(t)}{s}\frac{x'(t)}{x(t)}
 +\frac{s-b}{b}\frac{y'(t)}{x(t)+y(t)+a}\Big]dt.
\end{align*}
Further, since
\[
\oint_0^{T}trJ(x^{\ast},y^{\ast})dt=-\oint_0^{T}[x^{\ast}
+bry^{\ast}-\frac{1}{s}(s-b)(1-x^{\ast})(d+ry^{\ast})]dt,
\]
it is easy to show that
\begin{equation}
\begin{split}
 &\oint_0^{T}trJ(x(t),y(t))dt\\
&=\oint_0^{T}trJ(x^{\ast},y^{\ast})dt\\
&\quad +\oint_0^{T}\Big[\frac{s-b}{b}\frac{y'(t)}{x(t)+y(t)+a}
-\frac{s-b}{s}ry(t)\frac{x'(t)}{x(t)}\Big]dt \\
&\quad + \oint_0^{T}\big[\left(-(x(t)-x^{\ast})-br(y(t)-y^{\ast})\right)\big]dt \\
&\quad + \oint_0^{T}\Big[\frac{s-b}{s}\left((d+ry(t))(1-x(t))
 -(d+ry^{\ast})(1-x^{\ast})\right)\Big]dt\\
&=\oint_0^{T}trJ(x^{\ast},y^{\ast})dt \\
&\quad +\oint_0^{T}\Big[\frac{s-b}{b}\frac{y'(t)}{x(t)+y(t)+a}
 -\frac{s-b}{s}ry(t)\frac{x'(t)}{x(t)}\Big]dt \\
&\quad + \oint_0^{T}\Big[\left(-1-\frac{1}{s}(s-b)(d+ry(t))\right)
 (x(t)-x^{\ast})\Big]dt \\
&\quad + \oint_0^{T}\Big[r\left(\frac{1}{s}(s-b)(1-x^{\ast})-b\right)
(y(t)-y^{\ast})\Big]dt.
\end{split} \label{3.1}
\end{equation}
Clearly,
\begin{align*}
&\oint_0^{T}\Big[\frac{s-b}{b}\frac{y'(t)}{x(t)+y(t)+a}
-\frac{s-b}{s}ry(t)\frac{x'(t)}{x(t)}\Big]dt \\
&=\frac{s-b}{bs} \Big(\oint_{l}-\frac{bry}{x}dx+\frac{s}{x+y+a}dy \Big).
\end{align*}
Moreover, since $(x(t),y(t))$ satisfies system \eqref{1.3} and
 $(x^{\ast},y^{\ast})$ satisfies system
\eqref{2.1}, we can obtain that
\begin{align*}
\frac{x'(t)}{x(t)}
&= 1-x(t)-\frac{sy(t)}{x(t)+y(t)+a}\\
&= x^{\ast}+\frac{sy^{\ast}}{x^{\ast}+y^{\ast}+a}-x(t)-\frac{sy(t)}{x(t)+y(t)+a}\\
&= \Big[-1+\frac{sy^{\ast}}{(x^{\ast}+y^{\ast}+a)(x(t)+y(t)+a)}
 \Big](x(t)-x^{\ast})\\
&\quad -\frac{(x^{\ast}+a)s}{(x^{\ast}+y^{\ast}+a)(x(t)+y(t)+a)}(y(t)-y^{\ast}),
\end{align*}
and
\begin{align*}
\frac{y'(t)}{by(t)}
&= -d-ry(t)+\frac{x(t)}{x(t)+y(t)+a}\\
&= ry^{\ast}-\frac{x^{\ast}}{x^{\ast}+y^{\ast}+a}-ry(t)+\frac{x(t)}{x(t)+y(t)+a}\\
&= \frac{y^{\ast}+a}{(x^{\ast}+y^{\ast}+a)(x(t)+y(t)+a)}(x(t)-x^{\ast})\\
&\quad -\Big[r+\frac{x^{\ast}}{(x^{\ast}+y^{\ast}+a)(x(t)+y(t)+a)}\Big]
(y(t)-y^{\ast}).
\end{align*}
Thus, by Crammer's Rule, we have
\begin{equation}
\begin{gathered}
x-x^{\ast}
=\frac{[\frac{x^{\ast}}{x^{\ast}+y^{\ast}+a}+r(x(t)+y(t)+a)]
\frac{x'(t)}{x(t)}-[\frac{(x^{\ast}+a)s}{x^{\ast}+y^{\ast}+a}]\frac{y'(t)}{by(t)}}
 {\frac{rsy^{\ast}-x^{\ast}}{x^{\ast}+y^{\ast}+a}-r(x(t)+y(t)+a)
-\frac{as}{(x(t)+y(t)+a)(x^{\ast}+y^{\ast}+a)}},
\\
y-y^{\ast}
=\frac{[\frac{y^{\ast}+a}{x^{\ast}+y^{\ast}+a}]\frac{x'(t)}{x(t)}
-[\frac{sy^{\ast}}{x^{\ast}+y^{\ast}+a}-(x(t)+y(t)+a)]\frac{y'(t)}{by(t)}}
 {\frac{rsy^{\ast}-x^{\ast}}{x^{\ast}+y^{\ast}+a}-r(x(t)+y(t)+a)
-\frac{as}{(x(t)+y(t)+a)(x^{\ast}+y^{\ast}+a)}}.
\end{gathered}
\label{3.2}
\end{equation}

Let $D$ denote the region encloused by the periodic orbit $l(t)$ and
\begin{gather*}
M_1(x,y)=\frac{-1-\frac{1}{s}(s-b)(d+ry(t))}{-r(x(t)+y(t)+a)
+\frac{rsy^{\ast}-x^{\ast}}{x^{\ast}+y^{\ast}+a}
 -\frac{as}{(x^{\ast}+y^{\ast}+a)(x(t)+y(t)+a)}},
\\
M_2(x,y)=\frac{r[\frac{1}{s}(s-b)(1-x^{\ast})-b]}{-r(x(t)+y(t)+a)
+\frac{rsy^{\ast}-x^{\ast}}{x^{\ast}+y^{\ast}+a}
-\frac{as}{(x^{\ast}+y^{\ast}+a)(x(t)+y(t)+a)}},
\\
M_3(x,y)=\frac{-r+\frac{as}{(x(t)+y(t)+a)^2(x^{\ast}+y^{\ast}+a)}}
{bxy\big(-r(x(t)+y(t)+a)
 +\frac{rsy^{\ast}-x^{\ast}}{x^{\ast}+y^{\ast}+a}
-\frac{as}{(x^{\ast}+y^{\ast}+a)(x(t)+y(t)+a)}\big)^2}.
\end{gather*}
Then, from \eqref{3.1} and \eqref{3.2}, we obtain
\begin{align*}
 &\oint_0^{T}\Big[-1-\frac{1}{s}(s-b)(d+ry(t))\Big](x(t)-x^{\ast})dt\\
&\quad + \oint_0^{T}r\Big[\frac{1}{s}(s-b)(1-x^{\ast})-b\Big](y(t)-y^{\ast})dt\\
&=\oint_0^{T}M_1(x,y)\Big[\Big(\frac{x^{\ast}}{x^{\ast}+y^{\ast}+a}+r(x(t)+y(t)+a)
\Big)\frac{x'(t)}{x(t)}\Big]dt\\
&\quad - \oint_0^{T}M_1(x,y)\Big[\frac{(x^{\ast}+a)s}{x^{\ast}+y^{\ast}+a}
 \frac{y'(t)}{by(t)}\Big]dt \\
&\quad -\oint_0^{T}M_2(x,y)\Big[\Big(\frac{sy^{\ast}}{x^{\ast}+y^{\ast}+a}
 -(x(t)+y(t)+a)\Big)\frac{y'(t)}{by(t)}\Big]dt \\
&\quad +\oint_0^{T}M_2(x,y)\Big[\frac{y^{\ast}+a}{x^{\ast}
 +y^{\ast}+a}\frac{x'(t)}{x(t)}\Big]dt\\
&=\oint^{}_{l}\Big[M_1(x,y)\Big(\frac{x^{\ast}}{x^{\ast}+y^{\ast}+a}
 +r(x+y+a)\Big)\Big]\frac{dx}{x} \\
&\quad + \oint^{}_{l}\Big[M_2(x,y)\frac{y^{\ast}+a}{x^{\ast}+y^{\ast}+a}\Big]
 \frac{dx}{x} \\
&\quad +\oint^{}_{l}\Big[M_2(x,y)\Big((x+y+a)-\frac{sy^{\ast}}{x^{\ast}
+y^{\ast}+a}\Big)\Big]\frac{dy}{by} \\
&\quad - \oint^{}_{l}\Big[M_1(x,y)\frac{(x^{\ast}+a)s}{x^{\ast}+y^{\ast}+a}\Big]
\frac{dy}{by}.
\end{align*}
In the following, we denote that
\[
 K_1=-1-\frac{1}{s}(s-b)(d+ry),\quad
 K_2=r(\frac{1}{s}(s-b)(1-x^{\ast})-b).
\]
Then by applying Green's theorem, we  deduce that
\begin{align*}
 &\oint_0^{T}trJ(x(t),y(t))dt\\
&=\oint_0^{T}trJ(x^{\ast},y^{\ast})dt+\iint_{D}\frac{s-b}{bs}
(-\frac{s}{(x+y+a)^2}+\frac{br}{x})\,dx\,dy\\
&\quad +\iint_{D}K_1 M_3(x,y)\Big[\frac{sx(x^{\ast}+a)}{x^{\ast}+y^{\ast}+a}
 \Big]\,dx\,dy \\
&\quad +\iint_{D}K_1 M_3(x,y)\Big[by(\frac{x^{\ast}}{x^{\ast}
 +y^{\ast}+a}+r(x+y+a))\Big]\,dx\,dy \\
&\quad +\iint_{D}K_2 M_3(x,y)\Big[\frac{by(y^{\ast}+a)}{x^{\ast}+y^{\ast}+a}
\Big]\,dx\,dy \\
&\quad +\iint_{D}K_2 M_3(x,y)\Big[x(\frac{sy^{\ast}}{x^{\ast}+y^{\ast}+a}-(x+y+a))
 \Big]\,dx\,dy \\
&\quad +\iint_{D}\frac{\frac{r}{by}[\frac{1}{s}(s-b)(1-x^{\ast})-b]}
 {-r(x+y+a)+\frac{rsy^{\ast}-x^{\ast}}{x^{\ast}+y^{\ast}+a}
 -\frac{as}{(x^{\ast}+y^{\ast}+a)(x+y+a)}}\,dx\,dy \\
&\quad + \iint_{D}\frac{\frac{r}{sx}(s-b)[\frac{x^{\ast}}{x^{\ast}
 +y^{\ast}+a}+r(x+y+a)]}
 {-r(x+y+a)+\frac{rsy^{\ast}-x^{\ast}}{x^{\ast}+y^{\ast}+a}
 -\frac{as}{(x^{\ast}+y^{\ast}+a)(x+y+a)}}\,dx\,dy \\
&\quad +\iint_{D}\frac{\frac{r}{x}[1+\frac{1}{s}(s-b)(d+ry)]}
 {-r(x+y+a)+\frac{rsy^{\ast}-x^{\ast}}{x^{\ast}+y^{\ast}+a}
 -\frac{as}{(x^{\ast}+y^{\ast}+a)(x+y+a)}}\,dx\,dy.
\end{align*} %\label{3.3}

It is easy to see that
\begin{itemize}
\item $\oint_0^{T}trJ(x^{\ast},y^{\ast})dt<0$ since $(x^{*},y^{*})$ 
is locally asymptotically stable under
 conditions \eqref{H5}.

\item $\iint_{D}(-\frac{s}{(x+y+a)^2}+\frac{br}{x})\,dx\,dy<0$ because of
  the condition $r< \frac{s\underline{x}}{b(\overline{x}+\overline{y}+a)^2}$.

\item $M_3(x,y)>0$ because of $r<\frac{as}{(\overline{x}+\overline{y}+a)^2(x^{*}+y^{*}
+a)}$.

\item the denominator in $M_1$ is negative since condition \eqref{H5}
implies that $x^{*}-rsy^{*}>1+sd-s>0$.
\end{itemize}
Thus, if $(s-b)(1-x^{\ast})-sb=0$, then $s>b$ and we have 
$\oint_0^{T}trJ(x(t),y(t))dt<0$.
And if $(s-b)(1-x^{\ast})-sb>0$, then $s>b$ and we have
\[
\frac{by(y^{\ast}+a)}{x^{\ast}+y^{\ast}+a}
+x\Big(\frac{sy^{\ast}}{x^{\ast}+y^{\ast}+a}-(x+y+a)\Big)<0
\]
by the inequality $\underline{x}(x^{*}+\underline{x}
+\underline{y}+a-1)>b\overline{y}(1-d-ry^{*})$,
which also implies that $\oint_0^{T}trJ(x(t),y(t))dt<0$.

Consequently, under conditions \eqref{H5} and \eqref{H6},
the positive equilibrium $(x^{*},y^{*})$
is globally asymptotically stable.
\end{proof}


\begin{remark} \label{rmk3.9} \rm
From the proof of Theorem~\ref{Theorem3.4}, we can find that
when $r=0$, $\oint_0^{T}trJ(x(t),y(t))dt<0$ always holds.
Thus, for system \eqref{1.3}
with $r=0$, the local and global asymptotic stability
of the positive equilibrium coincide, which was also proved in \cite {first}.
\end{remark}

\begin{remark}  \label{rmk3.10} \rm
Similarly, it is also interesting to find conditions to ensure that
$\oint_0^{T}trJ(x(t),y(t))dt>0$; that is,
\begin{equation}
\oint_0^{T}\Big[-x(t)-bry(t)+\frac{(s-b)x(t)y(t)}{(x(t)+y(t)+a)^2}\Big]dt>0,
\label{3.4}
\end{equation}
such that the local and global asymptotic stability of the positive
equilibrium $(x^{*},y^{*})$ coincide.
\end{remark}

In Theorem~\ref{Theorem3.4}, we  obtained a sufficient condition
for global stability of the positive equilibrium based on the
permanence condition \eqref{H5} which is also the local asymptotic 
stability condition.
But, this condition contains $x^{*}$ and $y^{*}$
which we need to numerically solve the positive equilibrium first.
Intuitively, condition \eqref{H6} can be simplified
by replacing $x^{*}$ and $y^{*}$ with boundary values defined in $\Gamma$.
However, for the case of $s>b$, the sub-condition $(s-b)(1-x^{*})-sb>0$
in condition \eqref{H6} can not be simplified further since $\overline{x}=1$.

\section{Hopf bifurcations with respect to one parameter $r$ and all 
five parameters}
\label{sec:3}

Compared to the systems studied in \cite{cc,first,second}, system \eqref{1.3} 
involves the density dependence of the predators and thus has the extra term 
$ry^2$ with $r>0$. In this section, we will first focus on the Hopf bifurcation 
with respect to the parameter $r$, and then extend our discussions to all 
other parameters.

From \cite{cc}, we know that system \eqref{1.3} with $r=0$ has a unique 
positive equilibrium if and only if condition \eqref{H1} holds.
Thus, system \eqref{1.3} with $r\geq 0$ has a unique positive equilibrium
if and only if condition \eqref{H1} holds.
Following \cite{zhangzhifen}, one can see that this positive equilibrium 
$(x^{*}, y^{*})$ must smoothly depend on the parameters $a>0$, $d>0$, $s>0$,
$b>0$ and $r\geq 0$.

For the Hopf bifurcation, we perform two different choices of parameters 
to analyze the change of $\mathit{Sign}(a_1)$ by following the ideas in \cite{wi}.

\begin{lemma} \label{Lemma4.1}
For system \eqref{1.3} with $r=0$,
if \eqref{H1} holds and $s>\max\{b, \frac{bd}{1+d}+\frac{1}{1-d^2} \}$,
then we have $a_1<0$ if and only if 
\[
a<\frac{1}{s}\frac{s-b}{s+(s-b)d} \Big[\frac{(s-b)d}{s+(s-b)d}+s-1-sd \Big].
\]
\end{lemma}

\begin{proof}
Let $\rho=\frac{(s-b)d}{s+(s-b)d}$ and
$\hat{a}=\frac{\rho}{sd}(\rho+s-1-sd)$.
Then we see $\hat{a}>0$ if and only if $s>\frac{bd}{1+d}+\frac{1}{1-d^2}$.
As $r=0$ in system \eqref{2.1}, it gives $ads=x^{*}(x^{*}+s-1-sd)$.
Namely, when $a=\hat{a}$, it has $x^{*}=\rho$.

As $r=0$ in system \eqref{1.3}, if $s>b$, we have
$a_1=[ 1+(1-\frac{b}{s})](x^{*}-\rho)$.
Because $x^{*}=\frac{1}{2}[1+sd-s+\sqrt{(1+sd-s)^2+4ads}]$,
$x^{*}$ increases as $a$ increases.
So, we see that $a_1<0$ if and only if $a<\hat{a}$.
\end{proof}


\begin{lemma} \label{Lemma4.2}
Under condition \eqref{H1},
if $s>S_1$ and $b+\frac{b}{s}-1\geq 0$, then $\frac{\partial a_1}{\partial r}>0$.
\end{lemma}

\begin{proof}
According to the equation $x^{*}=rsy^{*}+as\frac{d+ry^{*}}{x^{*}}+1+sd-s$,
$ry^{*}$ increases as $x^{*}$ increases.
Moreover, with the condition $b+\frac{b}{s}-1\geq 0$,
$a_1$ increases as $x^{*}$ increases since
\[
a_1=[1+d(1-\frac{b}{s})]x^{*}+( b+\frac{b}{s}-1) ry^{*}
 +(1-\frac{b}{s}) rx^{*}y^{*}-d ( 1-\frac{b}{s} ).
\]
Hence, it suffices to prove that $x^{*}$ increases as $r$ increases.

Since $(x^{*}, y^{*})$ satisfies system \eqref{2.1}, under the condition 
$s>S_1$, $(x^{*}, y^{*})$ can be determined by the two hyperbolic equations,
whose corresponding locations can be roughly shown in Figure \ref{fig1} (left) 
and Figure \ref{fig1} (middle) by the classifications provided in \cite{L14}.
Notice that the hyperbolic curves in Figure \ref{fig1} (left) do not depend 
on the parameter $r$ and
$(\frac{ad}{1-d},0)$ lies on the right hyperbolic curve in Figure \ref{fig1} 
(middle).
Further, for the second equation in system \eqref{2.1},
$\frac{dy}{dx}=\frac{1-d-ry}{r(x+2y+a)+d}>0$ holds for the right curve in 
the first quadrant, especially for $r_1<r_2$, 
$(x,y)=(\frac{ad}{1-d},0)$, and 
$\frac{1-d-r_1y}{r_1(x+2y+a)+d}>\frac{1-d-r_2y}{r_2(x+2y+a)+d}$.
Moreover, in the first quadrant, when $r_1<r_2$, the curve corresponding to $r_1$
lies above the curve corresponding to $r_2$, which is shown in 
Figure \ref{fig1} (right).
Otherwise, the curve corresponding to $r_1$ must intersect with the 
curve corresponding to $r_2$
at a point $(x_0,y_0)$ in the first quadrant. However, since
$\frac{1-d-r_1y_0}{r_1(x_0+2y_0+a)+d}<\frac{1-d-r_2y_0}{r_2(x_0+2y_0+a)+d}$,
it arrives at a contradiction.
As shown in Figure \ref{fig1} (left) and Figure \ref{fig1} (right), 
as $r$ increases, we see that $x^*$ increases. Thus, we complete the proof.
\end{proof}


\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.32\textwidth]{fig1a} % curve1.eps
 \includegraphics[width=0.32\textwidth]{fig1b} % curve2.eps
 \includegraphics[width=0.32\textwidth]{fig1c} % comparison.eps
\end{center}
 \caption{Left: Curves for the hyperbolic equation $sy=(1-x)(x+y+a)$.
 Middle: Curves for the hyperbolic equation $x=(d+ry)(x+y+a)$.
 Right: Curves for $x=(d+ry)(x+y+a)$ with $r_2>r_1>r_0=0$.}
\label{fig1}
 \end{figure}

Based on the above two lemmas, let us first focus on the Hopf bifurcation
with respect to the parameter $r$, regarding other parameters as fixed constants.
Under condition \eqref{H1},
if $s>S_1$ and $b+\frac{b}{s}-1\geq 0$, then $\frac{da_1(r)}{dr}>0$ 
holds for system \eqref{1.3} with $r\geq 0$.
Furthermore, by Lemmas~\ref{Lemma4.1} and \ref{Lemma4.2},
we have the following result:

\begin{lemma} \label{Lemma4.3}
For system \eqref{1.3}, under the same conditions as shown in
Lemmas~\ref{Lemma4.1} and \ref{Lemma4.2}, there exists a unique $r^{*}$
such that $a_1(r^*)=0$. Moreover, $a_1(r)<0$ if $r<r^{*}$ and $a_1(r)>0$ 
if $r>r^{*}$.
\end{lemma}

\begin{proof}
Let $(x_0,y_0)$ be the point in Figure \ref{fig1} (left) 
with $x_0=\frac{d(1-\frac{b}{s})}{1+d(1-\frac{b}{s})}$.
Under the conditions shown in Lemmas \ref{Lemma4.1} and \ref{Lemma4.2},
we can obtain that for the sufficiently large $r$, $x^*\geq x_0$ holds,
and hence $a_1(r)=[1+d(1-\frac{b}{s})]x^{*}+(b+\frac{b}{s}-1)ry^{*}
 +(1-\frac{b}{s})rx^{*}y^{*}-d(1-\frac{b}{s})>0$ holds.
Otherwise, for all $r\geq 0$ and $x^*< x_0$,
due to the hyperbolic curves in Figure \ref{fig1} (left) and 
Figure \ref{fig1} (middle), $y^*>y_0$ as $r \to +\infty$ and
so $ry^*\to +\infty$ as $r \to +\infty$. 
This contradicts the fact that $ry^*<1-d$.

From Lemmas~\ref{Lemma4.1} and~\ref{Lemma4.2}, there exists one unique $r^{*}$
such that $a_1(r^{*})=0$ and $a_1(r)<0$ if and only if $r<r^{*}$.
\end{proof}

From Lemma~\ref{Lemma4.3}, $a_1(r^{*})=0$ under the conditions
shown in Lemmas~\ref{Lemma4.1} and \ref{Lemma4.2}.
Moreover, recalling from Section 2 that when $s\leq\frac{1}{1-d}$, 
we have $a_2(r)>0$, and thus $a_2(r^*)>0$.
Since $a_1$ and $a_2$ smoothly depend on the parameters 
$a>0$, $b>0$, $d>0$, $s>0$ and $r\geq 0$,
there exists a neighborhood $W$ of $r^*$ such that for all $r\in W$, 
there holds $a^2_1(r)<4a_2(r)$.
This implies that the characteristic equation \eqref{2.4} has conjugate 
complex roots for all $r\in W$, denoted as
$\lambda := \alpha(r)\pm i\omega(r)
=\frac{-a_1}{2}\pm \frac{\sqrt{4a_2-a_1^2}}{2}i$.


Notice that $\alpha(r^{*})=0$ and $\omega(r^{*})>0$.
Hence the characteristic equation \eqref{2.4} has
one pair of conjugated pure imaginary roots $\lambda=\pm i\omega(r^{*}) $.
By the center manifold theorem \cite{wi}, the orbit structure
near $(x, y, r)=(x^{*},y^{*},r^{*})$ can be determined by the vector
field \eqref{2.2} restricted to the center manifold, which has the following form
\[
 \begin{pmatrix}
 \frac{dx}{dt} \\
 \frac{dy}{dt}
 \end{pmatrix}
= \begin{pmatrix}
 \alpha(r) & -\omega(r) \\
 \omega(r) & \alpha(r)
 \end{pmatrix}
 \begin{pmatrix}
 x \\ y
 \end{pmatrix}+
 \begin{pmatrix}
 f^1(x,y,r) \\
 f^2(x,y,r)
 \end{pmatrix}
\]
where $f^1$ and $f^2$ are nonlinear in $x$ and $y$.
Then, by the Hopf bifurcation theory \cite{wi}, we can directly obtain 
the following result on Hopf bifurcation for system \eqref{1.3}.

\begin{theorem} \label{Theorem4.1}
For system \eqref{1.3}, if the following conditions
$ad < 1-d$, $s\leq\frac{1}{1-d}$, $b+\frac{b}{s}-1\geq 0$,
$s>\max\{b, 2a, 1+a, \frac{bd}{1+d}+\frac{1}{1-d^2}\}$ and
\[
a<\frac{1}{s}\frac{s-b}{s+(s-b)d}\Big[\frac{(s-b)d}{s+(s-b)d}+s-1-sd\Big]
\]
hold, then the positive equilibrium
$(x^{*},y^{*})$ is an unstable focus for $0<r^{*}-r\ll 1$ and
an asymptotically stable focus for $0<r-r^{*}\ll 1$.
Moreover, when $u(r^{*})>0$, there exists a neighborhood $U$ of the 
positive equilibrium $(x^{*},y^{*})$
such that the system has a unique unstable periodic orbit in $U$ 
for $0<r-r^{*}\ll 1$ (and hence a stable periodic orbit exists
outside this unstable periodic orbit). When $u(r^{*})<0$, there exists a 
neighborhood $U$ of the positive equilibrium $(x^{*},y^{*})$
such that the system has a unique stable periodic orbit in $U$ for $0<r^{*}-r\ll 1$,
where $u(r^{*})$ is given in \cite{wi} as:
 \begin{align*}
  u(r^{*})
&=\frac{1}{16}[f_{xxx}^1+f_{xyy}^1+f_{xxy}^2+f_{yyy}^2]
+ \frac{1}{16\omega(r^{*})}\Big[f_{xy}^1(f_{xx}^1+f_{yy}^1)\\
&\quad -f_{xy}^2(f_{xx}^2+f_{yy}^2)-f_{xx}^1f_{xx}^2+f_{yy}^1f_{yy}^2\Big].
\end{align*}
\end{theorem}

In Theorem~\ref{Theorem4.1}, we have discussed the existence of limit cycles 
for sufficiently small $|r-r^*|$.
In fact, combining Lemmas~\ref{Lemma4.3} and \ref{DCLemma} with the proof 
of Theorem~\ref{Theorem3.2}, it is straightforward to obtain the following 
corollary.

\begin{corollary} \label{Corollary4.10}
For system \eqref{1.3}, if the conditions shown in Theorem~\ref{Theorem4.1} hold,
then there is at least one stable limit cycle when $r<r^{*}$.
Moreover, the limit cycle is unique if condition \eqref{H4} is also satisfied.
\end{corollary}

Additionally, we have the following stability results on the positive equilibrium.

\begin{corollary} \label{Corollary4.100}
For system \eqref{1.3}, if the conditions shown in Theorem~\ref{Theorem4.1} hold,
then the positive equilibrium is unstable if $r<r^{*}$
and locally asymptotically stable if $r>r^{*}$. Moreover,
the positive equilibrium is globally asymptotically stable when $r>r^{*}$ 
and condition \eqref{H4} holds.
\end{corollary}

Up to now, we have discussed the Hopf bifurcation only
with respect to the parameter $r$.
It is also interesting to discuss the Hopf bifurcation with respect to all parameters.
In the rest of this section, we would like to set a special geometrical structure
such that the analysis process can be simplified by applying the classical Hopf
bifurcation theory on systems with one single parameter.

For this, let $\mu=\frac{(b-s)x^{*}y^{*}}{(x^{*}+y^{*}+a)^2}+x^{*}+bry^{*}$ 
be an auxiliary map, which can be also regarded as an auxiliary parameter and 
will be used to analyze the Hopf bifurcation.
Intuitively, by this parameter $\mu$, we restrict the five parameters of 
system \eqref{1.3} to a special geometrical structure such that the 
five-dimensional parameters can be simply mapped to one-dimensional parameter.
Thus, it enables us to consider the Hopf bifurcation along with
this special geometrical structure,
where the five parameters to some extent affect together as one parameter $\mu$.
Note that from Lemmas~\ref{Lemma4.1} and \ref{Lemma4.3},
there certainly exist conditions on the five parameters $a,b,d,r,s$
to assure the possibilities of $\mu >0,~\mu<0$ and $\mu=0$, respectively.


Moreover, since $a_1$ and $a_2$ smoothly depend on the five parameters,
$a_2>0$ if $s\leq\frac{1}{1-d}$ and
$\mu$ can be regarded as a continuous map of the five parameters.
For the parameters $a,b,d,r,s$ satisfying $\mu(a,b,d,r,s)=0$,
there exists a neighborhood $K\subsetneq \mathbb{R}^5$ of $(a,b,d,r,s)$
(and thus a neighborhood $\Omega\subsetneq \mathbb{R}$ of $\mu=0$)
such that for all $(a,b,d,r,s)\in K$ (and thus $\mu \in \Omega$), 
$a^2_1(a,b,d,r,s)<4a_2(a,b,d,r,s)$ holds.
This implies that the characteristic equation \eqref{2.4} has conjugate 
complex roots in $K$, denoted as $\lambda := \beta(\mu)\pm i\varphi(\mu)
=\frac{-a_1}{2}\pm \frac{\sqrt{4a_2-a_1^2}}{2}i$.

Clearly, when $\mu=0$, $\beta(0)=0$, $\beta'(0)=-\frac{1}{2}< 0$
and $\varphi(0)>0$, the characteristic equation \eqref{2.4} has
one pair of conjugated pure imaginary roots $\lambda=\pm i\varphi(0)$.
By the center manifold theorem \cite{wi}, we can see that the orbit structure
near $(x^*, y^*, a,b,d,r,s)$ with $\mu(a,b,d,r,s)=0$ (i.e. near $(x^*,y^*,\mu)$ 
with $\mu=0$) is determined by the vector field \eqref{2.2} restricted to 
the center manifold, which has the  form
\[
 \begin{pmatrix}
 \frac{dx}{dt} \\
 \frac{dy}{dt}
 \end{pmatrix}
= \begin{pmatrix}
 \beta(\mu) & -\varphi(\mu) \\
 \varphi(\mu) & \beta(\mu)
 \end{pmatrix}
 \begin{pmatrix}
 x \\ y
 \end{pmatrix}+
 \begin{pmatrix}
 f^1(x,y,\mu) \\
 f^2(x,y,\mu)
 \end{pmatrix}
\]
where $f^1$ and $f^2$ are nonlinear in $x$ and $y$.

Since $\beta'(0)=-\frac{1}{2}<0$, by the classical Hopf bifurcation
theory on dynamical systems with one single parameter \cite{wi},
we can obtain the following theorem on the Hopf bifurcation with respect
to all the five parameters along with the special geometrical structure
$\mu=\frac{(b-s)x^{*}y^{*}}{(x^{*}+y^{*}+a)^2}+x^{*}+bry^{*}$.

\begin{theorem} \label{Theorem4.2}
For system \eqref{1.3}, if $s\leq\frac{1}{1-d}$, then
the positive equilibrium $(x^{*},y^{*})$ is an unstable focus for $-1 \ll\mu <0$,
or an locally asymptotically stable focus for $0<\mu \ll1$.
Moreover, when $\psi(0)>0$, there exists a neighborhood $\bar{U}$ of the positive
equilibrium $(x^{*},y^{*})$ such that the system has a unique unstable periodic
orbit in $\bar{U}$ for $0<\mu \ll1$ (and hence a stable periodic orbit exists
outside this unstable periodic orbit).
When $\psi(0)<0$, there exists a neighborhood $\bar{U}$ of the positive
equilibrium $(x^{*},y^{*})$ such that the system has a unique
stable periodic orbit in $\bar{U}$ for $-1 \ll\mu <0$,
where the coefficient $\psi(0)$ is given as in \cite{wi}:
\begin{align*}
\psi(0)&=\frac{1}{16}[f_{xxx}^1+f_{xyy}^1+f_{xxy}^2+f_{yyy}^2]
+\frac{1}{16\varphi(0)}\Big[f_{xy}^1(f_{xx}^1+f_{yy}^1) \\
&\quad -f_{xy}^2(f_{xx}^2+f_{yy}^2)-f_{xx}^1f_{xx}^2+f_{yy}^1f_{yy}^2\Big].
 \end{align*}
\end{theorem}

\section{Numerical simulations}

\label{sec:4}

In this section, we illustrate some numerical examples.

\begin{example} \label{Example5.1} \rm
Let $a=3$, $b=5$, $s=1$, $d=0.2$ and $r=1$,
then system \eqref{1.3} becomes
\begin{equation}
\begin{gathered}
x'(t)=x(t)(1-x(t))-\frac{x(t)y(t)}{x(t)+y(t)+3},\\
y'(t)=5y(t) \big(-0.2-y(t)+\frac{x(t)y(t)}{x(t)+y(t)+3} \big).
\end{gathered} \label{5.1}
\end{equation}
Since $0.6=ad<1-d=0.8$ and $1=s<2a=6$, according to Theorem~\ref{Theorem2.2},
the positive equilibrium $(x^{*},y^{*})\approx(0.9888,0.0451)$ of
system \eqref{5.1} is locally asymptotically stable.
In fact, the positive equilibrium is globally asymptotically stable too,
since $1=bd \geq 1$ and $1=2a-b \geq 1$, and
condition \eqref{H3} in Theorem~\ref{Theorem3.1} holds.
The global asymptotic stability can be seen from Figure \ref{fig2} (left). 
Note that in Figure \ref{fig2} (left),
all six orbits go to $(0.9888,0.0451)$
as $t$ tends to $+\infty$, starting from initial points
$(0.1,0.2)$, $(0.5,0.3)$, $(0.7,0.6)$,
$(0.9,1.0)$, $(1.2,1.05)$ and $(1.5,1.35)$,
respectively.
\end{example}


\begin{example} \label{Example5.2} \rm
Let $a=1$, $b=1$, $s=0.5$, $d=0.1$ and $r=1$,
then system \eqref{1.3} becomes
\begin{equation}
\begin{gathered}
x'(t)=x(t)(1-x(t))-\frac{0.5x(t)y(t)}{x(t)+y(t)+1},\\
y'(t)=y(t) \Big(-0.1-y(t)+\frac{x(t)y(t)}{x(t)+y(t)+1} \Big).
\end{gathered}\label{5.2}
\end{equation}
Since $0.1=ad<(1-d)(1-s)=0.45$ and $0.5=s<1$, according to Theorem~\ref{Theorem3.3},
the positive equilibrium
$(x^{*},y^{*})\approx(0.9300, 0.3144)$ of system \eqref{5.2} is
locally asymptotically stable.
In fact, by Theorem~\ref{Theorem3.2},
the positive equilibrium is globally asymptotically stable too,
since $0.5=s<b=1$.
The global asymptotic stability can be seen from Figure \ref{fig2} (right).
Note that in Figure \ref{fig2} (right), all six orbits all go to $(0.9300, 0.3144)$
as t tends to $+\infty$, starting from initial points
$(0.1,0.2)$, $(0.5,0.3)$, $(0.7,0.6)$,
$(0.9,1.0)$, $(1.2,1.05)$ and $(1.5,1.35)$,
respectively.
\end{example}


\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.48\textwidth]{fig2a} % Fig1.eps
 \includegraphics[width=0.48\textwidth]{fig2b} % Fig2.eps
\end{center}
 \caption{Left: Evolutions of six orbits for system \eqref{5.1}.
 Right: Evolutions of six orbits for system \eqref{5.2}.}
\label{fig2}
 \end{figure}


\begin{example} \label{Example5.3} \rm
Let $a=2$, $b=1/20$, $s=1/2$, $d=1/6$, and $r=1/750$,
then system \eqref{1.3} becomes
\begin{equation}
\begin{gathered}
x'(t)=x(t)(1-x(t))-\frac{\frac{1}{2}x(t)y(t)}{x(t)+y(t)+2},\\
y'(t)=\frac{1}{20}y(t) \Big(-\frac{1}{6}-\frac{1}{750}y(t)
+\frac{x(t)y(t)}{x(t)+y(t)+2} \Big). \label{5.3}
\end{gathered}
\end{equation}
The positive equilibrium is $(x^{*},y^{*})\approx (0.7969, 1.9126)$.
A straightforward calculation gives that $\frac{1}{2}=s<1$,
$\frac{1}{3}=ad<(1-d)(1-s)=\frac{5}{12}$,
$(s-b)(1-x^{\ast})-sb \approx 0.0664 > 0$,
$\frac{1}{750}=r<\frac{as}{(a+1+\frac{1}{d})^3}=\frac{1}{729}$,
$\frac{1}{750}=r<\frac{s(1-s)}{(a+1+\frac{1}{d})^2}=\frac{1}{324}$, and
$\frac{1}{20}=b<d(a-1)(1-s)=\frac{1}{12}$. This
implies that conditions \eqref{H5} and \eqref{H6} in 
Theorem~\ref{Theorem3.4} are satisfied.
Thus, $(x^{*},y^{*})$ is globally attractive, which can be seen from
 Figure \ref{fig3}.
Note that in Figure \ref{fig3}, the six orbits starting from initial points
$(0.1,1.7)$, $(0.3,3.5)$, $(0.5,0.3)$,
$(1,2.7)$, $(1.2,1.05)$, $(1.5,3)$,
respectively, go to $(0.7969, 1.9126)$
as t tends to $+\infty$.
\end{example}

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.6\textwidth]{fig3} % Fig3.eps
\end{center}
 \caption{Evolutions of six orbits for system \eqref{5.3}.}
\label{fig3}
\end{figure}

\begin{example} \label{Example5.5} 
Let $a=\frac{1}{1000}$, $b=1$, $s=\frac{4}{3}$, and $d=1/4$,
then system \eqref{1.3} becomes
\begin{equation}
\begin{gathered}
x'(t)=x(t)(1-x(t))-\frac{\frac{4}{3}x(t)y(t)}{x(t)+y(t)+\frac{1}{1000}},\\
y'(t)=y(t) \Big(-\frac{1}{4}-ry(t)+\frac{x(t)}{x(t)+y(t)+\frac{1}{1000}} \Big).
\end{gathered}\label{5.5}
\end{equation}

A direct calculation shows that $\frac{1}{4000}=ad<1-d=\frac{3}{4}$,
$\frac{4}{3}=s>\max \{b, 1+a, 2a, \frac{bd}{1+d}+\frac{1}{1-d^2} \}
=\frac{19}{15}$,
\[
\frac{1}{1000}=a<\frac{1}{s}\frac{s-b}{s+(s-b)d} 
\big(\frac{(s-b)d}{s+(s-b)d}+s-1-sd \big)=\frac{3}{289},
\]
and $\frac{3}{4}=b+\frac{b}{s}-1\geq 0$.
This indicates that the conditions in Lemmas~\ref{Lemma4.1} and \ref{Lemma4.2} hold,
and the condition $\frac{4}{3}=s\leq\frac{1}{1-d}=\frac{4}{3}$ also holds.
Hence, there is one unique $r^{*}$, which can be calculated numerically as
$r^{*}\approx 0.2250395$, and then $u(r^{*})\approx 9929481>0$.
We can choose $r=0.228>r^{*}$, and find $0.0045=a_2>0.0001>\frac{1}{4}a^2_1$.
By Theorem~\ref{Theorem4.1}, the positive equilibrium
$(x^{*},y^{*})\approx( 0.0422, 0.1101)$ is locally asymptotically
stable and there is at least one unstable closed orbit and one stable closed
orbit in the first quadrant.
Note that in Figure \ref{fig4} (left), the orbit starting from initial
points $(0.0416,0.109)$ and $(0.0413,0.1081)$ tends to a stable limit cycle 
as $t$ approaches $+\infty$;
and in Figure \ref{fig4} (middle) and 4 (right), $(x(t), y(t))$
starting from initial points $(0.0416,0.109)$ tends to a periodic form as 
$t$ approaches $+\infty$.
Note that the unstable limit cycle in
the region enclosed by the stable limit cycle is very close to the positive 
equilibrium and thus has not been shown in Figure \ref{fig4} (left).
\end{example}


\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.32\textwidth]{fig4a} % Fig5.eps
 \includegraphics[width=0.32\textwidth]{fig4b} % Fig5a.eps
 \includegraphics[width=0.32\textwidth]{fig4c} % Fig5b.eps
\end{center}
\caption{Phase diagram (left) for system \eqref{5.5} with evolution of 
$x(t)$ (middle) and $y(t)$ (right) as $t\to +\infty$.}
\label{fig4}
\end{figure}

\subsection*{Conclusions and future works}

In this article, we considered the stability of the unique positive equilibrium
and Hopf bifurcation with respect to parameters in a density-dependent
 predator-prey system with the
Beddington-DeAngelis functional response.
We started with the existence and uniqueness of the positive equilibrium, 
which can not be a saddle, and provided first a weaker sufficient condition 
by the Lyapunov function method and then a concrete condition only depending 
on parameters by the Routh-Hurwitz criterion for local stability.

Moreover, we presented several sufficient conditions for global stability 
of the positive equilibrium by two classical criteria.
That is, by Dulac's criterion, we directly obtained \eqref{H1} and \eqref{H3} 
for global stability. By the divergency criterion, we established 
\eqref{H1}, \eqref{H2} and \eqref{H4} as the sufficient conditions of global 
attractiveness. By Grammer's rule and Green's theorem,
we derived the divergency integral and further obtained \eqref{H5} and \eqref{H6}
as the sufficient conditions of global attractiveness.

Afterwards, we analyzed the Hopf bifurcation with respect to the
parameter $r$ by exploring the monotonicity of $a_1(r)$ and successively 
obtaining a unique $r^{*}$ such that $a_1(r^*)=0$.
Furthermore, we introduced an auxiliary map 
$\mu=\frac{(b-s)x^{*}y^{*}}{(x^{*}+y^{*}+a)^2}+x^{*}+bry^{*}$
to restrict the five parameters to a special one-dimensional geometrical 
structure and analyzed the Hopf bifurcation with respect to all five 
parameters along with this geometrical restriction.

Note that Hwang verified that for system \eqref{1.1}, the local stability 
and global stability of the positive equilibrium coincide in \cite{first}. However,
from analysis in Section 4, we find that for system \eqref{1.3},
the coincidence between local stability and global stability does not hold 
because of the occurrence of the parameter $r$.
Consequently, the analysis results show that the predator density dependence 
rate $r$ has a significant effort on the global dynamics of system \eqref{1.3}.

Finally, some numerical simulations have been performed to
illustrate our analytical results.

Notice that from the analysis in Section 4,
system \eqref{1.3} has limit cycles under certain conditions.
However, the problem on the number of limit cycles is not involved at this stage.
So, it is interesting to further explore the number of the limit cycles
with their location estimate in our future work,
as well as the necessary and sufficient conditions for the uniqueness 
of limit cycles.

\subsection*{Acknowledgments}
 This work is supported by National Science Foundation of China under 11290141, 
11371047 and 11422111.



\begin{thebibliography}{00}

\bibitem{Bio1} P. A. Abrams, L. R. Ginzburg;
\emph{The nature of predation: prey-dependent, ratio-dependent or neither?},
Trends Ecol. Evol., 15, 337--341 (2000).

\bibitem{bs} D. D. Bainov, P. S. Simeonov;
\emph{Systems with impulse effect: stability
theory and applications}. Ellis Horwood Limited, Chichester (1989).

\bibitem{bss} D. D. Bainov, P. S. Simeonov;
\emph{Impulsive differential equations: periodic solutions and applications}.
Longman Scientific and Technical, New York (1993),

\bibitem{b} J. R. Beddington;
\emph{Mutual interference between parasites or predators and its effect 
on searching efficiency}, J. Animal. Ecol., 44, 331-340 (1975).

\bibitem{Method1} E. Beretta, Y. Kuang;
\emph{Geometric stability switch criteria in delay differential
systems with delay dependent parameters},
SIAM J. Math. Anal., 33, 1144--1165 (2002).

\bibitem{cc} R. S. Cantrell, C. Cosner;
\emph{On the dynamics of predator-prey models with
the Beddington-DeAngelis functional response},
J. Math. Anal. Appl., 257, 206-222 (2001).

\bibitem{Bio2} C. Cosner, D. L. DeAngelis, J. S. Ault, D. B. Olson;
\emph{Effects of spatial grouping on the functional response of predators},
Theor. Pop. Biol., 56, 65--75 (1999).

\bibitem{G1} J. M. Cushing;
\emph{Periodic time-dependent predator-prey systems},
SIAM J. Appl. Math., 32, 82--95 (1977).

\bibitem{d} D. L. DeAngelis, R. A. Goldstein  R. V. O'Neil;
\emph{A model for trophic interaction}, Ecology, 56, 881-892 (1975).

\bibitem{du1} Z. J. Du, X. Chen, Z. Feng;
\emph{Multiple positive periodic solutions to a predator-prey model
with Leslie-Gower Holling-type II functional response and harvesting terms},
Discrete Contin. Dyn. Syst. Ser. S, 7(6), 1203-1214 (2014).

\bibitem{du2} Z. J. Du, Z. Feng;
\emph{Periodic solutions of a neutral impulsive predator-prey model with
Beddington-DeAngelis functional response with delays},
J. Comput. Appl. Math., 258, 87-98 (2014).

\bibitem{PB} W. Hahn;
\emph{Stability of Motion}. Springer, (1967).

\bibitem{DC} J. Hale;
\emph{Ordinary Differential Equations}. Krieger, Malabar (1980).

\bibitem{Method2} J. K. Hale, P. Waltman;
\emph{Persistence in infinite-dimensional systems},
SIAM J. Math. Anal., 20, 388--395 (1989).

\bibitem{first} T. W. Hwang;
\emph{Global analysis of the predator-prey system with
{B}eddington-{D}e{A}ngelis functional response},
J. Math. Anal. Appl., 281, 395-401 (2003).

\bibitem{second} T. W. Hwang;
\emph{Uniqueness of limit cycles of the predator-prey
system with Beddington-DeAngelis functional responsee},
J. Math. Anal. Appl., 290, 113--122 (2004).

\bibitem{pmab} P. Kratina, M. Vos, A. Bateman, B. R. Anholt;
\emph{Functional response modified by predator density}, Oecologia,
159, 425-433 (2009).

\bibitem{L14} H. Li, Z. She;
\emph{A Density-Dependent Predator-Prey Model
with {B}eddington-{D}e{A}ngelis Type}, Electron. J. Differ. Eq.,
192, 1-15 (2014).

\bibitem{L15} H. Li, Z. She;
\emph{Uniqueness of periodic solutions of a nonautonomous
density-dependent predator-prey system}, J. Math. Anal. Appl.,
442, 886-905 (2015).

\bibitem{L1502} H. Li, Z. She;
\emph{Dynamics of a nonautonomous density-dependent predator-prey model with
Beddington-DeAngelis functional response}, Int. J. Biomath., 9(4), 
Article ID 1650050, 25 pages (2016).
doi:10.1142/S1793524516500509.

\bibitem{L11} H. Li, Y. Takeuchi;
\emph{Dynamics of the density dependent predator-prey system with
{B}eddington-{D}e{A}ngelis functional response},
J. Math. Anal. Appl., 374, 644--654 (2001).

\bibitem{G2} S. Liu, E. Beretta;
\emph{A stage-structured predator-prey model of {B}eddington-{D}e{A}ngelis type},
SIAM J. Appl. Math., 66, 1101--1129 (2006).

\bibitem{L13} Z. She, H. Li;
\emph{Dynamics of a desity-dependent stage-structured predator-prey system with
Beddington-DeAngelis functional response},
J. Math. Anal. Appl., 406, 188-202 (2013).

\bibitem{Method3} H. R. Thieme;
\emph{Persistence under relaxed point-dissipativity (with application 
to an endemic model)},
SIAM J. Math. Anal., 24, 407--435 (1993).

\bibitem{G3} Z. Wang J. Wu;
\emph{Qualitative analysis for a ratio-dependent predator-prey model 
with stage-structure and diffusion},
Nonlinear Anal. RWA, 9, 2270--2287 (2008).

\bibitem{wi} S. Wiggins;
\emph{Introduction to applied nonlinear dynamical systems
and chaos}. Spring-Verlag, New York (1990).

\bibitem{zhangzhifen} Z. Zhang, C. Li, Z. Zheng,  W. Li;
\emph{Bifurcation theory of vector fields} (in Chinese).
Higher Education Press, Beijing (1997).

\bibitem{G4} J. Zhao, J. Jiang;
\emph{Permanence in nonautonomous {L}otka-{V}olterra system with predator-prey}, 
Appl.~Math.~Comput., 152, 99--109 (2004).

\bibitem{G5} H. Zhu, S. Campbell, G. Wolkowicz;
\emph{Bifurcation analysis of a predator-prey system with nonmonotonic functional
response}, SIAM J. Appl. Math., 63, 636--682 (2002).


\end{thebibliography}

\end{document}