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

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2012 (2012), No. 150, pp. 1--13.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2012 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2012/150\hfil Multiple positive periodic solutions]
{Multiple positive periodic solutions for a food-limited two-species
ratio-dependent predator-prey patch system with delay and harvesting}

\author[H. Fang \hfil EJDE-2012/150\hfilneg]
{Hui Fang}

\address{Hui Fang \newline
Department of Mathematics\\
Kunming University of Science and Technology\\
Kunming, Yunnan 650500, China}
\email{kmustfanghui@hotmail.com}

\thanks{Submitted April 18, 2012. Published August 29, 2012.}
\thanks{Supported by grants 10971085 and 11061016 from the National Natural
Science Foundation \hfill\break\indent  of China}
\subjclass[2000]{34C25, 92D25}
\keywords{Coincidence degree; periodic solution; food-limited supply; 
\hfill\break\indent
predator-prey patch system; harvesting rate}

\begin{abstract}
 By using Mawhin's coincidence degree theory, this paper establishes some
 sufficient conditions on the existence of four positive periodic
 solutions for a food-limited two-species ratio-dependent
 predator-prey patch system with delay and harvesting. Some novel
 techniques are employed to obtain the appropriate {\it a priori}
 estimates. An example is given to illustrate our results.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\allowdisplaybreaks


\section{Introduction}

Since the exploitation of biological resources and the harvest of
population species are related to the optimal management of
renewable resources (see \cite{c3}), many researchers pay much attention
to the study of population dynamics with harvesting in mathematical
bioeconomics. For example, Brauer and Soudack \cite{b1,b2} analyzed the
global behaviour of some predator-prey systems under constant rate
harvesting and/or stocking of both species. Xiao and Jennings \cite{x2}
studied the Bogdanov-Takens bifurcations in predator-prey systems
with constant rate harvesting. But all the coefficients in the
system they studied are constants. Recently, some researchers
studied the existence of multiple positive periodic solutions for
some predator-prey systems under the assumption of periodicity of
the parameters by using Mawhin's coincidence degree theory 
(see \cite{c2,f1,x1,z}). These papers only focused on predator-prey systems without
diffusion. However, due to the spatial heterogeneity and unbalanced
food resources, the migration phenomena of biological species can
often occur between heterogeneous spatial environments or patches.
On the other hand, as pointed out by Huusko and Hyvarinen in \cite{h}:
"the dynamics of exploited populations are clearly affected by
recruitment and harvesting, and the changes in harvesting induced a
tendency to generation cycling in the dynamics of a freshwater fish
population". So, it is very important to describe the potential role
of harvesting as an external factor in producing and maintaining the
periodic fluctuation of population of species. Mathematically, this
is equivalent to the investigation of harvesting induced periodic
solutions.

In this paper, we consider the following food-limited two-species
ratio-dependent predator-prey patch system with delay and
harvesting:
\begin{equation}\label{e1.1}
\begin{aligned}
x_1'(t)
&=\frac{x_1(t)}{k_1(t)+c_1(t)x_1(t)}
\Big[a_1(t)-a_{11}(t)x_1(t)-\frac{a_{13}(t)x_1(t)y(t)}{m(t)y^2(t)+x_1^2(t)}\Big]\\
&\quad + D_1(t)[x_2(t)-x_1(t)]-H_1(t),\\
x_2'(t)
&=\frac{x_2(t)}{k_2(t)+c_2(t)x_2(t)}\big[
a_2(t)-a_{22}(t)x_2(t)\big]
 +D_2(t)[x_1(t)-x_2(t)]-H_2(t), \\
y'(t)
&=y(t)\Big[-a_3(t)+\frac{a_{31}(t)x_1^2(t-\tau)}{m(t)y^2(t-\tau)+x_1^2(t-\tau)}\Big],
\end{aligned}
\end{equation}
where $x_1$ and $y$ are the population numbers of prey species $x$
and predator species $y$ in patch 1, and $x_2$ is the population
number of species $x$ in patch 2. Prey dispersion is included in the
model to allow for prey refuge from predation. So, the predator
species $y$ is confined to patch 1, while the prey species $x$ can
diffuse between two patches. For biological relevance of allowing
for prey refuge from predation, see Yakubu \cite{y}. $a_1(t)(a_2(t))$ is
the natural growth rate of prey species $x$ in patch 1 (patch 2) in
the absence of predation, $a_3(t)$ is the natural death rate of
predator species $y$ in the absence of food, $a_{13}(t)$ is the
death rate per encounter of prey species $x$ due to predation,
$a_{31}(t)$ is the efficiency of turning predated prey species $x$
into predator species $y$. $k_i(t)$ $(i=1,2)$ are the population
numbers of prey species in patch 1 (patch 2) at saturation,
respectively.

When $c_i(t)\neq 0(i=1,2),\frac{a_i(t)}{k_i(t)c_i(t)}$ $(i=1,2)$ are the rate of
replacement of mass in the population at saturation (including the
replacement of metabolic loss and of dead organisms). The effect of
delay $\tau$ refers to the dynamics of a predator being related to
the predation in the past. $m(t)$ is the half capturing saturation
coefficient. $a_{ii}(t)$ $(i=1,2)$ are the density-dependent
coefficients. $D_i(t)$ $(i=1,2)$ are diffusion coefficients of species $x$.
 $H_i(t)$ $(i=1,2)$ denote the harvesting rates.
$$
\frac{a_{13}(t)x_1(t)y(t)}{m(t)y^2(t)+x_1^2(t)},\quad
\frac{a_{31}(t)x_1^2(t-\tau)}{m(t)y^2(t-\tau)+x_1^2(t-\tau)}
$$
 are ratio-dependent Holling type III functional responses. 
When $c_i(t)\equiv 0$  $(i=1,2)$, $H_i(t)\equiv 0$  $(i=1,2)$, Dong and Ge in
\cite{d} showed that \eqref{e1.1} has at least one positive periodic solution
under the appropriate conditions. When $c_i(t)\neq 0$ $(i=1,2)$, system
\eqref{e1.1} is a food-limited population model. For other food-limited
population models, we refer to \cite{g2,g3,g4,li} and the references cited
therein.

To our knowledge, a few papers have been published on the existence of
multiple periodic solutions for \eqref{e1.1}. In this paper, we study the
existence of multiple positive periodic solutions of \eqref{e1.1} by using
Mawhin's coincidence degree. Since system \eqref{e1.1} involves the
diffusion terms and the rates of replacement, the methods used in
\cite{c2,f1,x1,z} are not available to system \eqref{e1.1}.
Our method of defining the operator $N(u,\lambda)$ facilitates the
computation of Brouwer degree $\deg_B(JQN(\cdot,0)|_{\ker L},\Omega\cap \ker L,0)$.


\section{Existence of four positive periodic solutions}

To give the proof of the main results, we first make the
following preparations (see \cite{g1}).

Let $X, Z$ be normed vector spaces, $L: \operatorname{dom}L\subset X\to Z$ 
a linear mapping, $N: X\times [0,1]\to Z$ is a
continuous mapping. The mapping $L$ will be called a Fredholm
mapping of index zero if $dim \ker L=codim \operatorname{Im}L<+\infty$ and
 $\operatorname{Im}L$ is
