\section{Main Results, Singular Perturbation Analysis}

The key question is what are the possible long term dynamical
outcomes of two weak competitors when a top-predator is introduced
over the stronger one? We will use the $xyzw$-model
(\ref{ScaledSystem}) to study the problem assuming the drastic
trophic time diversification ($0<\zeta, \epsilon_2\ll 1$ for the
chain $xyw$, and a slower reproducing competitor $z$ with
$0<\epsilon_1\ll 1$ and a comparable one with $\epsilon_1=O(1)$.
Since $x_{\rm ynl}>x_{\rm xfd}|_{\{z=0\}}$ and $x_{\rm znl}>x_{\rm
xfd}|_{\{y=0\}}$, and because we will study qualitative changes
when the relative competitiveness $E_{xy}=x_{\rm ynl}-x_{\rm
znl}=-E_{xz}$ exchanges hand, it is necessary to assume
\begin{equation}\label{WeakCompetitors}
\min\{x_{\rm ynl}, x_{\rm znl}\}>\max\{x_{\rm xfd}|_{\{z=0\}},
x_{\rm xfd}|_{\{y=0\}}\},
\end{equation}
throughout our discussion.

We will use the same singular orbit analysis as in the previous
section but adapt it to the 4-dimensional case. Hence, our
analysis contains detailed geometrical analysis on the nullcline
manifolds, their fold turning points, and their transcritical
turning points. We will do this in an order dictated by the time
scales of singular orbits: from fastest variable down to the
slowest. This will allow us to capture typical long term singular
orbit attractors. We will only choose to comment on the singular
attractors' persistence, which almost always is the case in spite
of conspicuous absence of some rigorous proofs.

\subsection{Nullclines And Frozen Dynamics}
We now discuss singular orbits that are on the trivial
$x$-nullcline solid $x=0$, with initial conditions of
$(0,y_0,z_0,w_0)$. When on the $x=0$ solid, $\dot y =
y(-\frac{w}{\beta_3 + y} - \delta_1) < 0,\dot z = \epsilon_1z(-
\delta_2) < 0$. So the $yz$-flow will decrease towards $(y,z)=0$.
Moreover, since the transcritical surface
(\ref{TranscriticalLine}) ($f(0,y,z)=1-y/\beta_1-z/\beta_2=0$) has
a negative slope, the $yz$-vector field is transversal to the
transcritical line, and as a result the $yz$ flow will always
cross it to the unstable side ($f(0,y,z)>0$) of the $x=0$ solid.
By Pontryagin's delay of loss of stability, the orbit will leave
the $x=0$ surface and develop in the $x$ direction to the
attracting nontrivial $x$-nullcline manifold in $\{f(x,y,z)=0\}$

Because solutions will always leave the $x=0$ solid for the
nontrivial $x$-nullcline, and stay there as to be seen shortly
under the weak predator assumption on $y, z$, we will just
consider Eq.(\ref{ScaledSystem}) constraint to the $x$-manifold
$f(x,y,z)=0$ except in a few cases where such behavior is noted.

The attracting branch of the nontrivial $x$-nullcline must satisfy
$f_x(x,y,z)<0$ by definition. Therefore it can be solved from
$f(x,y,z)=1-x-\frac{y}{\beta_1+x}-\frac{z}{\beta_2+x}=0$ by
implicity function theorem in terms of $x$ as a function of $y,z$
or explicitly as in Sec.4. Because it is independent of $w$, the
manifold is a solid or hyper-surface in the $xyzw$-space, parallel
to the $w$-axis, which we denote by ${\mathcal S}: x=q(y,z)$ for
points $(y,z)$ interior to the fold surface (\ref{FoldLine}).
Since $f_x<0$ on ${\mathcal S}$ by the definition for the stable
branch of the $x$-nullcline, and $f_y<0,f_z<0$, we have that
$q_y=-f_y/f_x<0$ and $q_z=-f_z/f_x<0$. Hence, $q$ decreases in
both $y$ and $z$ variables, with $q(0,0)=1$. When projected to the
$yzw$-space, ${\mathcal S}$ is a cylindrical solid in the first
$yzw$-octant, bounded by the coordinate planes and the fold
cylindrical surface (\ref{FoldLine}). In other words, it is the
domain of definition in the positive $yzw$-octant for the function
$q$. We will use the same notation $S$ for its $yzw$-projection
solid.

Eq.(\ref{ScaledSystem}) when restricted onto the $x$-slow manifold
solid ${\mathcal S}$ is the $yzw$-equations with $x=q(y,z)$. It is
again singularly perturbed by parameters $0<\epsilon_i\ll 1$ with
$y$ being the fast variable. Hence, we need to consider its
$y$-nullcline in turn for the next level of singular orbit
construction. The $zw$-plane ($y=0$) is always the trivial
$y$-nullcline. The intersection of the nontrivial $y$-nullcline
within ${\mathcal S}$ is determined by the system of equations
$f(x,y,z)=0,g(x,y,w)=0$, which is expressed as
\[
w=p(y,z):=\left(\frac{q(y,z)}{\beta_1+q(y,z)}-\delta_1\right)(\beta_3+y),
\]
with $(y,z)$ constraint to the region $\Delta:=\{
\frac{q(y,z)}{\beta_1+q(y,z)}-\delta_1\ge 0, y\ge 0,z\ge 0\}$,
equivalently $\Delta=\{q(y,z)\ge
\beta_1\delta_1/(1-\delta_1)=x_{\rm ynl}, y\ge 0,z\ge 0\}$. It is
a surface in the $x$-slow manifold solid ${\mathcal S}$. Since $q$
is decreasing in $y,z$, and the boundary curve $x=q(y,z)=x_{\rm
ynl}$ is a line, obvious from the defining equation
$f(x,y,z)=(1-x-y/(\beta_1+x)-z/(\beta_2+x)=0$, therefore the
domain of definition $\Delta$ for the $y$-nullcline $w=p(y,z)$ in
${\mathcal S}$ is a triangle bounded by the axises $y=0,z=0$ and
the line $q(y,z)=x_{\rm ynl}$, see Fig.\ref{figxyzw}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}%[ht]
\label{figxyzw}
\parbox[b]{5cm}
{\centerline{\includegraphics[width=5cm,height=5cm]{figxyzwview.eps}}}
\hskip 1cm
\parbox[b]{5cm}
{\centerline{\includegraphics[width=5cm,height=5cm]{figxyzwyzview.eps}}}

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

\caption{(a) Nullcline surfaces inside the $x$-nullcline solid
${\mathcal S}$ which is bounded by the $x$-fold plane parallel to
the $w$-axis. (b) The same surfaces in the $yz$-projected view.
Left to the $y$-fold curve, the $y$-nullcline surface is unstable,
and to its right, it is stable. The $y$-nullcline surface is
bounded by the co-ordinate axises and the green $y$-boundary
curve.}
\end{figure}

The $y$-transcritical curve is given by
$w=p(0,z)=(\frac{x}{\beta_1+x}-\delta_1)\beta_3$ with $x=q(0,z)$
which satisfies $z=(1-x)(\beta_2+x)$, with $z$ between $0$ and
$z=(1-x_{\rm ynl})(\beta_2+x_{\rm ynl})$. It is monotone
decreasing in $z$. In fact,
$w_z=p_z(0,z)=\frac{\beta_1\beta_3}{(\beta_1+x)^2}x_z$ whose sign
is the same as $1/x_z=z_x=2(\frac{1-\beta_2}{2}-x)<0$ because
$x\ge x_{\rm ynl}>x_{\rm xfd}|_{\{y=0\}}=(1-\beta_2)/2$ by our
assumption (\ref{WeakCompetitors}).

As we mentioned earlier in Sec.5 that the $xyw$-system with $z=0$
has a unique $y$-nullcline fold point on the stable $x$-manifold
$f=0$, and it is a special case of Proposition \ref{propyfold}.
More specifically, a fold turning point for the $y$-flow in the
solid ${\mathcal S}$ is determined by $g_y(q(y,z),y,w)=0,
g(q(y,z),y,w)=0$. In terms of $w=p(y,z)$ for the $y$-nullcline in
${\mathcal S}$: $g(q(y,z),y,w)=0$, $g_y(q(y,z),y,w)=0$ is
equivalent to $w_y=p_y(y,z)=-g_y(q(y,z),y,w)/g_w(q(y,z),y,w)=0$.
Also, the stable branch of the $y$-nullcline is determined by
$g_y(q(y,z),y,w)<0$, which is equivalent to $w_y=p_y(y,z)<0$
because $g_w<0$ always.

Proposition \ref{propyfold} proves that the $y$-nullcline surface
$w=p(y,z)$ in ${\mathcal S}$ indeed has a unique $y$-fold curve
$w_y=p_y(y,z)=0$ on $w=p(y,z)$ running from $\Delta$'s boundary
$z=0$, which is the above special case, to another boundary $y=0$.
Denote this fold curve by $(y,z,w)_{\rm yfd}$. This curve divides
the $y$-nullcline into the stable branch and the unstable branch.
The unstable part contains the $y$-transcritical point, in
particular $x=1,y=z=0, w=q(0,0)$, and the stable part contains the
boundary $q(y,z)=x_{\rm ynl}$, see Proposition \ref{propyfold}.

The nontrivial $z$-nullcline $h(x)=0$ in the $x$-slow manifold
solid $f(x,y,z)=0$ is given by $x=x_{\rm znl},
f(x,y,z)=1-x-y/(\beta_1+x)-z/(\beta_2+x)=0$. It is a line in the
$yz$-plane and a plane parallel to the $w$-axis in the solid
${\mathcal S}$. Important properties about it include the
following:
\begin{enumerate}
\item If $z$ is competitive ($E_{xy}=x_{\rm ynl}-x_{\rm znl}>0$), then the $z$-nullcline $x=x_{\rm znl}$ does not
intersect the $y$-slow manifold
$w=p(y,z)=(\frac{x}{\beta_1+x}-\delta_1)(\beta_3+y)$ since $w<0$
with $x=x_{\rm znl}$. On the $yz$-plane in ${\mathcal S}$, the
line $x=x_{\rm znl}$ lies between the $x$-fold curve and the
boundary line $x=q(y,z)=x_{\rm ynl}$ for for the $y$-slow manifold
$w=p(y,z)$. See Fig.\ref{figxyzw}.
\item If $z$ is not competitive ($E_{xy}=x_{\rm ynl}-x_{\rm znl}<0$), then the $z$-nullcline $x=x_{\rm znl}$ intersects the
$y$-slow manifold
$w=p(y,z)=(\frac{x}{\beta_1+x}-\delta_1)(\beta_3+y)$ in the solid
${\mathcal S}:f(x,y,z)=1-x-y/(\beta_1+x)-z/(\beta_2+x)=0$ along a
line
\[
w=p(y,z)=\left(\frac{x}{\beta_1+x}-\delta_1\right)(\beta_3+y),\quad
z=(\beta_2+x)\left(1-x-\frac{y}{\beta_1+x}\right)
\]
parameterized by $y\in [0,(1-x)(\beta_1+x)]$, with $x=x_{\rm znl}$
above. Along it $w$ increases while $z$ decreases with increase in
$y$. We denote this line by $(y,z,w)_{\rm znl}$.
\item The $z$-transcritical line is given by the intersection of the
trivial $z$-nullcline $z=0$ and $x=x_{\rm znl}$ in the solid
${\mathcal S}: f(x,y,z)=1-x-y/(\beta_1+x)-z/(\beta_2+x)=0$, which
is the line $y_{\rm ztr}=(1-x_{\rm znl})(\beta_1+x_{\rm znl})$.
Point $(y_{\rm ztr},0,w_{\rm ztr})=(y_{\rm ztr},0, p(y_{\rm
ztr},0))$ is the $z$-transcritical point on the $y$-slow manifold
$w=p(y,z)$. See Fig.\ref{figxyzw}.
\end{enumerate}

Finally, the nontrivial $w$-nullcline ($k(y)=0$) is given by the
hyper-plane
\[
y=y_{\rm wnl}=\frac{\beta_3\delta_3}{1-\delta_3}.
\]
It intersects the $y$-slow manifold $w=p(y,z)$ in ${\mathcal S}$
along a curve $w=p(y_{\rm wnl},z)$ with $z$ from $0$ to $z_{\rm
wnl} =(1-x_{\rm ynl}-y_{\rm wnl}/(\beta_1+x_{\rm
ynl}))(\beta_2+x_{\rm ynl})$, which is the intersection of
$y=y_{\rm wnl}$ with the boundary $w=0, q(y,z)=x_{\rm ynl}$. This
$w$-nullcline will intersect the $y$-fold curve if $w$ is
efficient as defined in Sec.6, or may not if $w$ is weak as shown
in Fig.\ref{figxyzw}. The point $(y_{\rm wnl},z_{\rm wnl},w)$ is
the $w$-transcritical point on the $y$-slow manifold $w=p(y,z)$.

As some immediate consequences to our geometric nullcline analysis
above, we summarize some results below before carrying out our
main singular perturbation analysis.

\subsubsection{Invariance of Solid ${\mathcal S}$} We demonstrate
here that under the weak competitor condition
(\ref{WeakCompetitors}), no singular solutions can escape from the
$x$-slow manifold solid ${\mathcal S}$. In fact, because of the
range restriction in $x$ for the $x$-fold in (\ref{FoldLine}), it
always lie below both the $y$-nullcline $x=x_{\rm ynl}$ when $w=0$
and the $z$-nullcline $x=x_{\rm znl}$. Therefore, both $y$ and $z$
decrease on the $x$-fold. With respect to the $x$-slow manifold
solid ${\mathcal S}$, the vector field on the $x$-fold surface
always points into the solid, therefore no singular orbits can
escape.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{figydmtxyz}
\centerline{\includegraphics[width=5in,height=4.8in]{figydmtxyz.eps}}
\caption{(a) Nullclines surfaces for the $yzw$-system in the
$x$-slow manifold solid ${\mathcal S}$. (b,c,d) Competing
$yz$-phase planes for frozen $w$ with $\epsilon_2=0$.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsubsection{Dominant $Z$ Dynamics} The long term dynamics for
the simplest case that $z$ is $xyz$-dominant ($E_{xy}>0$) does not
depends on the singular perturbation structure for the constraint
$yzw$-slow system on the solid ${\mathcal S}$. From Sec.??? we
have for $w_0=0$, $w(t)\equiv 0$ and $y(t)\to 0,z(t)\to z_{\rm
znl}|_{\{y=w=0\}}$. Because $yg(x,y,0)>yg(x,y,w)$, we have
\[
y|_{w_0=0}(t)\ge y|_{w_0>0}(t)\Longrightarrow y|_{w_0>0}(t)\to 0,
\]
when all other conditions are equal. This leads eventually to
$w(t)\to 0$. Hence, we have proved the following statement

\begin{pro}\label{zdmt} If $z$ is $xyz$-dominant,$E_{xy}>0$, then
$z$ is $xyzw$-dominant, i.e., $y(t)\to 0, w(t)\to 0$
\end{pro}

For the remainder of the paper, we always assume that $y$ is
$xyz$-dominant ($E_{xy}<0$), in addition to our principle
assumption that both $y,z$ are weak competitors
(\ref{WeakCompetitors}).

\subsubsection{$W$-Induced Competition Exclusion} Under the
condition that $y$ is $xyz$-dominant we know the $z$-nullcline
must cut through the $y$-slow manifold $w=q(y,z)$ in ${\mathcal
S}$ along a curve which is decreasing in $w$ with increasing $z$
as shown in Fig.\ref{figydmtxyz}(a). If we hold each $w$ constant
in Eq.(\ref{ScaleSystem}) with $x=p(y,z)$, the system is a
2-dimensional competing species system in $y,z$. The phase
portrait can be obtained by $w$-sectioning for the portrait of
Fig.\ref{figydmtxyz}. In fact, let $w=w_0>0$ be any fixed
constant. There are 3 typical cases for the competing
$yz$-dynamics. (1) For $0<w_0<w_{\rm znl}|_{\{y=0\}}$, the
$xy$-equilibrium point is globally attracting in the solid
${\mathcal S}$ as in Fig.\ref{figydmtxyz}(b). (2) For $w_{\rm
znl}|_{\{y=0\}}<w_0<w_{\rm znl}|_{\{z=0\}}$, both $xy$-equilibrium
point and $xz$-equilibrium point are locally stable with their
basins of attraction separated by a separatrix through the $xyz$
equilibrium point which is unstable, as in
Fig.\ref{figydmtxyz}(c). (3) For $w_{\rm znl}|_{\{z=0\}}<w_0$, $z$
dominates, all nontrivial solutions converge to the
$xz$-equilibrium point as in Fig.\ref{figydmtxyz}(d). The
remaining cases are not typical, they are bifurcations from one of
the three cases (1,2,3) to another. This demonstrates that
predation on $y$ from $w$ enhances the competitiveness of $y$'s
competitor. Competition co-existence is excluded with fixed amount
of $w$. One of the main results of this paper demonstrates
surprisingly that this principle is not to be extended in for a
$w$ dynamically evolving with the system.

\subsubsection{$Z$-Induced $W$-Weakness} If we hold $z$ constant
with frozen $\epsilon_1=0$, then the reduced $xyw$-chain
(\ref{ScaleSystem}) in the solid ${\mathcal S}$ is a predator-prey
system in $y, w$. Figure \ref{figydmtxyzyzview}(a) shows a
projected view of the 3-dimensional portrait
Fig.\ref{figydmtxyz}(a) to the $yw$-plane. The green,
parabola-like curves are $z$-section curves on the $y$-slow
manifold $w=q(y,z)$. The $z$-section $y$-nullcline for the
$yw$-system. The size of the $z$-section $y$-nullcline decreases
with increasing $z$. The $w$-nullcline is always a line $y=y_{\rm
wnl}$ independent of the competition intensity from $z$. If it
cuts the $z$-section $y$-nullcline left of the $y$-fold, then $w$
is $xyw$-weak for the corresponding fixed competition intensity
from $z$. Otherwise, $w$ is efficient. Note that the
$yw$-equilibrium points decreases steadily toward extinction
($w=0$) with increasing competition intensity from $z$. Again this
simplistic scenario does not extend to dynamically evolving $z$ as
we will show later.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{figydmtxyzyzview}
\centerline{\includegraphics[width=5in,height=1.8in]{figydmtxyzyzview.eps}}
\centerline{(a)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(b)} \caption{(a) $yw$-phase portraits with fixed $z$ when
$\epsilon_1=0$. (b) The projected view of Fig.\ref{figydmtxyz}(a)
to the $zw$-plane.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsection{Dynamics of Fast $Y$, Slow $W$, and Slower $Z$}
We first consider the system for which the maximum reproductive
rates ranges fast to slow in the order of $x, y, w, z$:
$0<\zeta\l1 1, 0<\epsilon_1\ll \epsilon_2\ll 1$. Biologically, the
new comer competitor $z$ is not only not $xyz$-competitive but
also not productive relative to the $xyz$-chain. We will
demonstrate below that because of $w$'s predation on $y$, $z$ is
able to survive in the form of equilibrium co-existence and
another likely form of $y$-fast cycle co-existence. A cycle is
$y$-fast if it contains $y$-fast jumps between the two stable
$y$-slow manifold.

\subsubsection{Equilibrium Coexistence} Figure \ref{figydmtxyz}
illustrates the case that $y_{\rm wnl}<y_{\rm yfd}|_{\{z=0\}}$.
The top-left figure gives a 3-dimensional illustration for the the
nullclines of the $yzw$-system in the solid ${\mathcal S}$. The
remaining ones illustrate the co-ordinate sections of the 3-d
view. If colors are visible, green is for $y$, red is for $z$,
deep blue is for $w$. For example, green curves outline the
$y$-nullcline, with in particular solid curves outline its stable
branches and dotted ones its unstable part. Since $y$ is the
$yzw$-fast variable, all singular orbits must develop in its
stable slow-manifolds with $yzw$-fast jumps in between, we
distinguish further by painting its trivial ($y=0$) stable part
light green, its unstable part ($\dot y>0$ on $w=q(y,z)$) lighter
green, and its nontrivial stable part ($\dot y<0$ on $w=q(y,z)$)
blank. Apart from a good 3-dimensional view, the most useful view
is the project $zw$-view as shown in the bottom right figure of
Fig.{\ref{figydmtxyz}}. In it the nontrivial stable $y$-slow
manifold is outlined by the blank area and the lighter green area.
The former is in fact the overlapping projections of two $y$-slow
manifold: the nontrivial stable and the trivial unstable ones. The
latter is on the other hand the overlapping projections of three
$y$-slow manifold: the nontrivial stable, the nontrivial unstable,
and the trivial stable ones. The green area is part of the trivial
stable branch only. The top boundary curve of the lighter green
area is the $y$-fold curve $(y,z,w)_{\rm yfd}$ of $w=q(y,z)$. The
dashed lower boundary plus the remaining dash green consist of the
$y$-transcritical curve $(y,z,0)_{\rm ytr}$. Half-filled circles
denote transcritical points. Squares denote some corresponding
points of Pontryagin's delay of loss of stability. Both are color
coordinated: same color squares always pair with same color
half-filled circles. Open circles are unstable equilibrium points
and filled circles are stable equilibrium points. The green curves
in the bottom left figure of Fig.\ref{figydmtxyz} are some
$z$-section curves ($z=\text{constant}$) on the $y$-nullcline
$w=q(x,y)$. It is important to note that as we have proved earlier
that both nontrivial $z$-nullcline and $w$-nullcline traces on the
stable $y$-slow manifold (both open and lighter green areas in the
$zw$-view of Fig.\ref{figydmtxyz}) are monotone decreasing in $z$.
The $z$-nullcline always traces from $z=0$ to the
$y$-transcritical curve $(y,z)_{\rm ytr}$. The $w$-nullcline curve
always has its $w$-lower end point on $w=0$. The
$xyzw$-equilibrium point ${\mathcal E}$ is the intersection of
these $z$-, $w$-nullcline curves on the $y$-slow manifold. Because
the nontrivial $z$-, $w$-nullcline surfaces are both planes
parallel to the $w$-axis, they intersect along a line parallel to
the $w$-axis, if they do. Therefore this $zw$-nullcline line
intersects the $y$-slow manifold at a unique point ${\mathcal E}$.
The condition that $z$-, $w$-nullcline surfaces intersect is
\begin{equation}\label{cndtzwintsct}
y_{\rm wnl}\le y_{\rm znl}|_{\{z=0\}},
\end{equation}
see Fig.\ref{figydmtxyz}'s $yz$-projected illustration.

Assume the condition (\ref{cndtzwintsct}) for the existence of the
$xyzw$-equilibrium point ${\mathcal E}$. For it to be stable, it
is necessary to have the intersection to occur on the stable
branch of the $y$-slow manifold. For this to happen, the
$w$-higher end point of the $w$-nullcline trace must lie above the
$z$-nullcline. Figure \ref{figcoeq} illustrates the two cases: the
$w$-high $w$-nullcline point is not the fold Fig.\ref{figcoeq}(a)
and is. Condition for the former can be given explicitly as
\begin{equation}\label{xyzweqptcnd1}
w_{\rm wnl}=q(y_{\rm wnl},0)>w_{\rm znl}|_{\{z=0\}}.
\end{equation}
Condition for the latter is given implicitly as follows since we
do not have an explicit expression of the $y$-nullcline fold curve
$(y,z,w)_{\rm yfd}$:
\begin{equation}\label{xyzweqptcnd2}
w_{\rm yfd}|_{\{y=y_{\rm wnl}\}}>w_{\rm znl}|_{\{y=y_{\rm wnl}\}}.
\end{equation}
where $w_{\rm yfd}|_{\{y=y_{\rm wnl}\}}$ is the $w$-coordinate of
the intersection of the $y$-fold and $w$-nullcline $\{y=y_{\rm
wnl}\}$, if they intersect nonempty.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{figcoeqpt}
\centerline{\includegraphics[width=5in,height=1.8in]{figcoeqpt.eps}}
\centerline{(a)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(b)}
\caption{Typical $yzw$-singular orbits for cases of
$0<\epsilon_1\ll \epsilon_2\ll 1$ with stable $xyzw$-equilibrium
point ${\mathcal E}$: (a) The $w$-nullcline does not go through
the $y$-fold}; (b) The $w$-nullcline goes through the $y$-fold.
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{pro}\label{propcoexteqpt} Under the conditions (\ref{xyzweqptcnd1}) or
(\ref{xyzweqptcnd2}) and for sufficiently small $0<\epsilon_1\ll
\epsilon_2\ll 1$, the $xyzw$-equilibrium point ${\mathcal E}$ is
locally stable. Moreover, at the singular limit
$\epsilon_1/\epsilon_2=0$ all singular orbits converge to
${\mathcal E}$ for the case of (\ref{xyzweqptcnd1}).
\end{pro}

\begin{proof} We consider first the singular case that $w$ is fast
relative to $z$. As illustrated in Fig.\ref{figcoeqpt}(a) for the
case (\ref{xyzweqptcnd1}). Initial points 1, 3 are the nontrivial
stable $y$-slow manifold. For type 1 initial points, singular
orbits jump to the nontrivial $w$-nullcline first and then
converge to ${\mathcal E}$ following the $z$-slow flow. For type 3
initial points, singular orbits jump to the trivial $w$-nullcline
$w=0$, cross the $w$-transcritical point, jump up the nontrivial
$w$-nullcline at some $w$-PDLS(Pontryagin's delay of loss of
stability) points, and then converge to ${\mathcal E}$ like type 1
points do. Type 2 points start on the trivial $y$-slow manifold.
They are represented by dash arrows. They will either jump to the
nontrivial stable $y$-nullcline to become type 1 or 3 points, or
they will first cross the $y$-transcritical points, jump to the
nontrivial stable $y$-slow manifold at some $y$-PDLS points to
become type 1 or 3 points. Type 4 points have an additional
feature. They jump to the trivial $w$-nullcline first, then
following the $z$-slow flow to cross the $y$-transcritical point
and jump to the nontrivial stable $y$-nullcline at some $y$-PDLS
points before converging to ${\mathcal E}$. Hence, all singular
orbits converge to ${\mathcal E}$. The local stability of
${\mathcal E}$ at the singular limit $\epsilon_1/\epsilon_2=0$ is
similarly argued for both cases concerning only initial points of
type 1 as illustrated in Fig.\ref{figcoeqpt}(b). By Fenichel's
geometric theory of singular perturbation (\cite{Fen79}) the local
stability of ${\mathcal E}$ persists for small
$0<\epsilon_1/\epsilon_2\ll 1$. This completes the proof.
\end{proof}

For non-coexistence in equilibrium points, we are left with two
cases. Either there is no $xyzw$-equilibrium point if condition
(\ref{cndtzwintsct}) is violated: $y_{\rm wnl}\le y_{\rm
znl}|_{\{z=0\}}$, or ${\mathcal E}$ is unstable because
(\ref{xyzweqptcnd2}) is violated, in which case ${\mathcal E}$
lies on the unstable part of the nontrivial $y$-nullcline, see
Fig.\ref{figyfstcocycl}(a). These cases are under the current rate
assumption that $0<\epsilon_1\ll \epsilon_2\ll 1$. As we will see
below, ${\mathcal E}$ will become unstable if the rate condition
is reversed to $0<\epsilon_2\ll \epsilon_2\ll 1$ even if condition
(\ref{xyzweqptcnd1}) holds as in Proposition\ref{propcoexteqpt}.
For the former case, it is straightforward to check the following
statement.

\begin{pro}\label{propzdiesout} If condition (\ref{cndtzwintsct})
does not hold, i.e., $y_{\rm wnl}> y_{\rm znl}|_{\{z=0\}}$, then
all singular orbits converge to the $xyw$ equilibrium point with
$z$ dies out.
\end{pro}

\subsubsection{Coexistence in $Y$-fast Cycle} We can conclude from
Propositions \ref{propcoexteqpt} and \ref{propzdiesout} that one
necessary condition for possible existence of a periodic
$xyzw$-cycle is when condition (\ref{xyzweqptcnd2}) holds because
under condition (\ref{xyzweqptcnd1}) the equilibrium point
${\mathcal E}$ globally attracts all singular orbits with
non-vanishing initials. Under condition (\ref{xyzweqptcnd2}), the
long term $xyw$-dynamics is the $x$-slow, $y$-fast cycle discussed
in Sec.???. Although $z$ is not $xyz$-competitive by our primary
assumption that $E_{xy}<0$, it can be $xyzw$-competitive if the
$xyw$-cycle is unstable. For example, if the $xyw$-cycle lies
behind the $z$-nullcline, i.e., $w_{\rm znl}|_{\{z=0\}}<w_{\rm
ypd}$, with $w_{\rm ypd}$ denoting the $y$-PDLS point of the
$xyw$-cycle, then it is unstable because $z$ increases near it. As
we concluded from Sec.??? that the analytical condition for the
$xyzw$-stability of the $xyw$-cycle is given by the stability of
$z=0$ for the variational equation $\dot z=zh(x_c(t))$ with
$x_c(t)$ denoting the $x$-coordinate of the $xyw$-cycle. In
particular the $xyw$-cycle is stable if $C_{xyw}:=\int_0^T
h(x_c(t)dt<0$ with $T$ being the cycle's period. Since we have not
obtained an explicit expression for the $y$-fold, an explicit
expression for the stability criterion is not possible at this
time. Nonetheless, we can conclude the following statement
(already under the condition that both $y, z$ are weak, $y$ is
$xyz$-dominant).

\begin{pro}\label{proyfstcycl} If $z$ is
$xyzw$-competitive, i.e.,$C_{xyw}>0$, and ${\mathcal E}$ is
unstable because condition ({xyzweqptcnd2}) is not reversed: $
w_{\rm yfd}|_{\{y=y_{\rm wnl}\}}<w_{\rm znl}|_{\{y=y_{\rm
wnl}\}}$, then there exists a non-equilibrium co-existing
$xyzw$-attractor.
\end{pro}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{figyfstcocycl}
\centerline{\includegraphics[width=5in,height=1.8in]{figyfstcocycl.eps}}
\ \centerline{(a)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(b)} \caption{(a) Existence of a non-equilibrium point coexisting
$xyzw$-attractor. (b) Existence of two co-existing
$xyzw$-attractor.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In fact without loss of generality, assume all singular orbits
start out on the nontrivial stable $y$-slow manifold. For all
initials attracted to the nontrivial stable $w$-nullcline, they
move in the decreasing direction of $z$ when they hit
$w$-nullcline. For initials near the $xyw$-cycle on $z=0$,
perturbed singular orbits drift in the increasing direction of $z$
because the $xyw$-cycle is unstable. Since ${\mathcal E}$ is
unstable, all singular orbits must converge to a nonequilibrium
co-existing $xyzw$-attractor. See Fig.\ref{figyfstcocycl}(a). We
conjecture that this attractor is a $y$-fast cycle.

\subsubsection{Multiple Coexisting Attractors} Any one of the
conditions of Proposition \ref{proyfstcycl} can go through a
bifurcation. Assuming at the bifurcation point of condition
(\ref{xyzweqptcnd2}) $w_{\rm yfd}|_{\{y=y_{\rm wnl}\}}=w_{\rm
znl}|_{\{y=y_{\rm wnl}\}}$, i.e., the equilibrium point ${\mathcal
E}$ lies on the $y$-fold, the $xyw$-cycle is still unstable and
there is a non-equilibrium co-existing $xyzw$-attractor. If this
$xyzw$-attractor is an asymptotically stable limit cycle as
conjectured above, then under a small perturbation to condition
(\ref{xyzweqptcnd2}): $w_{\rm yfd}|_{\{y=y_{\rm wnl}\}}>w_{\rm
znl}|_{\{y=y_{\rm wnl}\}}$, ${\mathcal E}$ becomes stable. At the
same slightly perturbed parameter value, the co-existing
$xyzw$-cycle also persists. Thus it results in a dynamics for
which there are at least 2 $xyzw$-attractors: the stable
equilibrium ${\mathcal E}$ and the $xyzw$-cycle. Figure
\ref{figyfstcocycl}(b) illustrates such a situation, in which the
point $p$ at the $y$-fold and $w$-nullcline drifts leftwards when
the singular perturbation $0<\epsilon_1\ll \epsilon_2$ is on as
shown.

\subsection{Dynamics of Fast $Y$, Slow $Z$, and Slower $W$}
We now consider the system for which the maximum reproductive
rates ranges fast to slow in the order of $x, y, z, w$:
$0<\zeta\l1 1, 0<\epsilon_2\ll \epsilon_1\ll 1$. Biologically, the
new comer competitor $z$ is not $xyz$-competitive but more
productive than the top predator $w$. We will demonstrate below
that because of the rate reversal between $z$ and $w$,
co-existence in the form of equilibrium point is not possible.
Instead, for moderate rate comparison
$\epsilon_1/\epsilon_2=O(1)$, the stable equilibrium point
${\mathcal E}$ undergoes a Hopf bifurcation to give rise to a
$y$-slow cycle. It grows with the time rate more divergent
$0<\epsilon_2\ll \epsilon_1\ll 1$. It can either grow into a large
singular $y$-slow cycle at $\epsilon_2/\epsilon_1=0$ or chaotic
attractor.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{fighopf}
\centerline{\includegraphics[width=5in,height=1.8in]{fighopf.eps}}
\ \centerline{(a)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(b)} \caption{(a) A $y$-slow cycle born through Hopf bifurcation
at ${\mathcal E}$. (b) It grows into a large singular cycle
showing the case of $w$-nullcline going through the $y$-fold. In
the case of (a) and as $\epsilon_2/\epsilon_1\to 0$, the cycle
grows but will never surpass the level $w=w_{\rm wnl}$, at which
the $xyw$-equilibrium point shown is not stable.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsubsection{Coexisting $Y$-Slow Cycle via Hopf Bifurcation}
Assume the same condition (\ref{cndtzwintsct}) that guarantees the
existence of the nontrivial equilibrium point ${\mathcal E}$. By
Proposition \ref{propcoexteqpt} we know that for small
$0<\epsilon_1\ll \epsilon_2\ll 1$, ${\mathcal E}$ is locally
stable under the condition (\ref{xyzweqptcnd1}) or
(\ref{xyzweqptcnd2}). Then the following result holds.

\begin{pro}\label{prophopf} Under the conditions (\ref{xyzweqptcnd1}) or
(\ref{xyzweqptcnd2}), then for each $0<\epsilon_2\ll 1$, there is
a value $\epsilon_1=\theta(\epsilon_2)$ such that a pair of
eigenvalues of the linearization at ${\mathcal E}$ crosses the
imaginary axis.
\end{pro}

\begin{proof} By Fenichel's geometric theory of singular
perturbation, the invariant $y$-slow manifold in a compact
neighborhood near ${\mathcal E}$ persists for small
$0<\epsilon_2\ll 1$ treating $0\epsilon_1/\epsilon_2\ll 1$ as a
time scale independent of $\epsilon_2$. On this 2-dimensional
manifold the $yzw$-system can be projected onto the $zw$-plane. In
a neighborhood of the ${\mathcal E}$, the $zw$-phase plane is
qualitatively like Fig.\ref{fighopf}(a). In particular, when
linearized, the equation can be written as
\begin{equation*}
\begin{split}
\dot u&=\epsilon_1c_1(au+v)\\
\dot v&=-\epsilon_2c_2(bu+v)
\end{split}
\end{equation*}
with $(u,v)=(0,0)$ corresponding to the equilibrium point
${\mathcal E}$ in $(z,w)$. Here, $a,b,c_1,c_2$ are constants
satisfying the following conditions. Because $u$-equation's right
side is the linearization of the vector field $\epsilon_1zh$
restricted on the $y$-slow manifold, the $u$-nullcline $au+v=0$ is
precisely the tangent line to the $z$-nullcline at ${\mathcal E}$.
Therefore, it has negative slope, forcing $a>0$. Exactly same
argument results in $b>0$. In addition, the relative position of
these $u$-, $v$-nullclines preserves that of its $z$-,
$w$-nullclines, we must have $b>a>0$ because $w$ decreases faster
on the $w$-nullcline than on the $z$-nullcline. Since $\dot z>0$
for points above the $z$-nullcline, which include points of large
$w$ as biologically the presence of $w$ enhances the
competitiveness of $z$. Since the linearization in terms of $u$
mirrors the same qualitative property, we must have the constant
$c_1>0$ to be positive. Exactly the same argument leads to
$c_2>0$. Since $\epsilon_i$ and $\epsilon_ic_i$ have the same
order of magnitude for each $i=1,2$, we will let
$\epsilon_i=\epsilon_ic_i$ for simpler notations in the following
calculation. Now the eigenvalues $\lambda_j$ of the linear
$uv$-system are
\[
\lambda_j=\frac{1}{2}[\epsilon_1a-\epsilon_2\pm\sqrt{(\epsilon_1a-\epsilon_2)^2
-4\epsilon_1\epsilon_2(b-a)}].
\]
Then we see easily that the eigenvalues are pure imaginaries at
$\epsilon_1a-\epsilon_2=0$, and their real parts change from
negative to positive as $\epsilon_1$ increases above
$\epsilon_2/a$. This completes the proof.
\end{proof}

This result implies that as the equilibrium point ${\mathcal E}$
changes its stability from stable to unstable as $\epsilon_1$
crosses the Hopf point, there must be a stable periodic orbit
bifurcating from the equilibrium point ${\mathcal E}$ generically,
see Fig.\ref{fighopf}(a) for an illustration. In fact in the limit
$\epsilon_2/\epsilon_1\to 0$ with $z$ as the singular variable
faster than $w$, we can demonstrate below the existence of a large
singular cycle under certain conditions.

\subsubsection{Large $YZ$-Fast Singular Cycle} First we describe
the key condition that will allow us to trap such an orbit. Start
at the $y$-fold point $(y_{\rm yfd},0,w_{\rm yfd})$ on $z=0$, c.f.
Fig.\ref{fighopf}(b). Following the $y$-fast flow, the singular
orbit lands on the trivial $y$-nullcline surface $y=0$ with the
other coordinates unchanged. Then the $z$-fast dynamics takes
over, sending the orbit horizontally to the attracting
$z$-nullcline $y=0,z=z_{\rm znl}$, landing on it with the same
$w$-coordinate $w_{\rm yfd}$. See both
Figs.\ref{fighopf}(b),\ref{figchaos}(a). Once on the line
$y=0,z=z_{\rm znl}$, the $w$-dynamics takes over, and the slow
$w$-orbit develops downwards. It first crosses the
$y$-transcritical point $(y,z,w)=(0,z_{\rm znl},w_{\rm ytr})$, and
must then jump to the nontrivial $y$-slow manifold $w=q(y,z)$ at a
$y$-PDLS point $(0,z_{\rm znl},w_{\rm ypd})$ due to Pontryagin's
delay of loss of stability with respect to the fast $y$-flow. The
value of $w_0:=w_{\rm ypd}|_{\{z=z_{\rm znl}\}}$ is determined as
follows.

First we need to be more specific about the time scales of the
reduced $yzw$-system in the solid ${\mathcal S}$. Rescale the time
$t\to \epsilon_2 t$ to obtain the reduced system as
\[
\epsilon_2 y'=yg(p(y,z),y,w),\ \epsilon_0 z'=zh(p(y,z)),\ w'=wk(y)
\]
with $\epsilon_0=\epsilon_2/\epsilon_1\ll 1$ as an independent
singular parameter from $\epsilon_2$. Thus, under the condition
that $0<\epsilon_2\ll \epsilon_1\ll 1$, the $y$ is the fastest of
the three at a rate of order $O(1/\epsilon_2)$, the $z$ is the
second fastest at a rate of order $O(
1/\epsilon_0)=O(\epsilon_1/\epsilon_2)\ll O(1/\epsilon_2)$.

Take any plane $y=y_0>0$ sufficient near $y=0$. Then a perturbed
orbit starting from $(y_0,z_0,w_{\rm yfd})$ for $z_0>0$ is first
attracted to the trivial $y$-nullcline $y=0$, only to emerge
somewhere below the $y$-slow manifold $w=q(y,z)$ and move toward
the manifold. In doing so it will hit the plane $y=y_0$ again in
direction opposite to when it starts. Hence following the orbit,
the integral identity below holds
\[
0=\int_{y_0}{y_0}\frac{1}{\epsilon_2y}dy.
\]
Since we can choose $0<y_0<y_{\rm wnl}$, the $w$-component of this
portion of the orbit always decreasing in $w$. Thus, we can make
the following substitution of variable from $y$ to $w$
\[
0=\int_{y_0}{y_0}\frac{1}{\epsilon_2y}dy=\int_{w_{\rm
yfd}}^{w_{\rm ypd}}\frac{1}{\epsilon_2y}\frac{dy}{dw}dw.
\]
Simplify the integral to obtain
\[
0=\int_{w_{\rm yfd}}^{w_{\rm ypd}}\frac{g(p(y,z),y,w)}{wk(y)}dw.
\]
Taking the limit $\epsilon_2\to 0$ first, we have $y=0$ on the
equation above. Taking the limit $\epsilon_0\to 0$ next for the
reduced $zw$-system, we have
$z=z_{\rm znl}$ for $z$ of the integrant. %By the definition of function $p$, we have
%$p(0,z_{\rm znl})=x_{\rm znl}$.
Therefore at the limit $\lim_{\epsilon_0\to 0}\lim_{\epsilon_2\to
0}$, the integral equation above becomes
\[
0=\int_{w_{\rm yfd}}^{w_{\rm ypd}}\frac{g(p(0,z_{\rm
znl}),0,w)}{wk(0)}dw,
\]
defining the $w_{\rm ypd}$ value on $z=z_{\rm znl}$ as shown in
Fig.\ref{fighopf}(b) as well as \ref{figchaos}(a). Because
$g(p(0,z_{\rm znl}),0,w)<0$ for $w>w_{\rm ytr}$ and $g(p(0,z_{\rm
znl}),0,w)>0$ for $0<w<w_{\rm ytr}$, and the fact that the
indefinite integral
\[
\int_{w_{\rm ytr}}^{u}\frac{g(p(0,z_{\rm znl}),0,w)}{wk(0)}dw
\]
is $0$ at $u=w_{\rm ytr}$ and $+\infty$ at $u=0$, the integral
equation above for the $y$-PDLS point $w_{\rm ypd}$ has a unique
solution between $0$ and $w_{\rm ytr}$.


Once it is on the $y$-slow manifold $w=q(y,z)$, the $z$-dynamics
takes over, sending the singular orbit to the trivial
$z$-nullcline on $w=q(y,0)$ with $w=w_0=w_{\rm ypd}|_{\{z=z_{\rm
znl}\}}$. The slow $w$-dynamics takes over subsequently, moving
the orbit upward. Notice that since the non-trivial $z$-nullcline
is between the $y$-fold curve and the boundary $0=q(y,z)$, the
non-trivial $z$-nullcline strikes through the $y$-slow manifold
$w=q(y,z)$ from $y=0$ to $z=0$. As a result it intersects on
$w=q(y,0)$ at a $z$-PDLS point between $w=0$ and $w=w_{\rm yfd}$
as shown in Fig.\ref{fighopf}(b) labelled with $w_{\rm znl}$. Also
between $w=0$ and the $y$-fold point there is no
$xyzw$-equilibrium point because of the weak predator assumption
on $w$.  We can conclude now that after its passing the
$z$-transcritical point $w_{\rm znl}$ on $z=0,w=q(y,0)$, the
$z$-dynamics becomes unstable. Pontryagin's delay of loss of
stability takes place for the singular orbit so that soon or later
it will develop in the increasing direction of $z$.

To be more precise on how the singular orbit turns around in
variable $z$, we again use the PDLS argument. Let $w_0=w_{\rm
ypd}|_{\{y=0,z=z_{\rm znl}\}}$ as above. Let $z_0>0$ be
sufficiently near $z=0$. A perturbed orbit starting from the plane
$z=z_0$ at point $(y_0,z_0,w_0)$ with $y_0$ defined by
$w_0=q(y_0,z_0)$ will eventually emerge away from $z=0$ and
strikes the plane $z=z_0$ in opposition direction. Let $w_{\rm
zpd}$ denote the $w$-coordinate of the point of return. Then, the
identity $0=\int_{z_0}^{z_0}\epsilon_0 dz/z$ holds. Suppose the
orbit does not pass the $w$-nullcline $y=y_{\rm wnl}$, then one
substitution from variable $z$ to $w$ is sufficient to change the
integral equation to
\begin{equation}\label{eqzpdptofyfd}
0=\int_{w_0}^{w_{\rm zpd}}\frac{\epsilon_0}{z}\frac{dz}{dw}dw=
\int_{w_0}^{w_{\rm zpd}}\frac{h(p(y,z))}{wk(y)}dw.
\end{equation}
Taking the double limit $\lim_{\epsilon_0\to 0}\lim_{\epsilon_2\to
0}$, the equation above becomes
\[
0=\int_{w_0}^{w_{\rm
zpd}}\frac{h(p(y,0))}{wk(y)}dw=\int_{w_0}^{w_{\rm
zpd}}\frac{h(p(q^{-1}(w),0))}{wk(q^{-1}(w))}dw,
\]
with $y=q^{-1}(w)$ uniquely defined by the invertible function
$w=q(y,0)$.

\begin{pro}\label{proplargehopf} If the $z$-PDLS point $w_{\rm
zpd}$ defined by the equation above lies below the $y$-fold point
$w_{\rm yfd}$, then there exists a singular $yz$-fast cycle.
\end{pro}

\begin{proof} Take the interval $I=[w_{\rm ypd},w_{\rm ytr}]$ on
$z=0,w=q(y,0)$ as shown in Fig.\ref{fighopf}(b). A return map
following the singular orbits of the interval can be defined.
Under the assumption that $w_{\rm zpd}$ never surpasses the
$y$-fold point $w_{\rm yfd}$, the lower end point $w_0=w_{\rm
ypd}$ returns to the interval $I$. Because the $z$-nullcline
decreases in $w$ with increase in $z$, we see clearly from
Fig.\ref{fighopf}(b) that the upper end point $w_{\rm ytr}$ must
also returns to $I$ as well. The continuity of the this return map
is obvious. Also it is strictly increasing because of the way all
PDLS points are defined their respect integral equations whose
solutions are unique. Hence there must be a unique fixed point of
the return map, and the fixed point indeed corresponds to a
$yz$-fast cycle. This completes the proof.
\end{proof}

The condition of this result must holds if the $w$-nullcline does
not go through the $y$-fold as with the case of
Fig.\ref{fighopf}(a). Because in such a case, the $w$-nullcline
intersects $w=q(y,z),z=0$ at an unstable equilibrium point. As a
result the $z$-PDLS point defined in the proposition will never
pass the unstable equilibrium point. In fact, the integral
equation defining $w_{\rm zpd}$ above can be separated into two
parts as
\[
0=\int_{w_0}^{w_{\rm zpd}}\frac{h(p(y,0))}{wk(y)}dw=
\int_{w_0}^{w_{\rm znl}}\frac{h(p(y,0))}{wk(y)}dw + \int_{w_{\rm
znl}}^{w_{\rm zpd}}\frac{h(p(y,0))}{wk(y)}dw,
\]
with $q^{-1}(w)$. The first integral is negative and finite. The
second integral is positive. The denominator of the integrant of
the latter
$\frac{h(p(y,0))}{wk(y)}\frac{h(p(q^{-1}(w),0))}{wk(q^{-1}(w))}$
diverges as $w$ increases to $w_{\rm wnl}$ at which
$k(q^{-1}(w))=0$. Hence the integral equation above must have a
solution for $w_{\rm zpd}<w_{\rm yfd}$ as required by Proposition
\ref{proplargehopf}'s assumption.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{figchaos}
\centerline{\includegraphics[width=5in,height=4in]{figchaos.eps}}
\caption{(a) The $z$-PDLS point for $w_0$ from $\Sigma$ is on
$y=z=0$ at $w=w_{\rm zpd}$. (b) Due to the $y$-fold point $w_{\rm
yfd}$ the return map $\pi$ develops a local minimum at the
pre-PDLS point $w_1$ of $w_{\rm zpd}$. (c) The bifurcation of the
local minimum through $w_0=w_1$. (d) Possible multiple extrema for
small $w_0$. (e) Combination PDLS jumps on singular orbits with
$\epsilon_2=0$.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsubsection{Co-existence in Chaos}
We will demonstrate in this subsection that when the condition of
Proposition \ref{proplargehopf} is violated the system can
bifurcate into chaotic dynamics. To do this we start from the
integral equation $0=\int_{z_0}^{z_0}\epsilon_0 dz/z$ that defines
the $z$-PDLS point $w_{\rm zpd}$ starting at $w_0$. Recall that
the assumption that the perturbed orbit from $w_0$ on $w=q(y,z)$
never crosses the $w$-nulllcline plane $y=y_{\rm wnl}$ gives rise
to the integral equation (\ref{eqzpdptofyfd}). Now suppose instead
that the orbit crosses the $w$-plane only once with $w$-coordinate
at the crossing as $w^\ast$. Then the orbit is increasing in $w$
from $w_0$ to $w^\ast$ and decreasing from $w^\ast$ to $w_{\rm
zpd}$. Then two separate substitutions from $z$ to $w$ are needed
to bring the equation $0=\int_{z_0}^{z_0}\epsilon_0 dz/z$ into
\begin{equation*}
\begin{split}
0&=\int_{w_0}^{w^\ast}\frac{\epsilon_0}{z}\frac{dz}{dw}dw
+\int_{w^\ast}^{w_{\rm zpd}}\frac{\epsilon_0}{z}\frac{dz}{dw}dw\\
&=\int_{w_0}^{w^\ast}\frac{h(p(y,z))}{wk(y)}dw+
\int_{w^\ast}^{w_{\rm zpd}}\frac{h(p(y,z))}{wk(y)}dw.
\end{split}
\end{equation*}
Although the two integrals look the same in form, the limiting
behaviors with $\epsilon_2\to 0$ are different. In the first
integral all points lie on the non-trivial $y$-slow manifold
$w=q(y,z)$ as $\epsilon_2\to 0$. In the second all points lie on
the trivial $y$-slow manifold $y=0$ instead. Moreover, the
$y$-fast dynamics forces $w^\ast\to w_{\rm yfd}$. Follow up by
taking $\epsilon_0\to 0$, then $z=0$ in both integrals. Therefore,
with the double limit $\lim_{\epsilon_0\to 0}\lim_{\epsilon_2\to
0}$, the integral equation becomes
\begin{equation}\label{eqzpdptofyfdoveryfd}
\begin{split}
0&=\int_{w_0}^{w_{\rm yfd}}\frac{h(p(y,0))}{wk(y)}dw+ \int_{w_{\rm
yfd}}^{w_{\rm zpd}}\frac{h(p(0,0))}{wk(0)}dw\\
&=\int_{w_0}^{w_{\rm znl}}\frac{h(p(y,0))}{wk(y)}dw\\
&\quad +\left[\int_{w_{\rm znl}}^{w_{\rm
yfd}}\frac{h(p(y,0))}{wk(y)}dw+\int_{w_{\rm yfd}}^{w_{\rm
zpd}}\frac{h(p(0,0))}{wk(0)}dw\right]\\
&:=I_1+I_2,
\end{split}
\end{equation}
with $y=q^{-1}(w)$ defined from $w=q(y,0)$ as before, where
$I_1=\int_{w_0}^{w_{\rm znl}}\frac{h(p(y,0))}{wk(y)}dw$ and $I_2$
denote the two integrals in the bracket. See
Fig.\ref{figchaos}(a). Notice that $I_1$ is negative which can
diverge to $-\infty$ if $w_0\to 0$. Whereas $I_2$ is positive,
bounded above. This analysis implies that the condition of
Proposition \ref{proplargehopf} must bifurcates into such a
situation prescribed by equation (\ref{eqzpdptofyfdoveryfd}) as
$w_0=w_{\rm ypd}|_{\{y=0,z=z_{\rm znl}\}}$ drops lower.

If we let $\Sigma$ be an $w$-interval on $w=q(y,z)$ as shown in
Fig.\ref{figchaos}(a) whose $w$-end points are the same as the
interval $I=[w_{\rm ypd},w_{\rm ytr}]$ from the proof of
Proposition \ref{proplargehopf}. For simpler notation let
$w_2=w_{\rm ytr}|_{\{y=0,z=z_{\rm znl}\}}$ and let $w_1$ be the
pre-PDLS point of the $y$-fold $w_{\rm yfd}$, i.e.
$\int_{w_1}^{w_{\rm yfd}}\frac{h(p(y,0))}{wk(y)}dw=0$. Then the
same singular orbit induced return map used in the proof
Proposition \ref{proplargehopf} is no longer monotone increasing.
Let $\pi$ denote the return map. Then $w_1$ is a critical point in
$\Sigma=[w_0,w_2]$. It increases in $[w_1,w_2]$ and decreases in
$[w_0,w_1]$ as depicted in Fig.\ref{figchaos}(b).
Fig.\ref{figchaos}(c) shows the case at bifurcation when
$w_0=w_1$.

If $w_0$ is just slightly perturbed below the bifurcation point
$w_0=w_1$, then the $z$-PDLS point $w_{\rm zpd}$ is perturbed
slightly below $w_{\rm yfd}$, and the image $\pi(w_0)$ of the
return map lies a little bit above $w_0$. Further decreasing $w_0$
will result in significant increase in $\pi(w_0)$, which may give
rise to chaotic dynamics of the return map $\pi$. Decreasing $w_0$
further can make the return map more complex by creating more
critical points as shown in Fig.\ref{figchaos}(d). In fact,
decreasing $w_0$ makes the negative integral $I_1$ of
(\ref{eqzpdptofyfdoveryfd}) more negative. To counter balance
$I_1$, $w_{\rm zpd}$ must further move downward on the $w$-axis
and come to meet the $y$-PDLS point $w_{\rm ypd}$. At this point
the $y$-fast dynamics must take over to switch the singular orbit
to the stable $y$-slow manifold $w=q(y,z)$. Further decreasing
$w_0$ moves the $w_{\rm zpd}$ upwards on the nontrivial stable
$y$-nullcline $w=q(y,0)$, splitting the integral $I_2$ into three
parts
\[
I_2=\int_{w_{\rm znl}}^{w_{\rm
yfd}}\frac{h(p(y,0))}{wk(y)}dw+\int_{w_{\rm yfd}}^{w_{\rm
ypd}}\frac{h(p(0,0))}{wk(0)}dw+\int_{w_{\rm ypd}}^{w_{\rm
zpd}}\frac{h(p(y,0))}{wk(y)}dw.
\]
Under this situation, another critical point emerges between $w_0$
and $w_1$ to become a local maximum for the return map $\pi$. Now
it is clear to see that further decreasing $w_0$ will make the
$z$-PDLS point $w_{\rm zpd}$ to move up on $w=q(y,0)$ and down on
$y=z=0$, creating more local extrema. Fig.\ref{figchaos}(d) shows
such a possible situation.

One sufficient condition for chaotic dynamics is to have three
adjacent critical points that straddle the diagonal $\pi=w$ in the
sense that the value of the middle point under the return map is
on one side of the diagonal $\pi=w$ and those of the other two are
on the opposite side of the diagonal so that these critical points
also lie between the minimum and maximum values of the points.
More precisely we have the following partial result.

\begin{pro}\label{propchaos} If there exist three adjacent
critical points of $\pi$ that straddle the diagonal, then there is
a subset of $\Sigma$ on which $\pi$ is equivalent to the shift
dynamics on two symbols.
\end{pro}

\begin{proof} A proof for such a continuous map is given in
\cite{DH022}.
\end{proof}

Numerical simulations indeed suggest such a scenario of
Proposition \ref{propchaos} as well as chaos due to the existence
of a single critical point for a quadratic-like map $pi$.

Finally we conclude this subsection by explaining a special case
of Pontryagin's delay of loss of stability when the slow flow is
2-dimensional. The $\epsilon_2$-fast dynamics is that of $y$, and
the $\epsilon_2$-slow dynamics is that of $z,w$. Start a point on
the $y$-fold as shown in Fig.\ref{figchaos}(e). The singular
$\epsilon_2$-fast dynamics sends it down on $y=0$. The
$\epsilon_2$-slow $zw$-dynamics takes over. It moves down in $w$,
crosses the $y$-transcritical curve, and jumps back to $w=q(y,z)$
at some $y$-PDLS point $w_{\rm ypd}$ as shown, and zigzags its way
toward $y=0,z=z_{\rm znl}$ afterward. Since we always pick
non-zero singular values $0<\epsilon_2\ll \epsilon_1\ll 1$ for
numerical simulations, such zigzag orbits are bound to occur as
shown in Fig.\ref{?????}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{figfstyzcycle}
\centerline{\includegraphics[width=3in,height=3.3in]{figfstyzcycle.eps}}
\caption{Singular orbits consist of fast jumps between the
$yz$-equilibrium branches $y=0,z=z_{\rm znl}$ and $z=0,w=q(y,0),
y\ge y_{\rm yfd}$ on which slow motion proceeds down and up
respectively.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsection{Dynamics of Fast $Y$ and $Z$ but Slow $W$} In terms of the singular
parameters, we assume in this case that $0<\epsilon_2\ll 1,
\epsilon_1=O(1)$ for the trophic time scales. With this
assumption, the reduced $yzw$-system has only 2 time scales in the
$x$-slow manifold solid ${\mathcal S}$: the $\epsilon_2$-fast
$yz$-dynamics and the $\epsilon_2$-slow $w$-dynamics. Recall that
for each frozen value $w=w_0$ with $\epsilon_2=0$ in
Eq.(\ref{ScaledSystem}), there are three typical cases: (1) For
$0<w_0<w_{\rm znl}|_{\{y=0\}}$, the $xy$-equilibrium point is
globally attracting in the solid ${\mathcal S}$ as in
Fig.\ref{figydmtxyz}(b). (2) For $w_{\rm
znl}|_{\{y=0\}}<w_0<w_{\rm znl}|_{\{z=0\}} $, both
$xy$-equilibrium point and $xz$-equilibrium point are locally
stable with their basins of attraction separated by a separatrix
through the $xyz$ equilibrium point which is unstable, as in
Fig.\ref{figydmtxyz}(c). (3) For $w_{\rm znl}|_{\{z=0\}}<w_0$, $z$
dominates, all nontrivial solutions converge to the
$xz$-equilibrium point as in Fig.\ref{figydmtxyz}(d). The
remaining cases are not typical, they are bifurcations from one of
the three cases (1,2,3) to another. When the $yzw$-system is
considered, the asymptotic singular orbit is rather simple: a
cycle exactly alike the large $yz$-cycle discussed above. Actually
we can trap the limit cycle in a similar way as illustrated in
Fig.\ref{figfstyzcycle}. For example, let $I=[w_{\rm znl},w_{\rm
yfd}]$ be the interval on $y=0,z=z_{\rm znl}$ as shown. A singular
orbit starting at the top end point $w_{\rm yfd}$ moves down
driven by the $\epsilon_2$-slow $w$-system. It passes the $y$-PDLS
point $w_{\rm znl}$, entering the $w$-region (1). The
$xz$-equilibrium point of the $yz$-system loses its stability with
respect to the $y$-direction. The $w$-slow orbit switches to the
$yz$-fast orbit only at the $y$-PDLS level $w=w_{\rm ypd}$.
Following the $yz$-flow, it flies to the dominant $xy$-equilibrium
point on $z=0$. Once it is on the $yz$-slow equilibrium manifold
$w=q(y,0),z=0$, the $w$-slow flow takes over again, developing
upward because of the partition by the $w$-nullcline $y=y_{\rm
wnl}$. Upon passing the $z$-transcritical point $w_{\rm ztr}$, the
$xy$-equilibrium point loses its stability in its $z$-direction.
It must then switch to the $yz$-fast dynamics at a level not
higher than the $y$-fold point $w_{\rm yfd}$, in the $w$-region
(3). In fact, the $y$-fold point is a $yz$-fold point. The
$yz$-fast system loses its stability in both directions. Because
$y,z$ have equal magnitude in time scale, any $w$-slow orbit must
immediately switch to the $yz$-fast system if it reaches the fold
point at all. Therefore, all singular orbits from the interval $I$
will not surpass the fold point and must return to $I$. Notice
that the lower end point $w_{\rm znl}$ obviously returns to $I$ as
shown since the nontrivial $z$-nullcline on $w=q(y,z)$ decreases
in $w$ with increasing $z$. In other words, the $w$-region (2) is
inaccessible for points of $I$. Also, by the way the PDLS points
are calculated from their respective integral equations discussed
above, the map is strictly increasing, hence having a unique fixed
point. This proves the following statement.

\begin{pro}\label{propfstyzcycle} Under the condition that
$0<\epsilon_2\ll 1, \epsilon_1=O(1)$, the asymptotic singular
orbit is a unique large $yz$-cycle.
\end{pro}

\subsection{Dynamics of Fast $Z$, Slow $Y$, and Slower $W$} The
last case we consider satisfies that the maximum reproductive
rates range from high to low in the order of $X,Z,Y,W$:
$0<\epsilon_2\ll 1\ll \epsilon_1\ll 1/\zeta$. In this case, all
nontrivial singular orbits eventually enter the $x$-slow manifold
solid ${\mathcal S}$ as before. The key difference from all
previous cases is in which direction and where the singular orbits
go next. Instead of going to either of the attracting $y$-slow
manifolds as before, singular orbits will move in the direction of
$z$ and converge to the attracting branches of the $z$-nullcline
because of the reproductive rate reversal between $y$ and $z$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{figzdynmc}
\centerline{\includegraphics[width=5in,height=2.7in]{figzdynmcab.eps}}
\ \centerline{(a)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(b)}

\centerline{\includegraphics[width=5in,height=2.7in]{figzdynmccd.eps}}
\ \centerline{(c)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(d)} \caption{In all the cases, the attracting $z$-manifolds
consist of the nontrivial $z$-nullcline, which is the vertical
plane parallel to the $w$-axis, and part of the trivial $z=0$
nullcline with $y\ge y_{\rm znl}$. (a) The case of weak $W$, $XYW$
non-competitive $Z$. (b) The case of weak $W$, $XYW$-competitive
$Z$. (c) The case of efficient $W$, $XYW$-competitive $Z$. (d) The
case of efficient $W$, $XYW$ non-competitive $Z$.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We have derived above that the nontrivial $z$-nullcline is a
vertical plane parallel to the $w$-axis in the solid ${\mathcal
S}$ as shown in Fig.\ref{figzdynmc}, through the line connecting
$(y_{\rm znl},0)$ to $(0,z_{\rm znl})$ on the $w=0$ plane. We have
also demonstrated above that to the left side of the plane $z$
increases because $y$ is relatively small, and to the right sige
of the plane $z$ decreases because $y$ is relatively large giving
$z$ a fierce competition. Also, part of the trivial $z$-nullcline
$z=0$ is also attracting, for which $y\ge y_{\rm znl}$ again
because of a relatively greater concentration of its competitor
$y$.

The principal ideas and techniques for analyzing the singular
orbits of this case are the same as the case where $y$ is faster
than $z$ as long as we replace $z$ for $y$ as the fast variable
second only to $x$ in all the arguments above.

For subcases are in order: (a) weak $w$ ($y_{\rm wnl}>y_{\rm
yfd}$) and $xyw$ non-competitive $z$ ($y_{\rm znl}<y_{\rm wnl}$);
(b) weak $w$ ($y_{\rm wnl}<y_{\rm yfd}$) and $xyw$-competitive
$z$; (c) efficient $w$ and $xyw$-competitive $z$ ($wyzw$-unstable
$xyw$-cycle); (d) efficient $w$ and $wyw$ non-competitive $z$. See
illustrations of these four cases in Fig.\ref{figzdynmc}.

Descriptions of these 4 cases are as follows. For the first case
(a), the $xyw$ equilibrium point is a global attractor for all
nontrivial singular orbits. Competitor $z$ dies out eventually, no
$xyzw$-coexistence state. For the case of (b), there is an
attracting, large $xyzw$-coexisting limit cycle. In fact, let
$\Sigma$ be the interval on the line $y=y_{\rm znl}, z=0$ as shown
with end points $w_{\rm ypd}<w_{\rm ytr}$ as shown, for which
point labelled $w_{\rm ypd}$ is the PDLS point coming from the
unstable equilibrium point $y=y_{\rm wnl},z=0,w=w_{\rm wnl}$ on
the parabola-like $y$-nullcline surface $w=q(y,z)$. Then it is
easy to see geometrically that the singular orbit of every point
from the interval $\Sigma$ returns to it. The singular induced
return map is monotone increasing. Therefore there must be a
unique fixed point of the return map similar to the case of
Fig.\ref{figchaos}(a), corresponding to an attracting limit cycle.
For the case of (c), chaotic dynamics can occur. Following the
same set-up as (b), the return map in $\Sigma$ is on longer just
monotone. As in the case of Fig.\ref{figchaos}(b), let $w_1$ be
the pre-PDLS point of the $y$-fold point $(y_{\rm yfd},0,w_{\rm
yfd})$ with respect to the fast $z$-flow. Then the singular orbit
induced return map on $\Sigma$ is qualitatively the same as that
of Fig.\ref{figchaos}(b) and (d) if $w_{\rm ypd}$ is pushed low
enough to the trivial plane $w=0$. For the last case (d), all
nontrivial singular orbits are attracted to the $xyw$-limit cycle
if it is just locally attracting with respect to the
$xyzw$-system.

In summary, if competitor $Z$ is out-reproducing its competitor
$Y$ and the respective trophic diversification hypothesis holds
for $XZ$ and $XYW$, then the only possible singular coexisting
states are either large $YZ$-limit cycles or chaotic attractors.
Co-existing equilibrium state is not possible.
