\section{Appendix A}

In general, let $\{X_1,X_2,\dots,X_n\}$ be a community of species
ordered in chains from bottom to top. It is said to be a {\em
direct} chain if $X_i$ is the prey of $X_{i+1}$ for all $1\le i\le
n-1$. We call $X_k$ a {\em chain predator} of any $X_i$ for $i<k$
and simply a predator if $i=k-1$ in a direct chain. We call $X_i$
a {\em chain prey} of $X_k$ for $i<k$ and simply a prey for
$k=i+1$ in a direct chain. The community is said to form a {\em
web} if any $X_k, k\ge 1$ is either a chain predator or a chain
prey and there is an $i$ such that $X_i$ is not a prey of
$X_{i+1}$. Two species $X_i, X_j$ are said to be {\em connected}
in a web if either one is a chain predator of the other in a
direct subchain or both share a common chain predator or chain
prey. A web is said to be {\em connected} if any two species in it
are connected and it is always meant to be connected unless
otherwise noted whenever a web is considered. A direct chain is a
special case of connected web.

A time continuous model for a community usually take the form:
$\dot X_i=X_if_i(X)$, $X=(X_1,X_2,\dots,X_n), i=1,2,\dots, n$. Any
hyperplane $X_{i_1}=X_{i_2}=\dots =X_{i_j}=0$ is invariant. In
other words, any collection of $X_{n_1},X_{n_1},\dots,X_{n_k}$
with $X_i=0,i\ne n_j, j=1,2,\dots,k$ forms a subsystem. $X_{n_k}$
is said to be $X_{n_1}X_{n_1}\cdots X_{n_k}$-{\em weak} if
$\{X_{n_1},X_{n_1},\dots, X_{n_k}\}$ is a direct chain and it has
a stable $X_{n_1}X_{n_1}\cdots X_{n_k}$ equilibrium point with
nontrivial values $X_{n_i}>0$. Let $X_{n_{k+1}}$ be a competitor
of $X_{n_{k}}$ such that both $X_{n_k}$ and $X_{n_{k+1}}$ share
the same subchain, i.e. there are species $X_{n_1},X_{n_2},\dots,
X_{n_{k-1}}$ so that both $X_{n_1}X_{n_2}\dots
X_{n_{k-1}}X_{n_{k}}$ and $X_{n_1}X_{n_2}\dots
X_{n_{k-1}}X_{n_{k+1}}$ are direct chains. Then $X_{n_{k+1}}$ is
said to be $X_{n_1}X_{n_2}\cdots X_{n_{k-1}}X_{n_{k}}$-{\em
competitive} if any $X_{n_1}X_{n_2}\cdots X_{n_{k-1}}X_{n_{k}}$
attractor is asymptotically unstable with respect to the expanded
$X_{n_1}X_{n_2}\cdots X_{n_{k-1}}X_{n_{k}}X_{n_{k+1}}$-subsystem.
By definition it is equivalent to
\[
\lim_{T\to\infty}\frac{1}{T}{\int_0^Tf_{n_{k+1}}(X_a(t))dt}>0
\]
with $X_a(t)$ denoting any dense solution on the attractor. This
criterion is derived from the variational equation of the
$X_{n_{k+1}}$-equation
\[
\dot V=f_{n_{k+1}}(X_a(t))V
\]
with $V=DX_{n_{k+1}}$ being the variation of the
$X_{n_{k+1}}$-variable because $X_{n_{k+1}}=0$ on any
$X_{n_1}X_{n_2}\cdots X_{n_{k-1}}X_{n_{k}}$ attractor. Solving $V$
we have
\[
V(t)=V_0\exp(\int_0^tf_{n_{k+1}}(X_a(s))ds):=V_0\exp(\lambda_t t)
\]
where $\lambda_T=\frac{1}{T}\int_0^Tf_{n_{k+1}}(X_a(s))ds$. Hence
$\lim_{T\to\infty}\lambda_T>0$ implies the asymptotic instability
of the $X_{n_1}X_{n_2}\cdots X_{n_{k-1}}X_{n_{k}}$ attractor.