closed in $Z$. If $L$ is a Fredholm mapping of index zero, there
then exist continuous projectors $P:  X\to X$ and $Q:
Z\to Z$ such that $ImP=\ker L, \operatorname{Im}L=\ker Q=Im(I-Q)$. If we
define $L_P: \operatorname{dom}L\cap \ker P\to \operatorname{Im}L$ as the restriction
$L|_{\operatorname{dom}L\cap \ker P}$ of $L$ to $\operatorname{dom}L\cap \ker P$, 
then $L_P$ is invertible. We denote the inverse of
that map by $K_P$. If $\Omega$ is an open bounded subset of $X$,
the mapping $N$ will be called $L$-compact on $\bar{\Omega}\times [0,1]$ if
$QN(\bar{\Omega}\times [0,1])$ is bounded and $K_P(I-Q)N:\bar{\Omega}\times [0,1]\to X$ is
compact; i.e., continuous and such that $K_P(I-Q)N(\bar{\Omega}\times [0,1])$ 
is relatively compact.  Since $\operatorname{Im}Q$ is isomorphic to $\ker L$,  
there exists  isomorphism $J:  \operatorname{Im}Q\to \ker L$.

For convenience, we introduce Mawhin's continuation theorem  \cite[p.29]{g1}
as follows.

\begin{lemma} \label{lem2.1}
Let $L$ be a Fredholm mapping of index zero and let $N:\bar{\Omega}\times
[0,1]\to Z$ be $L$-compact on $\bar{\Omega}\times [0,1]$.
Suppose
\begin{itemize}
\item[(a)] $Lu\neq \lambda N(u,\lambda)$ for every $u\in \operatorname{dom}L\cap\partial \Omega$ and every $\lambda\in (0,1)$;

\item[(b)] $QN(u,0)\neq 0$ for every $u\in \partial\Omega\cap \ker L$;

\item[(c)] Brouwer degree $\deg_B(JQN(\cdot,0)|_{\ker L},\Omega\cap \ker
L,0)\neq 0$.
\end{itemize}
Then $Lu=N(u,1)$ has at least one solution in
$\operatorname{dom}L\cap\bar{\Omega}$.
\end{lemma}

For the sake of convenience and simplicity, we denote
 \begin{gather*}
 \bar{g}=\frac{1}{T}\int_0^Tg(t)dt,\quad
g^{l}=\min_{{t}\in[0,T]}g(t),\quad
g^u=\max_{t\in[0,T]}g(t),
 \end{gather*}
where $g$ is a nonnegative continuous $T$-periodic function.
Set
\begin{gather*} N_{1}=\max\big\{(\frac{a_{1}}{a_{11}})^u,
(\frac{a_{2}}{a_{22}})^u\big\},\quad
b_1=\frac{(\frac{a_{13}}{k_1})^u}{2\sqrt{m^l}},\quad b_2=0.
\end{gather*}
For the rest of this article, we assume the following:
\begin{itemize}
\item[(A1)] $\tau>0$, $k_i(t)$ $(i=1,2)$, $a_i(t)$ $(i=1,2,3)$,
$a_{ii}(t)$ $(i=1,2)$, $a_{13}(t)$, $a_{31}(t)$, $m(t)$,
$D_i(t)$ $(i=1,2)$, $H_i(t)$ $(i=1,2)$ are positive continuous $T$-periodic
functions, $c_i(t)$ $(i=1,2)$
are nonnegative continuous $T$-periodic functions.
\item[(A2)] $a_{31}^{l}-\bar{a}_3>0$. 
\item[(A3)] $\frac{k_i^l}{k_i^l+c_i^uN_1}(\frac{a_i}{k_i})^l>D_i^u
+b_i+2\sqrt{(\frac{a_{ii}}{k_i})^uH_i^u}$ $(i=1,2)$.
\item[(A4)] $H_i^{l}>D_i^uN_1$ $(i=1,2)$.
\end{itemize}

For further convenience, we introduce the 12 positive numbers:
\begin{gather*}
l_i^{\pm}=\frac{(\frac{a_i}{k_i})^u\pm\sqrt{[(\frac{a_i}{k_i})^u]^2
-\frac{4 k_i^l}{k_i^l+c_i^uN_1}(\frac{a_{ii}}{k_i})^{l}(H_i^l-D_i^uN_1)}}
{\frac{2k_i^l}{k_i^l+c_i^uN_1}(\frac{a_{ii}}{k_i})^{l}},
\\
u_i^{\pm}=\frac{[\frac{k_i^l}{k_i^l+c_i^uN_1}(\frac{a_i}{k_i})^l-D_i^u-b_i]
\pm\sqrt{[\frac{k_i^l}{k_i^l+c_i^uN_1}(\frac{a_i}{k_i})^l-D_i^u-b_i]^2
-4(\frac{a_{ii}}{k_i})^uH_i^u}}{2(\frac{a_{ii}}{k_i})^u},
\\
x_i^{\pm}=\frac{\overline{(\frac{a_i}{k_i})}\pm
\sqrt{[\overline{(\frac{a_i}{k_i})}]^2-4\overline{(\frac{a_{ii}}{k_i})}\bar
H_{i}}}{2\overline{(\frac{a_{ii}}{k_i})}},  i=1,2.
\end{gather*}
It is not difficult to show that
\begin{equation} \label{e2.1}
l_i^{-}<x_i^{-}<u_i^{-}<u_i^{+}<x_i^{+}<l_i^{+}\quad (i=1,2)
\end{equation}
Now, we are ready to state the main result of this article.

\begin{theorem} \label{thm2.2}
 Assume that {\rm (A1)-(A4)} hold, then
 \eqref{e1.1} has at least four positive $T$-periodic solutions.
\end{theorem}

\begin{proof}
Since we are concerning with positive
solutions of \eqref{e1.1}, we make the change of variables,
\[
 x_j(t)=e^{u_j(t)}\;(j=1,2),\quad y(t)=e^{u_3(t)}.
\]
Then \eqref{e1.1} is rewritten as
\begin{equation} \label{e2.2}
\begin{aligned}
u_1'(t)
&=\frac{1}{k_1(t)+c_1(t)e^{u_1(t)}}
\Big[a_1(t)-a_{11}(t)e^{u_1(t)}-\frac{a_{13}(t)e^{u_1(t)}e^{u_3(t)}}
{m(t)e^{2u_3(t)}+e^{2u_1(t)}}\Big]\\
&\quad + D_1(t)[\frac{e^{u_2(t)}}{e^{u_1(t)}}-1]-\frac{H_1(t)}{e^{u_1(t)}},\\
u_2'(t)&=\frac{1}{k_2(t)+c_2(t)e^{u_2(t)}}
[a_2(t)-a_{22}(t)e^{u_2(t)}]
 + D_2(t)[\frac{e^{u_1(t)}}{e^{u_2(t)}}-1]-\frac{H_2(t)}{e^{u_2(t)}},\\
