\section{Predator-Prey System--Introduction to Singular Orbit Analysis}
There are 2 important ecologically relevant issues: long term
dynamics and short term outbreaks and crashes. Various methods can
be used to analyze these problems, not all equal in effectiveness.
To motivate the geometric method of this paper, we first give a
brief illustration and comparison of these methods. We will do so
in the context of the simplest case of one predator and one prey
system
\begin{equation}\label{eq2dPredatorPrey}
%\begin{array}{c}
\displaystyle{\zeta \frac{dx}{dt} = x\left(1 - x -
\frac{y}{\beta_1 +x}\right), \quad \frac{dy}{dt} = y
\left(\frac{x}{\beta_1+x} - \delta_1\right)},
%\end{array}
\end{equation}
which is the $XY$-subsystem of (\ref{ScaledSystem}) with the
absence of the competitor $Z$ ($Z=0$) and the top-predator $W$
($W=0$).

With regard to the simplest long term dynamics, the system may
have a unique $xy$-equilibrium point, expressed explicitly at
$x=x_{e}=(\beta_1\delta_1)/(1-\delta_1)>0,y=(1-x_e)(\beta_1+x_2)$.
To determine its local stability, one linearizes the system at the
equilibrium point, finds the eigenvalues, which quickly lead to
the following conclusion
\begin{enumerate}
\item The equilibrium state is local stable if and only if the
predator $y$ is weak (\ref{hypoweak}).
\item If $y$ is efficient,
\[
0<\frac{\beta_1\delta_1}{1-\delta_1}<\frac{1-\beta_1}{2},
\]
then a small limit cycle is created surround the unstable
equilibrium point via Hopf bifurcation.
\end{enumerate}
There are a few obvious drawbacks. First, it is only local.
Second, it is only about equilibrium point. Third, little can be
said about short-term temporal outbreaks and crashes. Last,
although theoretically possible, it is practically difficult to
impossible to explicitly solve for higher dimensional systems the
equilibrium points, or their eigenvalues, or to apply the
Routh-Hurwitz criterion (e.g. \cite{Mu89}) or other stability
criteria.

To circumvent some of these algebraic difficulties associated with
stability determination and at the same time to extract just the
right amount of information, there is another elementary but
effective technique that is, surprisingly, rarely seen in the
literature. We give an illustration of this method in Appendix B
because of the additional reason, other than the reason for
exposition, that it will be used in various places later.

The second method is to apply Kolmogoroff's Theorem (c.f.
\cite{M72}) to the system to conclude
\begin{enumerate}
\item either the $xy$-equilibrium state is {\em globally} stable,
\item or there is a {\em globally} stable limit cycle.
\end{enumerate}
By combining with the local stability result above we conclude
that the existence of a globally stable limit cycle occurs if and
only if the predator $y$ is efficient. Although this result gives
a complete qualitative description on the long-term dynamics of
the system, nothing can be said about the temporal outbreaks and
crashes. In addition, the method, as well as the closely related
Poincar\'e-Bendixson Theorem, is obviously 2-dimensional. We are
yet to see any extension of these methods to higher dimensional
webs and chains.

The third method is the phase plane analysis, more precisely, the
method of vector field analysis using the nullclines. The
$x$-nullcline, $\zeta dx/dt=xf(x,y,0)=0$, of
(\ref{eq2dPredatorPrey}) consists of the trivial branch, $x=0$,
and the nontrivial branch, a parabola $y=(1-x)(\beta_1+x)$ solved
from $f(x,y,0)=1 - x - \frac{y}{\beta_1 +x}=0$. The importance of
this method is immediately apparent---it tells where a population
decline turns into a recovery, and vise visa. More specifically,
if $dx/dt>0$ at a given set of populations $x,y$, then the prey
population $x(t)$ continues to increase at the time. Otherwise it
decreases if $dx/dt<0$. Therefore, the set of conditions at which
$dx/dt=0$ usually marks the transition between increase and
decrease in population. Exactly the same remarks apply to the
$y$-nullcline, $dy/dt=0$, which consists of the trivial branch
$y=0$ and the nontrivial branch $\frac{x}{\beta_1+x} - \delta_1=0$
or $x=x_{\rm ynl}:=\frac{\beta_1\delta_1}{1-\delta_1}$. Also
apparent is that the intersections of both variables nullclines
give rise to equilibrium points of the systems.