If the attractor is an equilibrium point, this condition reduces
to $f_{n_{k+1}}(X_a(t))\equiv f_{n_{k+1}}(X_a(0)) >0$. If it is a
$T_c$ periodic cycle, it is reduced to
\[
\frac{1}{T_c}{\int_0^{T_c}f_{n_{k+1}}(X_a(t))dt}>0.
\]
If it is an irreducible $X_{n_1}X_{n_2}\cdots
X_{n_{k-1}}X_{n_{k}}$ chaotic attractor, then it is the
$X_{n_{k+1}}$-directional positive Lyapunov exponent in the
expanded $X_{n_1}X_{n_2}\cdots X_{n_{k}}X_{n_{k+1}}$ system.

\section{Appendix B}

The technique is also based on a geometric aspect of the
problem---the nullcline analysis presented in Sec.4. Let
\[
\dot u=\frac{1}{\zeta}F(u,v),\quad \dot v=G(u,v)
\]
denote the linearization of Eq.(\ref{eq2dPredatorPrey}) after we
have translated the corresponding equilibrium point $(x_e,y_e)$ to
the origin $(0,0)$. Since it is the linearization of $xf(x,y,0)$,
$F(u,v)$ must take this form $F(u,v)=au+bv$ for some $a,b$. To
determine the signs of $a$ and $b$, we take a closer look at the
circumstances from which it arises. First the nullcline $F(u,v)=0$
is precisely the tangent line to the nullcline $xf(x,y,0)=0$. By
inspecting the nullcline $xf(x,y,0)=0$ as shown in
Fig.\ref{figintrosing1}(a), we conclude that the line
$F(u,v)=au+bv=0$ has a negative slope iff $y$ is weak: $x_{\rm
xfd}<x_{\rm ynl}$. Hence, $a, b$ must have the same sign. Second,
since $xf(x,y,0)<0$ iff the point $(x,y)$ lies above the nullcline
$f(x,y,0)=0$, $F(u,v)$ preserves the same property: $au+bv<0$ iff
$(u,v)$ lies above the line $au+bv=0$, which implies that
$a<0,b<0$. The exactly same argument can be applied to conclude
that $G(u,v)=cu$ for some $c>0$. Now the eigenvalues of the
lineariation can be explicitly derived as
\[
\lambda_{1,2}=\frac{1}{2}\left[\frac{a}{\zeta}\pm
\sqrt{\frac{a^2}{\zeta^2}+\frac{4bc}{\zeta}}\right].
\]
We can then conclude the following:
\begin{enumerate}
\item The eigenvalues have negative real parts iff $y$ is weak, i.e.
$a<0$ for $a<0,b<0,c>0$.
\item The equilibrium point undergoes a Hopf bifurcation at the
fold point when $a=0$.
\end{enumerate}
The same amount of qualitative information without actually
finding the equilibrium point and the linearization explicitly.



\section{Appendix C}

The key idea for singular orbit analysis is to break down, if
possible, a system into various time scales so that we only have
to deal with one or fewer variables at a time. The scaled system
(\ref{ScaledSystem}) indeed evolves at various different time
scales, provided $0 < \zeta, \epsilon_1, \epsilon_2 \ll 1$. Such
systems are called {\em singularly perturbed systems} and the
small time rate parameters are called the {\em singular
parameters}.

By expressing the derivatives explicitly in (\ref{ScaledSystem})
as follows
%\begin{equation*}
%\begin{array}{cc}
\[
\dot x = \frac{1}{\zeta} x f(x,y,z), \quad \dot y = y g(x,w),
\quad \dot z = \epsilon_1 z h(x), \quad \dot w = \epsilon_2 w
k(y),
\]
%\end{array}
%\end{equation*}
we see that the derivative $\dot x$ is extremely large due to the
$\frac{1}{\zeta}$ term. Thus all solutions evolve very quickly to
the nearest attracting $x$-nullcline surface, where $\dot x = 0$,
with other variables varying little near their initial values. (If
there is no attracting $x$-nullcline surface in that $x$-flow
direction, which is not the case for our system, the solution
diverges to infinity very quickly.) At the limit, $\zeta=0$, the
other variables are frozen at their initial values and the
$x$-initial instantaneously jumps in the direction of $x$ to the
attracting $x$-nullcline surface, which consists of stable
equilibrium points for the $x$-equation with other variables fixed
as parameters.