u_3'(t)&=-a_3(t)+\frac{a_{31}(t)e^{2u_1(t-\tau)}}{m(t)e^{2u_3(t-\tau)}
+e^{2u_1(t-\tau)}}.
\end{aligned}
\end{equation}
Take
\[ 
X=Z=\{u=(u_1, u_2, u_3)^T\in C(R, R^3):u_i(t+T)=u_i(t), i=1,2,3\}
\]
and define
\[
\|u\|=\max_{t\in[0,T]}|u_1(t)|+\max_{t\in[0,T]}|u_2(t)|+\max_{t\in[0,T]}|u_3(t)|,
\]
for $u=(u_1, u_2, u_3)^T$ in $X$ or in $Z$.
 Equipped with the above norm $\|\cdot\|$, it is easy to verify that
 $X$ and $Z$ are Banach  spaces.

Set
\begin{gather*}
\begin{aligned}
\Delta_{1}(u,t,\lambda)
&= \Big[\frac{k_1(t)+(1-\lambda)
c_1(t)e^{u_1(t)} }{k_1(t)+c_1(t)e^{u_1(t)}}\Big]
\Big[\frac{a_1(t)}{k_1(t)}-\frac{a_{11}(t)e^{u_1(t)}}{k_1(t)}\\
&\quad -\frac{\lambda
k_1(t)^{-1}a_{13}(t)e^{u_1(t)}e^{u_3(t)}}{m(t)e^{2u_3(t)}+e^{2u_1(t)}}\Big] 
+\lambda D_1(t)\big[\frac{e^{u_2(t)}}{e^{u_1(t)}}-1\big]-\frac{H_1(t)}{e^{u_1(t)}},
\end{aligned}
\\
\begin{aligned}
\Delta_{2}(u,t,\lambda)
&=\Big[\frac{k_2(t)+(1-\lambda)c_2(t)e^{u_2(t)}}{k_2(t)+c_2(t)e^{u_2(t)}}\Big]
\Big[\frac{a_2(t)}{k_2(t)}-\frac{a_{22}(t)e^{u_2(t)}}{k_2(t)}\Big]\\
&\quad +\lambda D_2(t)\big[\frac{e^{u_1(t)}}{e^{u_2(t)}}-1\big]
-\frac{H_2(t)}{e^{u_2(t)}},
\end{aligned}
\\
\Delta_{3}(u,t,\lambda)=
-a_3(t)+\frac{a_{31}(t)e^{2u_1(t-\tau)}}{m(t)e^{2u_3(t-\tau)}+\lambda
e^{2u_1(t-\tau)}}.
\end{gather*}
For any  $u\in X$, because of the
periodicity, we can easily check that
 $\Delta_{i}(u,t,\lambda)\in C(R^2,R)\;(i=1,2,3)$ are
$T$-periodic in $t$.

Let
\begin{gather*}
 L:\operatorname{dom}L=\{u\in X:u\in C(R,R^3)\}\ni u\mapsto u'\in Z,\\
 P:X\ni u\mapsto\frac{1}{T}\int_0^Tu(t)dt\in X,\\
 Q:Z\ni u\mapsto\frac{1}{T}\int_0^Tu(t)dt\in Z, \\
 N:X\times[0,1]\ni(u,\lambda)\mapsto(\Delta_1(u,t,\lambda),
\Delta_2(u,t,\lambda), \Delta_3(u,t,\lambda))^T\in Z.
\end{gather*}
 Here, for any $k\in R^3$, we also identify it as
the constant function in $X$ or $Z$ with the constant value $k$. It
is easy to see that
    \begin{gather*}
\ker L=R^3,\quad \operatorname{Im}L=\{u\in X:\int_0^Tu_i(t)dt=0,  i=1,2,3\}
    \end{gather*}
is closed in $Z$, $\dim\ker L=\operatorname{codim}\operatorname{Im}L=3$, and
 $P, Q$ are continuous
projectors such that 
$$
ImP=\ker L, \operatorname{Im}L=\ker Q=Im(I-Q).
$$ 
Therefore, $L$ is a Fredholm mapping of index zero. On the other
hand, $K_p:\operatorname{Im}L\mapsto \operatorname{dom}L\cap KerP$ has the form
\begin{gather*}
K_p(u)=\int_0^{t}u(s)ds-\frac{1}{T}\int_0^T\int_0^{t}u(s)\,ds\,dt.\end{gather*}
Thus,
\begin{gather*}
QN(u,\lambda)=\Big(\frac{1}{T}\int_0^T\Delta_1(u,t,\lambda)dt,
\frac{1}{T}\int_0^T\Delta_2(u,t,\lambda)dt,
\frac{1}{T}\int_0^T\Delta_3(u,t,\lambda)dt\Big)^T,
\\
K_p(I-Q)N(u,\lambda)=(\Phi_1(u,t,\lambda), \Phi_2(u,t,\lambda),
\;\Phi_3(u,t,\lambda))^T,
\end{gather*}
where
\begin{align*}
\Phi_j(u,t,\lambda)
&=\int_0^{t}\Delta_j(u,s,\lambda)ds-
\frac{1}{T}\int_0^T\int_0^{t}\Delta_j(u,s,\lambda)\,ds\,dt\\
&\quad -(\frac{t}{T}-\frac{1}{2})\int_0^T\Delta_j(u,s,\lambda)ds,\quad
j=1,2,3.\end{align*}
Obviously, $QN$ and $K_p(I-Q)N$ are continuous. By the Arzela-Ascoli
theorem, it is not  difficult to show that the closure of
$K_p(I-Q)N(\overline{\Omega}\times[0,1])$ is compact for
any open bounded set $\Omega\subset X$. Moreover,
$QN(\bar\Omega\times[0,1])$ is bounded. Thus $N$ is $L$-compact on
$\bar\Omega\times[0,1]$ with any open bounded set $\Omega\subset X$.

 To apply Lemma \ref{lem2.1}, we need to find at least four
appropriate open, bounded subsets
$\Omega_1,\Omega_2,\Omega_3,\Omega_4$ in $X$.


Corresponding to the operator equation $Lu=\lambda
N(u,\lambda)$, $\lambda\in(0,1)$, we have
\begin{equation} \label{e2.3}
\begin{aligned}
u_1'(t)
&= \lambda\Big[\frac{k_1(t)+(1-\lambda)
c_1(t)e^{u_1(t)} }{k_1(t)+c_1(t)e^{u_1(t)}}\Big]
\Big[\frac{a_1(t)}{k_1(t)}-\frac{a_{11}(t)e^{u_1(t)}}{k_1(t)}\\
&\quad -\frac{\lambda
k_1(t)^{-1}a_{13}(t)e^{u_1(t)}e^{u_3(t)}}{m(t)e^{2u_3(t)}+e^{2u_1(t)}}\Big]
+\lambda^2 D_1(t)\big[\frac{e^{u_2(t)}}{e^{u_1(t)}}-1\big]-\frac{\lambda
H_1(t)}{e^{u_1(t)}},
\end{aligned}
\end{equation}
\begin{equation} \label{e2.4}
\begin{aligned}
u_2'(t)&=\lambda\Big[\frac{
k_2(t)+(1-\lambda)c_2(t)e^{u_2(t)}}{k_2(t)+c_2(t)e^{u_2(t)}}\Big]
\Big[\frac{a_2(t)}{k_2(t)}-\frac{a_{22}(t)e^{u_2(t)}}{k_2(t)}\Big]\\
&\quad +\lambda^2 D_2(t)\Big[\frac{e^{u_1(t)}}{e^{u_2(t)}}-1\Big]
-\frac{\lambda H_2(t)}{e^{u_2(t)}},
\end{aligned}
\end{equation}
\begin{equation} \label{e2.5}
u_3'(t)=\lambda\Big[-a_3(t)+\frac{a_{31}(t)e^{2u_1(t-\tau)}}{m(t)e^{2u_3(t-\tau)}
+\lambda e^{2u_1(t-\tau)}}\Big].
\end{equation}