Features important for future analysis are: both $x$-branches
intersect at $(x_{\rm xtr},y_{\rm xtr})=(0,\beta_1)$, which is a
{\em transcritical point}, the maximum point of the parabola is
$(x_{\rm xfd},y_{\rm
xfd})=(\frac{1-\beta_1}{2},\frac{(1+\beta_1)^2}{4})$, which is a
{\em fold} point of the parabola $f(x,y,0)=0$, $dx/dt>0$ for
points below the parabola and $dx/dt<0$ for points above the
parabola. (More information is forthcoming on both transcritical
and fold points.) The $y$-transcritical point is $y=0,x=x_{\rm
ynl}$. Also, $dy/dt>0$ for points right of the nontrivial
$y$-nullcline, and $dy/dt<0$ for points left of it. Therefore a
phase portrait is given in Fig.\ref{figintrosing1}(a). Because of
an apparent vagueness of this equation's vector field, the method
does not always tell the stability of all the equilibrium points,
nor the existence of limit cycles. In addition, although the
nullcines separate population declines from recovery, they do not
foretell the magnitudes by which these events take place, which
often are of great practical importance because population
outbreaks and crashes are exclusively spoken of in quantitative
terms of magnitude. The method's qualitative nature is also its
major drawback.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}%[ht]
{\centerline{\includegraphics[width=5cm,height=4.5cm]{figintrosing1a.eps}}}

{\centerline{(a)}}

\smallskip

\parbox[b]{5cm}
{\centerline{\includegraphics[width=5cm,height=4.5cm]{figintrosing1b.eps}}}
\hskip 1cm
\parbox[b]{5cm}
{\centerline{\includegraphics[width=5cm,height=4.5cm]{figintrosing1c.eps}}}

\centerline{(b) \hskip 6cm  (c)}

\caption{(a) A typical vector field of the $xy$-system. The
$x$-component of the vector field points towards the solid
branches of the $x$-nullcline and away from the dashed ones.
Similar convention is used for branches of the $y$-nullcline.
Species $Y$ is weak iff $x_{\rm xfd}<x_{\rm ynl}$. (b) A typical
case of weak $y$ and the stability of the equilibrium point by
singular orbit analysis. The green curve is a relaxed orbit for
$0<\zeta\ll 1$. See text for the derivation of $y_{xpd}$ on the
phenomenon of Pontryagin's delay of loss of stability at which a
boom in $x$ population occurs although the recover starts after
its crossing the dashed part of the parabolic $x$-nullcline near
$y_{\rm xtr}$. (c) A typical case of efficient $y$, the existence
of a singular limit cycle and its relaxed cycle for $0<\zeta\ll
1$. \label{figintrosing1}}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This is where the singular perturbation method comes in to play
under the time diversification condition: $0<\zeta\ll 1$. The
method builds onto the last method of phase plane analysis where
the latter lacks: stability of equilibrium points, existence of
limit cycles, population booms and busts. All these questions can
be answered for the parameter range $0<\zeta\ll 1$.