Once the solution arrives on an attracting branch of the
$x$-nullcline, the next fastest part of the system takes over.
With the assumption that $0<\epsilon_2, \epsilon_1 \ll 1$ it is
the $y$ variable. Holding $z$ and $w$ as constant, the solution
may proceed to the nearest attracting $y$-nullcline (following the
1-directional $y$-flow) or may leave the $x$-nullcline via fold
points or transcritical points. If it goes to the $y$-nullcline,
we then look for the next fast variable which may be either $z$ or
$w$ or both depending on the relative sizes of $\epsilon_1$ and
$\epsilon_2$, and repeat the process during which the singular
solution may bounce among various nullclines or attractors before
settling down on an equilibrium point, a limit cycle, or a chaotic
attractor. Jumping off a nullcline occurs usually at fold or
transcritical turning points. The fast variable at a given time
scale can be more than 1 variable. For instance, since $y,z$ are
competitors which occupy the same trophic level above $x$, they
may reproduce at comparable rates. In such a case $(y,z)$ is the
next joint fast variable after $x$. As a consequence the singular
orbit may be attracted to not only the combined
$yz$-equilibrium/nullcline states but also other types of
$yz$-attractors such as invariable tori. The singular orbit
analysis still applies in that the next corresponding singular
orbit stops on such attractors before further developing.
Presumably the dynamical behaviors vary with particular models.

The short-term temporal transitions illustrated in Sec.4 are the
simplest and most elementary kinds of
predator-decline-prey-recovery-to-prey-outbreak and
predator-recovery-prey-decline-to-prey-collapse that one
encounters in ecological models. Mathematically, the derivation
involves only a simple substitution from the integration variable
$x_\zeta$ to $y_\zeta$ once since $y_\zeta(t)$ is a monotone
function of $t$ from $y_1$ to $y_2$. In generalizations to this
argument for some more elaborate outbreak-collapse population
dynamics, the fast $x$-nullcline may be surfaces or hyper-surfaces
and the slow $y$-flow may be multi-dimensional having once again
more than 1 time scales with fold and transcritical turning points
of its own. As a result the substitution in $y_\zeta$ may require
more-than-one monotone substitutions and the integral equation may
contain multiple pieces of integrals. However complicated the
models may be, the technique of orbit-integral-to-$\zeta$-limit
will always be an attractive tool to try, which always worked
whenever it is needed throughout the presentation. Many authors
have effectively used this method of which the elementary
derivation given above follows the original idea of \cite{Pon57},
see also \cite{Sch85,Mis94}.

The method of singular perturbations is to categorize asymptotic
singular orbits, which attracts other singular orbits near by, and
to demonstrate that the attracting singular orbits persist for
small but non-vanishing singular parameters. The importance of
studying singular attractors lies in the fact that almost
inevitably such singular attractors persist for the perturbed full
systems.


\section{Appendix D}

\begin{pro}\label{propyfold} Under conditions:
\[
\max\{\frac{1-\beta_i}{2}\}<\min\{\frac{\beta_i\delta_i}{1-\delta_i}\}\hbox{
and }
\beta_3<\frac{(\beta_1+1)^3}{\beta_1}\left(\frac{1}{\beta_1+1}-\delta_1\right)
%\beta_3<\frac{1-\delta_1}{\beta_1}\left(1-\frac{\beta_1\delta_1}{1-\delta_1}\right)(\beta_1+1)^2
\]
there is a unique $y$-fold curve $p_y(y,z)=0$ in $\Delta$ from
$z=0$ to $y=0$, and $p_y(0,0)>0,\ p_y(y,z)<0$ for $q(y,z)=x_{\rm
ynl}$.
\end{pro}

\begin{proof} The idea is to show that along any radius line
$y=s,z=ms$ in the first quadrant $m\ge 0,s\ge 0$ from the origin
$(0,0), s=0$ to the $\Delta$ boundary $q(s,ms)=x_{\rm ynl}$, there
is a unique zero for the function $w_y=p_y(s,ms)$.