Suppose that $(u_1(t), u_2(t), u_3(t))^T$ is a $T$-periodic
solution of \eqref{e2.3}, \eqref{e2.4} and \eqref{e2.5} for some $\lambda\in(0,1)$.
Choose $t_i^M, t_i^m\in[0,T]$, $i=1,2,3$, such that
$$
u_i(t_i^M)=\max_{t\in[0,T]}u_i(t),\quad
u_i(t_i^m)=\min_{t\in[0,T]}u_i(t), \quad i=1,2,3.
$$
Then, it is clear that
\[
u_i'(t_i^M)=0, \quad u_i'(t_i^m)=0, \quad i=1,2,3.
\]
From this and \eqref{e2.3}, \eqref{e2.4}, we obtain that
\begin{equation} \label{e2.6}
\begin{aligned}
0&=\Big[\frac{k_1(t_1^M)+(1-\lambda) c_1(t_1^M)e^{u_1(t_1^M)}
}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\Big]
\Big[\frac{a_1(t_1^M)-a_{11}(t_1^M)e^{u_1(t_1^M)}}{k_1(t_1^M)} \\
&\quad -\frac{\lambda k_1(t_1^M)^{-1}a_{13}(t_1^M)
 e^{u_1(t_1^M)}e^{u_3(t_1^M)}}{m(t_1^M)e^{2u_3(t_1^M)}+e^{2u_1(t_1^M)}}\Big]
+\lambda D_1(t_1^M)\big[\frac{e^{u_2(t_1^M)}}{e^{u_1(t_1^M)}}-1\big]
-\frac{H_1(t_1^M)}{e^{u_1(t_1^M)}},
\end{aligned}
\end{equation}
\begin{equation} \label{e2.7}
\begin{aligned}
0&=\Big[\frac{k_2(t_2^M)+(1-\lambda) c_2(t_2^M)e^{u_2(t_2^M)}
}{k_2(t_2^M)+c_2(t_2^M)e^{u_2(t_2^M)}}\Big]
 \Big[\frac{a_2(t_2^M)}{k_2(t_2^M)}
 -\frac{a_{22}(t_2^M)e^{u_2(t_2^M)}}{k_2(t_2^M)}\Big]
\\
&\quad +\lambda D_2(t_2^M)\big[\frac{e^{u_1(t_2^M)}}{e^{u_2(t_2^M)}}-1\big]
-\frac{H_2(t_2^M)}{e^{u_2(t_2^M)}},
\end{aligned}
\end{equation}
\begin{equation} \label{e2.8}
\begin{aligned}
0&=\Big[\frac{k_1(t_1^m)+(1-\lambda) c_1(t_1^m)e^{u_1(t_1^m)}
}{k_1(t_1^m)+c_1(t_1^m)e^{u_1(t_1^m)}}\Big]
\Big[\frac{a_1(t_1^m)-a_{11}(t_1^m)e^{u_1(t_1^m)}}{k_1(t_1^m)} \\
&\quad -\frac{\lambda
k_1(t_1^m)^{-1}a_{13}(t_1^m)e^{u_1(t_1^m)}
 e^{u_3(t_1^m)}}{m(t_1^m)e^{2u_3(t_1^m)}+e^{2u_1(t_1^m)}}\Big]
+\lambda D_1(t_1^m)\big[\frac{e^{u_2(t_1^m)}}{e^{u_1(t_1^m)}}-1\big]
-\frac{H_1(t_1^m)}{e^{u_1(t_1^m)}},
\end{aligned}
\end{equation}
\begin{equation} \label{e2.9}
\begin{aligned}
0&=\Big[\frac{k_2(t_2^m)+(1-\lambda) c_2(t_2^m)e^{u_2(t_2^m)}
}{k_2(t_2^m)+c_2(t_2^m)e^{u_2(t_2^m)}}\Big]
\Big[\frac{a_2(t_2^m)}{k_2(t_2^m)}-\frac{a_{22}(t_2^m)e^{u_2(t_2^m)}}{k_2(t_2^m)}
\Big] \\
&\quad +\lambda D_2(t_2^m)\big[\frac{e^{u_1(t_2^m)}}{e^{u_2(t_2^m)}}-1\big]
-\frac{H_2(t_2^m)}{e^{u_2(t_2^m)}}.
\end{aligned}
\end{equation}
For $u(t_i^M)\;(i=1,2)$, there are two cases to consider.

\noindent  \textbf{Case 1.} Assume that $u_1(t_1^M)\geq u_2(t_2^M)$, then
  $u_1(t_1^M)\geq u_2(t_1^M)$.
From this and \eqref{e2.6}, we have
\[
a_1(t_1^M)-a_{11}(t_1^M)e^{u_1(t_1^M)}>0,
\]
which implies
\[
e^{u_1(t_1^M)}<\frac{a_1(t_1^M)}{a_{11}(t_1^M)}\leq (\frac{a_1}{a_{11}})^u\leq N_1.
\]
That is
\[
u_2(t_2^M)\leq u_1(t_1^M)<\ln(\frac{a_1}{a_{11}})^u\leq\ln
N_1.
\]

\noindent \textbf{Case 2.} Assume that $u_1(t_1^M)<u_2(t_2^M)$, then
$u_2(t_2^M)>u_1(t_2^M)$.
From this and \eqref{e2.7}, we have
\[
a_2(t_2^M)-a_{22}(t_2^M)e^{u_2(t_2^M)}> 0,
\]
which implies
\[
e^{u_2(t_2^M)}<\frac{a_2(t_2^M)}{a_{22}(t_2^M)}\leq
(\frac{a_2}{a_{22}})^u\leq N_1.
\]
That is,
\[
u_1(t_1^M)<u_2(t_2^M)<\ln(\frac{a_2}{a_{22}})^u\leq\ln N_1.
\]
Therefore,
\begin{equation} \label{e2.10}
\max\{u_1(t_1^M),u_2(t_2^M)\}<\ln N_1.
\end{equation}
It follows from \eqref{e2.6} that
\begin{align*} 
&\Big[\frac{k_1(t_1^M)+(1-\lambda) c_1(t_1^M)e^{u_1(t_1^M)}
}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\Big]
\Big[\frac{a_1(t_1^M)}{k_1(t_1^M)}-\frac{a_{11}(t_1^M)e^{u_1(t_1^M)}}{k_1(t_1^M)}
\Big]-\frac{H_1(t_1^M)}{e^{u_1(t_1^M)}}\\
&=  \lambda \Big[\frac{k_1(t_1^M)+(1-\lambda)c_1(t_1^M)e^{u_1(t_1^M)}
}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\Big]
\frac{k_1(t_1^M)^{-1}a_{13}(t_1^M)e^{u_1(t_1^M)}e^{u_3(t_1^M)}}
 {m(t_1^M)e^{2u_3(t_1^M)}+e^{2u_1(t_1^M)}}\\