First the order of magnitude in the course of evolution is
explicitly expressed in the dimensionless form
(\ref{eq2dPredatorPrey}). For small $0<\zeta\ll 1$, population $x$
changes at a fast order of $O(1/\zeta)$ if it is not near
equilibrium $xf(x,y,0)\approx 0$) already, comparing to an
ordinary order of a constant magnitude $O(1)$ for variable $y$.
All non-equilibrium solutions are quickly attracted to a small
neighborhood of the {\em stable} branches of the $x$-nullcline:
$x=0,y> y_{\rm xtr}$ and $f(x,y,0)=0, x> x_{\rm xfd}$. Once they
are there, the large magnitude $O(1/\zeta)$ is neutralized by
$x$'s being near the nullcline state $xf(x,y,0)\approx 0$. Then
the $y$ dynamics, which is negligible when $x$ is away from its
nullcline state, cannot be neglected further. The subsequent
development now evolves according to the ordinary time scale of
variable $y$. There are clearly two phases: the fast development
in $x$ followed by the slow evolution in $y$. In ecological terms,
if $(x,y)$ lies above the parabola $f(x,y,0)=0$, then $x$
undergoes a sudden decline or population bust or crash during the
$x$-fast phase. Otherwise, if $(x,y)$ lies below the parabola,
then population outbreak or boom occurs for $x$. In contrast,
during the $y$-slow phase the $y$ population undergoes a slow
decline if it follows from a bust in $x$, and a slow recovery if
it follows from a boom in $x$. These fast and slow dynamics can be
captured both qualitatively and quantitatively at the limit
$\zeta=0$ to equation (\ref{eq2dPredatorPrey}).

More precisely, setting $\zeta=0$ in Eq.(\ref{eq2dPredatorPrey})
results in
\begin{equation}\label{eq2dPredatorPreyslow}
%\begin{array}{c}
0= x\left(1 - x - \frac{y}{\beta_1 +x}\right), \quad \frac{dy}{dt}
= y \left(\frac{x}{\beta_1+x} - \delta_1\right).
%\end{array}
\end{equation}
It is a system of algebraic and differential equations. It can
also be viewed as a differential equation on the $x$-nullcline
manifolds: $x=0$ and $f(x,y,0)=0$. It captures the dynamics of the
slow phase in the sense that if $(x_\zeta,y_\zeta)(t)$ is a
segment of a solution of Eq.(\ref{eq2dPredatorPrey}) for $\zeta>0$
that is near the $x$-nullcline during the $y$-slow phase, then the
limit $(x_0,y_0)(t)=\lim_{\zeta\to 0}(x_\zeta,y_\zeta)(t)$ must
satisfy the limiting equation (\ref{eq2dPredatorPreyslow}). For
this reason, we call (\ref{eq2dPredatorPreyslow}) the {\em slow
subsystem} of Eq.(\ref{eq2dPredatorPrey}) and view
Eq.(\ref{eq2dPredatorPrey}) in the slow time scale of variable
$y$, which in fact was originally scaled against the maximum
reproductive rate of species $Y$.

Similarly, if we change the time variable of
Eq.(\ref{eq2dPredatorPrey}) from $t$ to $\tau=t/\zeta$, then the
equation becomes
\begin{equation}\label{eq2dPredatorPreyfast}
%\begin{array}{c}
\frac{dx}{d\tau}= x\left(1 - x - \frac{y}{\beta_1 +x}\right),
\quad \frac{dy}{d\tau} =\zeta y \left(\frac{x}{\beta_1+x} -
\delta_1\right).
%\end{array}
\end{equation}
The time scale is set according to that of variable $x$. Setting
$\zeta=0$ in the equation above gives rise to
\begin{equation}\label{eq2dPredatorPreyfastsub}
%\begin{array}{c}
\frac{dx}{d\tau}= x\left(1 - x - \frac{y}{\beta_1 +x}\right),
\quad \frac{dy}{d\tau} =0.
%\end{array}
\end{equation}
It is a 1-dimensional system in the fast variable $x$ with $y$
frozen as a parameter. Again, the equation captures the dynamics
during the $x$-fast phase in the same sense as equation
(\ref{eq2dPredatorPreyslow}) does to the $y$-slow phase: the limit
$(x_0,y_0)(\tau)=\lim_{\zeta\to 0}(x_\zeta,y_\zeta)(\tau)$ must
satisfy equation (\ref{eq2dPredatorPreyfastsub}) if
$(x_\zeta,y_\zeta)(\tau)$ is a solution segment of
Eq.(\ref{eq2dPredatorPreyfast}) during the $x$-fast phase. For the
obvious reason, we call Eq.(\ref{eq2dPredatorPreyfastsub}) the
{\em fast subsystem} of Eq.(\ref{eq2dPredatorPrey}) and view
Eq.(\ref{eq2dPredatorPreyfast}) in the fast reproductive scale of
the prey $x$.