With $x=q(y,z)\in [x_{\rm ynl},1]$ below, the $y$-partial
derivative of $w$ is
\begin{equation*}
\begin{split}
w_y=& p_y(y,z)=\frac{\beta_1}{(\beta_1+x)^2}q_y(\beta_3+y)+
\left(\frac{x}{\beta_1+x}-\delta_1\right)\\
=&\frac{\beta_1}{(\beta_1+x)^2}\left(-\frac{f_y}{f_x}\right)(\beta_3+y)+
\left(\frac{x}{\beta_1+x}-\delta_1\right)\\
=&\left(-\frac{f_y(\beta_1+x)}{f_x}\right)\left[\frac{\beta_1}{(\beta_1+x)^3}(\beta_3+y)+
\left(-\frac{f_x}{f_y(\beta_1+x)}\right)\left(\frac{x}{\beta_1+x}-\delta_1\right)\right]\\
=&\frac{1}{f_x}\left[\frac{\beta_1}{(\beta_1+x)^3}(\beta_3+y)+
f_x\left(\frac{x}{\beta_1+x}-\delta_1\right)\right]\hbox{ (since
$f_y=-\frac{1}{\beta_1+x}$)}
\end{split}
\end{equation*}
Since $f_x=-1+\frac{y}{(\beta_1+x)^2}+\frac{z}{(\beta_3+x)^2}=-1$
at $x=1,y=z=0, s=0$, we have immediately
\begin{equation*}
\begin{split}
w_y|_{\{x=1,y=z=0\}}&=-\left[\frac{\beta_1\beta_3}{(\beta_1+1)^3}-
\left(\frac{1}{\beta_1+1}-\delta_1\right)\right]\\
&=-\frac{\beta_1}{(\beta_1+1)^3}\left[\beta_3-
\frac{(\beta_1+1)^3}{\beta_1}\left(\frac{1}{\beta_1+1}-\delta_1\right)\right]>0
\end{split}
\end{equation*}
by the assumption on $\beta_3$. Also, at $x=q(s,ms)=x_{\rm
ynl}=\beta_1\delta_1/(1-\delta_1)$,
\[
w_y|_{\{x=x_{\rm
ynl}\}}=\left(\frac{1}{f_x}\right)\frac{\beta_1}{(\beta_1+x)^3}(\beta_3+y)<0
\]
because $f_x<0$ in ${\mathcal S}$ by definition. Therefore there
must a zero of $w_y=p(s,ms)$ between $s=0$ and $q(s,ms)=x_{\rm
ynl}$. So it is only left to show that the zero is unique.

To this end, we make a simple change of variable $x=q(s,ms)$ so
that $x$ runs from $1$ to $x_{\rm ynl}$ with increasing $s$ from
$s=0$. Since $x_s=q_y+mq_z<0$ because both $q_y, q_z<0$ and $m\ge
0$, this is a valid substitution. Moreover, the inverse $s=s(x)$
is well defined, and is monotone decreasing: $s_x<0$. Denote
\[
w_y=\frac{1}{f_x}\left[\frac{\beta_1}{(\beta_1+x)^3}(\beta_3+y)+
f_x\left(\frac{x}{\beta_1+x}-\delta_1\right)\right]:=\frac{1}{f_x}Q(x).
\]
Then we know $Q(1)<0,Q(x_{\rm ynl})>0$, and we only need to show
$Q'(x)<0$ for $x\in[x_{\rm ynl},1]$. Writing $Q$ out
\[
Q(x)=\frac{\beta_1}{(\beta_1+x)^3}(\beta_3+s(x))+
\left(-1+\frac{s(x)}{(\beta_1+x)^2}+\frac{ms(x)}{(\beta_3+x)^2}\right)\left(\frac{x}{\beta_1+x}-\delta_1\right)
\]
we see that $Q'(x)<0$. In fact, the first term is obviously
decreasing in $x$. The second term is the product of two factors:
$a(x)b(x):=f_x(\frac{x}{\beta_1+x}-\delta_1)$. By product rule, we
have $a'(x)b(x)+a(x)b'(x)$. Since $a=f_x$ is negative, decreasing,
and $b$ is positive, increasing, we see clearly that it is
decreasing for the product. This completes the proof.
\end{proof}