&\quad -\lambda D_1(t_1^M)\frac{e^{u_2(t_1^M)}}{e^{u_1(t_1^M)}}+\lambda
D_1(t_1^M).
\end{align*}
Therefore,
\begin{align*}  
&\Big[\frac{k_1(t_1^M)+(1-\lambda) c_1(t_1^M)e^{u_1(t_1^M)}
}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\Big]
\Big[(\frac{a_1}{k_1})^l-(\frac{a_{11}}{k_1})^u
e^{u_1(t_1^M)}\Big]-\frac{H_1^u}{e^{u_1(t_1^M)}}\\
&<\Big[\frac{k_1(t_1^M)+(1-\lambda) c_1(t_1^M)e^{u_1(t_1^M)}
}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\Big]
\frac{a_{13}(t_1^M)}{2k_1(t_1^M)\sqrt{m(t_1^M)}}+D_1(t_1^M),
\end{align*}
which implies
\begin{align*}  
&\Big[\frac{k_1(t_1^M)+(1-\lambda) c_1(t_1^M)e^{u_1(t_1^M)}
}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\Big]
\Big[(\frac{a_1}{k_1})^l-(\frac{a_{11}}{k_1})^u
e^{u_1(t_1^M)}\Big]-\frac{H_1^u}{e^{u_1(t_1^M)}}\\
&<\Big[\frac{k_1(t_1^M)+(1-\lambda) c_1(t_1^M)e^{u_1(t_1^M)}
}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\Big]
\frac{(\frac{a_{13}}{k_1})^u}{2\sqrt{m^l}}+D_1^u.
\end{align*}
From this and noticing that
\[
\frac{k_1(t_1^M)}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\leq
\frac{k_1(t_1^M)+(1-\lambda) c_1(t_1^M)e^{u_1(t_1^M)}
}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\leq 1,
\]
we have
\[
\frac{k_1(t_1^M)}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}
(\frac{a_1}{k_1})^l-(\frac{a_{11}}{k_1})^u
e^{u_1(t_1^M)}-\frac{H_1^u}{e^{u_1(t_1^M)}}
<\frac{(\frac{a_{13}}{k_1})^u}{2\sqrt{m^l}}+D_1^u.
\]
Therefore, we have
\[
\frac{k_1^l}{k_1^l+c_1^uN_1}(\frac{a_1}{k_1})^l-(\frac{a_{11}}{k_1})^u
e^{u_1(t_1^M)}-\frac{H_1^u}{e^{u_1(t_1^M)}}
<\frac{(\frac{a_{13}}{k_1})^u}{2\sqrt{m^l}}+D_1^u;
\]
that is,
\[
(\frac{a_{11}}{k_1})^u e^{2u_1(t_1^M)}
-\Big[\frac{k_1^l}{k_1^l+c_1^uN_1}(\frac{a_1}{k_1})^l-D_1^u-b_1\Big]
e^{u_1(t_1^M)}+H_1^u>0.
\]
From (A3) and the above inequality, we have
 \begin{equation} \label{e2.11}
u_1(t_1^M)>\ln u_1^{+}\quad \mbox{or}\quad  u_1(t_1^M)<\ln u_1^{-}
 \end{equation}
Similarly, we have
\begin{equation} \label{e2.12}
u_1(t_1^m)>\ln u_1^{+}\quad\mbox{or}\quad  u_1(t_1^m)<\ln u_1^{-}
\end{equation}
 By a similar argument, from \eqref{e2.7}, it follows that
\[
(\frac{a_{22}}{k_2})^u e^{2u_2(t_2^M)}
-\Big[\frac{k_2^l}{k_2^l+c_2^uN_1}(\frac{a_2}{k_2})^l-D_2^u-b_2\Big]
e^{u_2(t_2^M)}+H_2^u>0.
\]
By (A3) and the above inequality, we have
 \begin{equation} \label{e2.13}
u_2(t_2^M)>\ln u_2^{+}\quad \mbox{or}\quad   u_2(t_2^M)<\ln u_2^{-}.
 \end{equation}
Similarly, we have
\begin{equation} \label{e2.14}  
u_2(t_2^m)>\ln u_2^{+}\quad \mbox{or}\quad u_2(t_2^m)<\ln u_2^{-}.
\end{equation}
Again, from \eqref{e2.6} it follows that
\begin{align*}
&\Big[\frac{k_1(t_1^M)+(1-\lambda) c_1(t_1^M)e^{u_1(t_1^M)}
}{k_1(t_1^M)+c_1(t_1^M)e^{u_1(t_1^M)}}\Big]
\Big[\frac{a_1(t_1^M)}{k_1(t_1^M)}-\frac{a_{11}(t_1^M)e^{u_1(t_1^M)}}{k_1(t_1^M)}
\Big]\\
&>\frac{H_1(t_1^M)}{e^{u_1(t_1^M)}}-
D_1(t_1^M)\frac{e^{u_2(t_1^M)}}{e^{u_1(t_1^M)}}.
\end{align*}
Hence, we have
\[
\frac{k_1^l}{k_1^l+c_1^uN_1}\frac{a_{11}(t_1^M)}{k_1(t_1^M)}e^{2u_1(t_1^M)}
-\frac{a_1(t_1^M)}{k_1(t_1^M)}e^{u_1(t_1^M)}+H_1(t_1^M)- D_1(t_1^M)e^{u_2(t_1^M)}<0,
\]
which implies
\[
\frac{k_1^l}{k_1^l+c_1^uN_1}(\frac{a_{11}}{k_1})^{l} e^{2u_1(t_1^M)}
-(\frac{a_1}{k_1})^u e^{u_1(t_1^M)}+H_1^l-D_1^u N_1<0.
\]
From (A3), (A4) and the above inequality, we have
  \begin{equation} \label{e2.15}
\ln l_1^{-}<u_1(t_1^M)<\ln l_1^{+}.
  \end{equation}
Similarly, we have
 \begin{equation} \label{e2.16}
\ln l_1^{-}<u_1(t_1^m)<\ln l_1^{+}.
 \end{equation}
By a similar argument, it follows from \eqref{e2.7} that
\[
\frac{k_2^l}{k_2^l+c_2^uN_1}(\frac{a_{22}}{k_2})^{l} e^{2u_2(t_2^M)}
-(\frac{a_2}{k_2})^u e^{u_2(t_2^M)}+H_2^l-D_2^u N_1<0.
\]
From (A3), (A4) and the above inequality, we have
\begin{equation} \label{e2.17}
 \ln l_2^{-}<u_2(t_2^M)<\ln l_2^{+}.