In terms of terminology, solutions as well as orbits of
Eq.(\ref{eq2dPredatorPreyslow}) are described by adjective {\em
slow} whereas those of Eq.(\ref{eq2dPredatorPreyfastsub}) by {\em
fast}. These orbits are also referred to as singular orbits as
well as their natural concatenations. By natural concatenation it
means the following. A fast orbit of
Eq.(\ref{eq2dPredatorPreyfastsub}) must asymptotically approach a
points on the $x$-nullcline $x=0$ or $f(x,y,0)=0$. From that
point, a slow orbit of Eq.(\ref{eq2dPredatorPreyslow}) develops on
the $x$-nullcline. The union of these two orbits oriented in the
common sense of time is one case of the natural concatenation of
fast and slow singular orbits. In practical terms, it captures the
transitional events from a population boom in $x$ to its slow
decline and $y$'s slow rise, or a bust in $x$ to its slow recovery
and $y$'s slow decline. See Fig.\ref{figintrosing1}(b,c) for
illustrations of concatenated singular orbits and their perturbed
orbits.

The remaining case of natural concatenation is associated with
mechanisms by which singular orbits jump away, rather than into,
nullclines. In practical terms, it deals with sudden temporal
transitions of population: a slow decline and a slow recovery
respectively in population $y$ and $x$ followed by a sudden
outbreak in $x$, or a slow build up and a slow decline
respectively in $y$ and in $x$ followed by a sudden crash in $x$.
Because of their theoretical importance to the singular orbit
analysis both qualitatively and quantitatively, and because of
their practical importance to sharp temporal shifts in population,
we give a detailed illustration to each scenario below.

The first case takes place near the transcritical point $(0,y_{\rm
xtr})=(0,\beta_1)$ at which two branches of the $x$-nullcline,
$x=0$ and $f(x,y,0)=0$, intersect. The purpose of the analysis
below to understand how outbreak in $x$ takes place in relation to
a beginning amount of predator present.