\end{equation}
Similarly, we have
      \begin{equation} \label{e2.18}  
\ln l_2^{-}<u_2(t_2^m)<\ln l_2^{+}.
      \end{equation}
It follows from \eqref{e2.11}, \eqref{e2.12}, \eqref{e2.15}, \eqref{e2.16} that
\begin{gather*}
u_1(t_1^M)\in(\ln l_1^{-},\ln u_1^{-})\cup(\ln u_1^{+},\ln l_1^{+}),
\\
u_1(t_1^m)\in(\ln l_1^{-},\ln u_1^{-})\cup(\ln u_1^{+},\ln l_1^{+}).
 \end{gather*}
It follows from \eqref{e2.13}, \eqref{e2.14}, \eqref{e2.17}, \eqref{e2.18} that
 \begin{gather*}
u_2(t_2^M)\in(\ln l_2^{-},\ln u_2^{-})\cup(\ln u_2^{+},\ln l_2^{+}),
\\
u_2(t_2^m)\in(\ln l_2^{-},\ln u_2^{-})\cup(\ln u_2^{+},\ln l_2^{+}).
 \end{gather*}
From \eqref{e2.5}, we have
\begin{equation} \label{e2.19}
\int_0^T a_3(t) dt
=\int_0^T\frac{a_{31}(t)e^{2u_1(t-\tau)}}{m(t)e^{2u_3(t-\tau)}+\lambda
e^{2u_1(t-\tau)}}dt.
\end{equation}
Therefore,
\begin{equation} \label{e2.20}
\begin{aligned}
 \int _0^T|u_3'(t)|dt
 &<\int_0^Ta_3(t)dt+\int_0^T\frac{a_{31}(t)e^{2u_1(t-\tau)}}{m(t)e^{2u_3(t-\tau)}+\lambda
 e^{2u_1(t-\tau)}}dt\\
 &=\int_0^Ta_3(t)dt+\int_0^Ta_3(t)dt\\
 &<2T \bar {a}_3:=d_3.
\end{aligned}
\end{equation}
By \eqref{e2.19}, there must be $\eta\in(0,T)$ such that
 \begin{equation} \label{e2.21}
\bar {a}_3=\frac{a_{31}(\eta)e^{2u_1(\eta-\tau)}}{m(\eta)e^{2u_3(\eta-\tau)}
+\lambda e^{2u_1(\eta-\tau)}}.
 \end{equation}
From \eqref{e2.21}, we have
\[ 
e^{2u_3(\eta-\tau)}=\frac{a_{31}(\eta)
-\lambda \bar {a}_3}{\bar {a}_3 m(\eta)}e^{2u_1(\eta-\tau)}.
\]
Therefore,
 \begin{gather} \label{e2.22} 
e^{2u_3(\eta-\tau)}<\frac{a_{31}^u}{\bar {a}_3m^l}N_1^2:=p, \\
 \label{e2.23} 
e^{2u_3(\eta-\tau)}>\frac{a_{31}^l-\bar {a}_3}{\bar {a}_3m^u}(l_1^{-})^2:=q.
 \end{gather}
Since $u_3(t)$ is a $T$-periodic function, there exists
$\eta_1\in[0,T]$ such that
\[
u_3(\eta_1)=u_3(\eta-\tau).
\]
Again, it follows from \eqref{e2.20}, \eqref{e2.22}, \eqref{e2.23} that for any
$t\in[0,T]$, we have
\begin{gather*}
u_3(t)=u_3(\eta_1)+\int_{\eta_1}^{t}u_3'(s)ds<\frac{1}{2}\ln p+d_3, \\
u_3(t)=u_3(\eta_1)+\int_{\eta_1}^{t}u_3'(s)ds>\frac{1}{2}\ln q-d_3.
\end{gather*}
Set 
\[
H=\max\{|\frac{1}{2}\ln p+d_3|,|\frac{1}{2}\ln q-d_3|\}.
\]
Then
\begin{equation} \label{e2.24}
\max_{t\in[0,T]}|u_3(t)|<H.
\end{equation}
 Clearly, $l_i^{\pm},u_i^{\pm}(i=1,2),H$ are independent of  $\lambda$. 
Now, let us consider $QN(u,0)$ with
$u=(u_1, u_2, u_3)^T\in R^3$. Note that
\[
QN(u,0)=\begin{pmatrix}
 \overline{(\frac{a_1}{k_1})}-\overline{(\frac{a_{11}}{k_1})}e^{u_1}-\frac{\bar
H_1}{e^{u_1}}\\
\overline{(\frac{a_2}{k_2})}-\overline{(\frac{a_{22}}{k_2})}e^{u_2}-\frac{\bar H_2}{e^{u_2}}\\
-\bar
{a}_3+\overline{(\frac{a_{31}}{m})}\frac{e^{2u_1}}{e^{2u_3}}
\end{pmatrix}.
\]
Letting $QN(u,0)=0$, we have 
\begin{gather} \label{e2.25}
\overline{(\frac{a_1}{k_1})}-\overline{(\frac{a_{11}}{k_1})}e^{u_1}-\frac{\bar
H_1}{e^{u_1}}=0, \\
 \label{e2.26}
\overline{(\frac{a_2}{k_2})}-\overline{(\frac{a_{22}}{k_2})}e^{u_2}
-\frac{\bar H_2}{e^{u_2}}=0, \\
 \label{e2.27}
-\bar {a}_3+\overline{(\frac{a_{31}}{m})}\frac{e^{2u_1}}{e^{2u_3}}=0.
\end{gather}
From \eqref{e2.25}, we have
\begin{equation} \label{e2.28}
\overline{(\frac{a_{11}}{k_1})}e^{2u_1}-\overline{(\frac{a_1}{k_1})}e^{u_1}
+\bar H_1=0,
\end{equation}
which implies that \eqref{e2.28} has two distinct roots $\ln x_1^{\pm}$.
From \eqref{e2.26}, we have
\begin{equation} \label{e2.29}
\overline{(\frac{a_{22}}{k_2})}e^{2u_2}-\overline{(\frac{a_2}{k_2})}e^{u_2}
+\bar H_2=0,
\end{equation}
which implies that \eqref{e2.29} has two distinct roots $\ln x_2^{\pm}$.
Again, it follows from \eqref{e2.27} that
\[
x_3^{\pm}=\sqrt{\frac{\overline{(\frac{a_{31}}{m})}}{\bar{a}_3}}x_1^{\pm}.
\]
Therefore, $QN(u,0)=0$ has four distinct solutions.
\begin{gather*}
\tilde{u}_1=(\ln x_1^{+},\ln x_2^{+},\ln x_3^{+})^T,\quad
\tilde{u}_2=(\ln x_1^{+},\ln x_2^{-},\ln x_3^{+})^T, \\
\tilde{u}_3=(\ln x_1^{-},\ln x_2^{+},\ln x_3^{-})^T,\quad 
\tilde{u}_4=(\ln x_1^{-},\ln x_2^{-},\ln x_3^{-})^T.
\end{gather*}
Choose $C>0$ such that
\begin{equation} \label{e2.30}
C>\max\Big\{|\ln \sqrt{\frac{\overline{(\frac{a_{31}}{m})}}{\bar{a}_3}}x_1^{+}|, 
|\ln \sqrt{\frac{\overline{(\frac{a_{31}}{m})}}{\bar{a}_3}}x_1^{-}|\Big\}.
\end{equation}
Let
\begin{align*}
\Omega_1=\Big\{&u=(u_1,u_2,u_3)^T\in X: 
\max_{t\in[0,T]}u_1(t)\in(\ln u_1^{+}, \ln l_1^{+}), \\
&\min_{t\in[0,T]}u_1(t)\in(\ln u_1^{+},  \ln l_1^{+}),\;
\max_{t\in[0,T]}u_2(t)\in(\ln u_2^{+}, \ln l_2^{+}), \\
&\min_{t\in[0,T]}u_2(t)\in(\ln u_2^{+}, \ln l_2^{+}),\;
\max_{t\in[0,T]}|u_3(t)|<H+C\Big\},
\end{align*}
\begin{align*}
\Omega_2=\Big\{&u=(u_1,u_2,u_3)^T\in X: 
\max_{t\in[0,T]}u_1(t)\in(\ln u_1^{+}, \ln l_1^{+}), \\
&\min_{t\in[0,T]}u_1(t)\in(\ln u_1^{+},  \ln l_1^{+}),\;
\max_{t\in[0,T]}u_2(t)\in(\ln l_2^{-}, \ln u_2^{-}), \\
&\min_{t\in[0,T]}u_2(t)\in(\ln l_2^{-}, \ln u_2^{-}),\;
\max_{t\in[0,T]}|u_3(t)|<H+C\Big\},
\end{align*}
\begin{align*}
\Omega_3=\Big\{&u=(u_1,u_2,u_3)^T\in X:
\max_{t\in[0,T]}u_1(t)\in(\ln l_1^{-}, \ln u_1^{-}), \\
&\min_{t\in[0,T]}u_1(t)\in(\ln l_1^{-},  \ln u_1^{-}),\;
\max_{t\in[0,T]}u_2(t)\in(\ln u_2^{+}, \ln l_2^{+}), \\
&\min_{t\in[0,T]}u_2(t)\in(\ln u_2^{+}, \ln l_2^{+}),\;
\max_{t\in[0,T]}|u_3(t)|<H+C\Big\},
\end{align*}
\begin{align*}
\Omega_4=\Big\{&u=(u_1,u_2,u_3)^T\in X:
\max_{t\in[0,T]}u_1(t)\in(\ln l_1^{-}, \ln u_1^{-}), \\
&\min_{t\in[0,T]}u_1(t)\in(\ln l_1^{-},  \ln u_1^{-}),\;
\max_{t\in[0,T]}u_2(t)\in(\ln l_2^{-}, \ln u_2^{-}), \\
&\min_{t\in[0,T]}u_2(t)\in(\ln l_2^{-}, \ln u_2^{-}),\;
\max_{t\in[0,T]}|u_3(t)|<H+C\Big\}.
\end{align*}
Then $\Omega_1,\dots,\Omega_4$ are bounded open subsets of $X$.
It follows from \eqref{e2.1} and \eqref{e2.30} that
$\tilde{u}_i\in\Omega_i$ $(i=1,2,3,4)$. It is easy to see that
$\bar{\Omega}_i\cap\bar{\Omega}_j=\emptyset$ $(i,j=1,2,3,4,i\neq j)$ and
 $\Omega_i$ satisfies (a) in Lemma \ref{lem2.1} for $i=1,2,3,4$.
Moreover, $QN(u,0)\neq0$ for $u\in \partial\Omega_i\cap \ker L$.

Set
\[
h(x)=b-ax-\frac{c}{x},\quad  x\in (0, +\infty),
\]
where $a,b,c$ are positive constants.
When $b>2\sqrt{ac}$, it is easy to see
that $h(x)=0$  has exactly two positive solutions
\[
x^{-}=\frac{b-\sqrt{b^2-4ac}}{2a},\quad
x^{+}=\frac{b+\sqrt{b^2-4ac}}{2a}
\]
such that $h'(x^{-})>0$, $h'(x^{+})<0$.
From this, a direct computation gives
\begin{gather*}
\deg\{JQN(\cdot,0),\Omega_1\cap \ker L,0\}=-1, \quad
\deg\{JQN(\cdot,0),\Omega_2\cap \ker L,0\}=1, \\
\deg\{JQN(\cdot,0),\Omega_3\cap \ker L,0\}= 1, \quad
\deg\{JQN(\cdot,0),\Omega_4\cap \ker L,0\}=-1.
\end{gather*}
Here, $J$ is taken as the identity mapping since $\operatorname{Im}Q=\ker L$. 
So far we have proved that $\Omega_i$ satisfies all the assumptions in
Lemma \ref{lem2.1}. Hence, \eqref{e2.2} has at least four $T$-periodic solutions
$(u_1^{i}(t), u_2^{i}(t),u_3^{i}(t))^T$ $(i=1,2,3,4)$ and
 $(u_1^{i}, u_2^{i},u_3^{i})^T\in \operatorname{dom}L\cap\bar{\Omega}_i$.
Obviously, $(u_1^{i}, u_2^{i},u_3^{i})^T$ $(i=1,2,3,4)$ are different. Let
$x_j^{i}(t)=e^{u_j^{i}(t)}$ $(j=1,2)$, $y^{i}(t)=e^{u_3^{i}(t)}$ $(i=1,2,3,4)$. Then
$(x_1^{i}(t), x_2^{i}(t),y^i(t))^T$ $(i=1,2,3,4)$ are four different positive 
$T$-periodic solutions of \eqref{e1.1}.  The proof is complete. 
\end{proof}

\begin{theorem} \label{thm2.3}
In addition to {\rm (A1), (A2), (A4),}
   assume further that  \eqref{e1.1} satisfies
\begin{itemize}
\item[(A3)*] $H_i^u<\frac{N_1}{N_1+1}
\big[\frac{k_i^l}{k_i^l+c_i^uN_1}(\frac{a_i}{k_i})^l-(\frac{a_{ii}}{k_i})^u-b_i
\big]$ $(i=1,2)$.
\end{itemize}
Then \eqref{e1.1} has at least four positive $T$-periodic
solutions.
\end{theorem}

\begin{proof} 
It suffices to verify (A3) in Theorem \ref{thm2.2}. Indeed,
it follows from (A3)* and (A4) that
\[
D_i^u+b_i+2\sqrt{(\frac{a_{ii}}{k_i})^uH_i^u}
<\frac{H_i^l}{N_1}+b_i+(\frac{a_{ii}}{k_i})^u+H_i^u
<\frac{k_i^l}{k_i^l+c_i^uN_1}(\frac{a_i}{k_i})^l,
\]
for $i=1,2$.
\end{proof}