Let $x=a, 0<a<x_{\rm xfd}=\frac{1-\beta_1}{2}$ be any line
sufficiently near the $y$-axis. Let $(a,y_1)$ be any initial point
that is on the line and above the parabola. Let
$x_\zeta(t),y_\zeta(t)$ be the corresponding solution of the
perturbed system (\ref{eq2dPredatorPrey}) with $0<\zeta\ll 1$ and
with the initial point. The solution must move down because $\dot
y<0$ for it is to the left of the non-trivial $y$-nullcline. It
moves leftwards above the unstable parabola $x$-nullcline until it
crosses the parabola at a point $(x_p,y_p), 0<x_p<a,\beta_1=y_{\rm
xtr}<y_p$ at which the solution curve is vertical. It then moves
down but rightwards since now the $x$ is in the recovery mode
$\dot x>0$. A finite time $T_\zeta>0$ later it intersects the
cross section line $x=a$ at a point denoted by $(a,y_2(\zeta))$.
In general the second intercept $y_2(\zeta)$ depends on $\zeta$.
Now integrating along the solution curve we have
\begin{equation*}
\begin{split}
0&=\zeta\ln x|_a^a=\int_0^{T_\zeta}\frac{\zeta\dot
x_\zeta}{x_\zeta}dt=\int_{y_1}^{y_2(\zeta)}\frac{\zeta\dot
x_\zeta}{x_\zeta\dot
y_\zeta}dy_\zeta\\
&=\int_{y_1}^{y_2(\zeta)}\frac{f(x_\zeta,y_\zeta)}{y_\zeta
g(x_\zeta,y_\zeta)}dy_\zeta.
\end{split}
\end{equation*}
Taking limit $\zeta\to 0$ on both sides of the equation above,
using the facts that $x_\zeta\to 0, x_p\to 0, y_p\to y_{\rm
xtr}=\beta_1$ and the notations that $y_1=y, y_2(\zeta)\to y_{\rm
xpd}, y_\zeta\to s$, we obtain the integral equation
\[
\int_{y_{\rm xpd}}^y \frac{f(0, s, 0)}{s g(0,s,0)} ds = 0
\]
for variables $y$ and $y_{\rm xpd}$. In practical terms, $y$
approximates the value followed by a collapse in $x$ starting at
the initial $(a,y)$ and $y_{\rm xpd}$ approximates the value that
immediately precedes a outbreak in $x$. Because the integrant
$f(0,s,0)/sg(0,s,0)$ has opposite signs for $s>y_{\rm xtr}$ and
for $s<y_{\rm xtr}$, we immediately conclude the following
qualitative fact: the larger the value $y$ after the $x$-collapse
the smaller the value $y_{\rm xpd}$ before the $x$-outbreak, as
well as its quantification: $y_{\rm xpd}$ is the value so that the
two areas bounded by the integrant graph are equal: $\int_{y_{\rm
xtr}}^y \frac{f(0, s, 0)}{s g(0,s,0)} ds = \int_{y_{\rm
xpd}}^{y_{\rm xtr}} -\frac{f(0, s, 0)}{s g(0,s,0)} ds$. Notice
that the branch $x=0,y<y_{\rm xtr}$ is effectively unstable for
the fast $x$-subsystem Eq.(\ref{eq2dPredatorPreyfastsub}).

The above phenomenon that singular orbits develop beyond a
transcritical point, continue along the unstable part of a
nullcline before jumping away is called {\em Pontryagin's delay of
loss of stability}. The transcritical point is referred to as a
{\em transcritical turning point}. The case illustrated is for the
type of transcritical points at which the fast variable goes
through a phase of crash-recovery-outbreak. In other cases of
generalization, they may be responsible for a reversal phase for
the fast variable. However, all the known PDLS cases of our
$XYZW$-model are of the crash-recovery-outbreak type described
above.

The last remaining case is associated with a slow build-up and a
slow decline in $y$ and $x$ respectively followed by a sudden
crash in $x$. For Eq.(\ref{eq2dPredatorPreyslow}) it happens only
when the $y$ species is efficient: $x_{\rm ynl}<x_{\rm xfd}$ as in
Fig.\ref{figintrosing1}(c) and near the fold point of the
parabolic $x$-nullcline $f(x,y,0)=0$. In this case we follows a
similar set-up and steps as for the transcritical turning point
above with a few modifications. First let $x=a$ with $x_{\rm
ynl}<a<x_{\rm xfd}$ instead and let $(a,y_1)$ be any initial point
from $x=a$ so that $y_1<y_{\rm xfd}=(1-x_{\rm xfd})(\beta_1+x_{\rm
xfd})$ instead. Because it is below the parabola, the perturbed
solution $(x_\zeta,y_\zeta)(t)$ moves right-up, crosses the stable
parabola $x$-nullcline vertical, moves left-up and hits the
cross-section line $x=a$ again at a point denoted by
$(a,y_2(\zeta))$. Use the exact orbit-integral-to-$\zeta$-limit
argument above to get a similar integral equation
\[
\int_{y_1}^{y_2(0)} \frac{f(x_0, y_0, 0)}{y_0 g(x_0, y_0,0)} dy_0
= 0,
\]
with $\lim_{\zeta\to 0}(x_\zeta, y_\zeta)=(x_0,y_0)$ being the
$y$-slow singular orbit on the parabola $f(x_0, y_0,0)=0$. Since
$y_2(\zeta)>y_{\rm xfd}$ for $\zeta>0$, we must have $y_2(0)\ge
y_{\rm xfd}$. The equation above does not hold if $y_2(0)> y_{\rm
xfd}$ because the integrant would be strictly negative for
integration interval $y_{\rm xfd}<y_0<y_2(0)$. Hence we conclude
that $y_2(0)=y_{\rm xfd}$, independent of the any initial point
$(a,y_1)$ with $y_1<y_{\rm xfd}$. In practical terms, the
collapses in $x$'s population occurs invariably whenever $y$'s
population reaches the fold value $y_{\rm xfd}$.