Finally, we describe the biological meaning of (A1)--(A4)
and Theorem \ref{thm2.2}. In realistic world, the environment is always
varying periodically with time. This motivates us to consider 
system \eqref{e1.1}  under the condition (A1).
(A2) indicates  that the efficiency of turning predated prey species 
$x$ into predator species $y$ is higher than the natural death rate of predator
species $y$, which is the necessary condition for the survival of
predator species $y$. Note that
$\sqrt{(\frac{a_{ii}}{k_i})^uH_i^u}$ is the geometric mean of
$(\frac{a_{ii}}{k_i})^u$ and $H_i^u$ $(i=1,2)$, which describes
the mean effect of the intra-species competition and harvesting.
Hence, (A3) indicates that the natural growth rate of prey
species $x$ with a food-limited supply in patch 1 is higher than the
decay rate due to the intra-species competition, predation,
dispersion and harvesting, and that the natural growth rate of prey
species $x$ with a limited food supply in patch 2 is higher than the
decay rate due to the intra-species competition, the dispersion and
harvesting. Since $N_1$ is the environmental carrying capacity of
prey species $x$, (A4) indicates that the effect of harvesting
on populations is stronger than the effect of the dispersion during
each time period. From the ecological viewpoint, Theorem \ref{thm2.2}
indicates that a food-limited two-species ratio-dependent
predator-prey patch systems with delay and harvesting can lead to
four different periodic fluctuations in a periodic environment if
the regulated harvesting on prey populations is made according to
(A1)--(A4).


\section{An Example}

Take $\tau=2$, $T=2$, $k_1(t)=4+\sin (\pi t)$, 
$k_2(t)=3+\sin (\pi t)$, $c_i(t)=0.2+0.05\sin(\pi t)$ $(i=1,2)$,
\begin{gather*}
H_1(t)=\frac{1+\sin^2(\pi t)}{24},\quad
H_2(t)=\frac{1+\sin^2(\pi t)}{10}, \quad
D_1(t)=D_2(t)=\frac{1+\sin^2(\pi t)}{200}, \\
a_1(t)=(4+\sin(\pi t))^2, \quad
a_{11}(t)=\frac{(4+\sin(\pi t))^2}{4},\quad 
a_{13}(t)=\frac{(4+\sin(\pi t))^2}{4},\\
m(t)=1+\sin^2(\pi t), \quad
a_2(t)=(3+\sin(\pi t))^2,\quad
a_{22}(t)=\frac{(3+\sin(\pi t))^2}{4}, \\
a_{31}(t)=4+\sin(\pi t),\quad
a_3(t)=2+\sin(\pi t).
\end{gather*}
Then we have
$k_1^l=3$, $k_2^l=2$, $c_i^u=0.25$ $(i=1,2)$,
\begin{gather*} 
H_1^u=\frac{1}{12},\quad H_2^u=\frac{1}{5},\quad
H_1^l=\frac{1}{24},\quad H_2^l=\frac{1}{10},\quad D_1^u=D_2^u=\frac{1}{100},\\
N_1=4, \quad 
(\frac{a_1}{k_1})^l=3,\quad 
(\frac{a_{11}}{k_1})^u=\frac{5}{4},\quad
(\frac{a_{13}}{k_1})^u=\frac{5}{4}, \quad
m^l=1,\\
b_1=\frac{5}{8}, \quad 
(\frac{a_2}{k_2})^l=2,\quad
(\frac{a_{22}}{k_2})^u=1, \quad 
a_{31}^l=3,\quad \bar{a}_3=2.
\end{gather*}
It is easy to see that the conditions in
Theorem \ref{thm2.3} are satisfied. By Theorem \ref{thm2.3}, 
system \eqref{e1.1} has at
least four positive $2$-periodic solutions.


\begin{thebibliography}{10}

\bibitem{b1} F. Brauer, A. Soudack; 
\emph{Coexistence properties of some predator-prey systems under constant 
rate harvesting and stocking},  J. Math. Biol. 12 (1981) 101-114.

\bibitem{b2} F. Brauer, A. Soudack; 
\emph{On constant effort harvesting and stocking in a class of predator-prey
systems}, J. Theoret. Biol. 95 (1982) 247-252.

\bibitem{c1} F. Chen, D. Sun, J. Shi;
 \emph{Periodicity in a food-limited population model with toxicants and
state dependent delays}, J. Math. Anal. Appl. 288 (2003) 136-146.

\bibitem{c2} Y. Chen; 
\emph{Multiple periodic solutions of delayed predator-prey systems
with type IV functional responses}, Nonlinear Anal. Real World Appl.
5 (2004) 45-53.

\bibitem{c3} C. Clark; 
\emph{Mathematical Bioeconomics, The optimal
management of renewable resources}, 2nd Ed. John Wiley \& Sons, New
York-Toronto, USA, 1990.

\bibitem{d} S. Dong, W. Ge; 
\emph{Periodic solutions for two-species ratio-dependent predator-prey patch 
systems with delay}, Acta Math. Appl. Sin. 27 (2004) 132-141.

\bibitem{f1} J. Feng, S. Chen; 
\emph{Four periodic solutions of a generalized delayed predator-prey
system}, Appl. Math. Comput.  181 (2006) 932-939.

\bibitem{g1} R. Gaines, J. Mawhin;
\emph{Coincidence Degree and Nonlinear Differetial Equitions}, 
Springer Verlag, Berlin, 1977.

\bibitem{g2} K. Gopalsamy, M. Kulenovic, G. Ladas; 
\emph{Environmental periodicity and time delays in a food-limited population model}, 
J. Math. Anal. Appl. 147 (1990) 545-555.

\bibitem{g3} S. Gourley, M. Chaplain; 
\emph{Travelling fronts in a food-limited population model with time delay}, 
Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 132A (2002) 75-89.

\bibitem{g4} S. Gourley, J. So;
\emph{Dynamics of a food-limited population model incorporating non local delays 
on a finite domain}, J. Math. Biol. 44 (2002) 49-78.

\bibitem{h} A. Huusko, P. Hyvarinen;
\emph{A high harvest rate induces a tendency to generation cycling in a
freshwater fish population}, Journal of Animal Ecology, 74(3), 2005,
525-531.

\bibitem{li} Y. Li; 
\emph{Periodic solutions of periodic generalized food limited
model}, Int. J. Math. Math. Sci. 25 (2001) 265-271.

\bibitem{x1} Y. Xia, J. Cao, S. Cheng;
\emph{Multiple periodic solutions of a delayed stage-structured
predator-prey model with non-monotone functional responses}, Appl.
Math. Modelling  31  (2007) 1947-1959.

\bibitem{x2} D. Xiao and L. Jennings; 
\emph{Bifurcations of a ratio-dependent predator-prey system with
constant rate harvesting}, SIAM J. Appl. Math. 65 (2005) 737-753.


\bibitem{y} A. Yakubu; 
\emph{Searching predator and prey dominance in discrete predator-prey systems with
dispersion}, SIAM J. Appl. Math. 61 (2000) 870-888.

\bibitem{z} Z. Zhang, T. Tian; 
\emph{Multiple positive periodic solutions for a generalized predator-prey
system with exploited terms}, Nonlinear Anal. Real World Appl. 9 (2007) 26-39.

\end{thebibliography}

\end{document}