With all the necessary ingredients in place for the method of
singular orbit analysis we now complete our showcase for the
predator-prey model Eq.(\ref{eq2dPredatorPrey}). For the $y$-weak
case, $x_{\rm xfd}<x_{\rm ynl}$, all non-equilibrium singular
orbits converge to the co-existing steady state as shown in
Fig.\ref{figintrosing1}(b). For the $y$-efficient case, $x_{\rm
xfd}>x_{\rm ynl}$, all non-equilibrium singular orbits converge to
the limiting singular cycle $ABCD$ as shown in
Fig.\ref{figintrosing1}(c).

The importance of this analysis at the limit $\zeta=0$ lies in the
fact that all the singular asymptotic structures will persist for
small (often moderate $\zeta$ in practice) $0<\zeta\ll 1$! The
theory of persistence has been well-developed, c.f.
\cite{Pon57,OMa74,Fen79,Eck83,Bon87,Sch85,Mis94,Den99}, which
enables us to focus our attention primarily on singular orbit
structures in applications if the main purpose is to understand
the underlining dynamics. For this reason, we will only make
sparse comments on persistence questions throughout the paper.

Comparing the singular orbit analysis to the other methods
described above, it is several magnitude superior! Its only
drawback seems to be its conspicuously lack of treatment for large
$\zeta\gg 1$. A closer examination shows however the drawback is
not as significant as it first seems. First, the practical
interpretation of large $\zeta>1$ implies that the predator
out-reproduces its prey. Such rate-reversal systems are less
common than their rate non-reversal counterparts. Second, if the
predator has the potential to multiply faster than its prey, then
Verhulst's logistic growth assumption must be incorporated into
predator's model by assuming that $Y$'s per-capita death rate is
density dependent: $d+d_0Y$. Upon non-dimensionalizing, the $y$
takes this form $dy/dt=y(x/(\beta_1+x)-\delta_1-\delta_0y)$. For
large $\zeta\gg 1$, the resulting model is again a singular
perturbed system for which $y$ is fast and $x$ is slow. The same
singular orbit analysis can be applied again. Without Verhulst's
assumption, application of the singular orbit analysis will fail.
More specifically, since the nontrivial $y$-nullcline $x=x_{\rm
ynl}$ would be a vertical line parallel to all $y$-fast singular
orbits, at the limit $\zeta=\infty$, all fast $y$-orbits for
$x>x_{\rm ynl}$ would fly to infinity without bound, which would
leads to an ecological absurdity that a fixed amount of prey
sustains any ever increasing amount of predator. The failure does
not come from the method rather than the non-realistic assumption
of non-logistic growth on the fast reproducing predator $y$. We
will not treat the rate-reversal models any further in this paper.

We conclude this section by pointing out that a summary of the
singular orbit analysis methodology, in particular, analysis of
sharp temporal transitions involving multiple fold and
transcritical turning points, is given in Appendix C for future
reference as we include more and more species to our exposition in
subsequence exposition.
