\documentclass{article}
\title{``This Town Ain't Big Enough for the Both of Us'':\\How big would be big enough? Can it be too big?}

\author{Tony H. Kim}
\date{}

\usepackage{amsmath}

\pdfpagewidth 8.5in
\pdfpageheight 11in
\setlength\topmargin{-0.25in}
\setlength\headheight{0in}
\setlength\headsep{0in}
\setlength\textheight{9in}
\setlength\textwidth{6.5in}
\setlength\oddsidemargin{0in}
\setlength\evensidemargin{0in}
\setlength\parskip{0.0in}
\setlength\parindent{0.0in}

\begin{document}
\maketitle


\setlength\parskip{0.2in}

\section{Preliminaries}
In this assignment we investigate the law of competitive exclusion. Stated in 1934 by G.F. Gause, competitive exclusion states that two competing species can coexist only if they exploit their environment differently. Our laboratory is the following dynamical system in $N_1$, $N_2$ and $P$:

$$ \begin{array}{l}
\dot{N_1} = N_1 \cdot (r - \frac{r}{k} N_1 - \frac{r}{k} N_2 - bP)  \\
\dot{N_2} = N_2 \cdot (r - \frac{r}{k} \alpha N_1 - \frac{r}{k} N_2 - (b - \epsilon)P) \\
\dot{P} = P \cdot (cbN_1 + c(b-\epsilon)N_2 - d)
\end{array} $$

where $N_1$ and $N_2$ represent the two prey populations, and $P$ the population of the predator. A brief discussion of the parameters $(r, k, \alpha, b, \epsilon, c, d)$ follows:

Clearly, the parameter $r>0$ signifies the growth rate of the prey populations as evidenced by the linear growth terms $\dot{N_i} = r N_i$. The significance of the parameter $k>0$ can be easily inferred by rewriting the equation in $\dot{N_1}$ as:
\[\displaystyle\dot{N_1} = N_1 \cdot (r - \frac{r}{k} \cdot (N_1 + N_2) - bP) \]
where it is evident that, in the absence of the predator, if $N_1 + N_2 = k$ then $\dot{N_1} = 0$. Hence, we may interpret the parameter k as the natural ``carrying capacity'' for the prey populations. As hinted by the title, the focus of this study is to observe the effects of k on the dynamics of the ($N_1$, $N_2$, $P$) system.

In the above manner we also interpret the parameter $\alpha \geq 1$. We see that $\dot{N_2} = 0$ when the term $\alpha N_1 + N_2$ equals $k$. However, since the parameter $\alpha$ magnifies the ``effective'' population of $N_1$ in smothering the growth of $N_2$, $\alpha$ gives a measure of the inherent advantage that $N_1$ has over $N_2$.

The coefficient b in the overlap terms $N_i P$ gives the rate at which $N_i$ are predated by $P$. However, in the dynamics of $N_2$ we find that the actual overlap coefficient is $(b-\epsilon)$ rather than $b$. Hence the parameter $\epsilon$ measures the degree to which $N_2$ can more successfully avoid predation in comparison to $N_1$. Since this system seeks to model the dynamics of two prey species and their common predator, we impose an upper bound on $\epsilon$, such that the term $(b-\epsilon)$ is positive. If this condition were not satisfied, we would find that $N_2$ predates upon $P$, which is a system of entirely different character.

The parameters $c$ and $d$ occur only in the differential equation for $P$ and gives the growth and death rates for the predator, respectively.

[Table summarizing the parameters]

\section{Dynamics in the absence of the predator}

We consider the cases $\alpha > 1$ and $\alpha = 1$ separately.

\subsection{The funnel trap ($\alpha > 1$)}

By initializing the system at $P_0 = 0$ we obtain the reduced system:

$$ \begin{array}{l}
\dot{N_1} = N_1 \cdot (r - \frac{r}{k} N_1 - \frac{r}{k} N_2)  \\
\dot{N_2} = N_2 \cdot (r - \frac{r}{k} \alpha N_1 - \frac{r}{k} N_2) \\
\end{array} $$

We can visualize this system and its nullclines on the plane:

[Enter nullcline equations]

[Insert diagram here]

As can be seen, there are only three fixed point solutions which correspond to:

\begin{itemize}
\item $(N_1,N_2)^{*} = (0,0)$: Mutual extinction.
\item $(N_1,N_2)^{*} = (k,0)$: $N_1$ drives $N_2$ into extinction.
\item $(N_1,N_2)^{*} = (0,k)$: $N_2$ drives $N_1$ into extinction.
\end{itemize}

This result immediately shows that in the reduced $N_1$, $N_2$ system there is no fixed point that allows for mutual coexistence for $N_1$ and $N_2$. Thus, for this model, we have verified Gause's comment that when a single species posesses a sole advantage over its competitor, there cannot be a steady state solution involving \emph{both} $N_1$ and $N_2$.

As the $\alpha > 1$ case confers on $N_1$ competitive advantage over $N_2$, we may expect that the only \emph{stable} solution among the three is $(N_1, N_2)^* = (k,0)$. Linear analysis is little help in this case, however, as the two non-origin fixed points reveal to be marginal with the Jacobian matrices:
\begin{equation}
J^{*} = \left(
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,1} & a_{2,3} \\
a_{3,1} & a_{3,1} & a_{3,3} \\
\end{array}
\right)
\label{}
\end{equation}

for $(k,0)$ and $(0,k)$ respectively. We are luckier at the origin, however, where the linearized stability turns out to be:
\begin{equation}
\left(
\begin{array}{c}
\dot{\delta N_1} \\
\dot{\delta N_2} \\
\end{array}
\right) = \left(
\begin{array}{cc}
r & 0 \\
0 & r \\
\end{array}
\right) \cdot \left(
\begin{array}{c}
\delta N_1 \\
\delta N_2 \\
\end{array}
\right)
\label{}
\end{equation}

which shows that near the origin the components of the displacement vector $\delta N_1$ and $\delta N_2$ are decoupled and are unstable in both in the so-called ``star'' geometry.

For the stability of the remaining fixed points, we consider the phase plane instead from a geometric perspective. Taking into account the signs of $\dot{N_1}$ and $\dot{N_2}$ in each region of the plane defined by the nullclines gives:

[Insert diagram]

In each region and along each nullcline, we expect the direction of the flow to be as described in the representative vectors above. We can immediately interpret the space between the two sloped nullclines as a funnel that traps all trajectories that begin in the first quadrant and leads them to the fixed point at $(k,0)$. To see this, consider a trajectory that begins in the region marked A in the above figure. The trajectory is continuously pushed to the right until it enters the trap and finds itself unable to cross the nullcline $N_2 = k - N_1$; and since everywhere in the trapping region $\dot{N_1}$ is positive and $\dot{N_2}$ is negative, this trajectory must approach the fixed point at $(N_1, N_2)^{*} = (k,0)$ asymptotically.

Likewise, a trajectory that begins in region B is initially attracted to the trapping region but finds itself unable to escape past the nullcline $N_2 = k - \alpha N_1$, and in time dwindles to $(N_1, N_2)^{*} = (k,0)$ as well. Representative cases of these two trajectories are shown below.

[Insert diagram]

The above argument shows that for the $N_1$, $N_2$ system in the absence of predators, there is absolutely \emph{no} possible coexistence between $N_1$ and $N_2$. Firstly, we find that variations in r have no effect in the qualitative geometric features of the $N_1,N_2$ phase space. Furthermore, variations in k and $\alpha > 1$ within their parameter ranges only stretch and shrink the funnel \emph{without} changing its general behavior. In this way $N_1$ is quite the unpleasant neighbor, unwilling to share living space under any conditions. 

\framebox[6.5in][c]{\parbox[t]{5.5in}{In conclusion, in the absence of the predator, the only possible long term solution is $(N_1, N_2)^{*} = (k,0)$, i.e. complete domination by $N_1$, \emph{regardless} of the system parameters, including system size.}}

\subsection{Symmetric case ($\alpha = 1$)}

For completeness we consider now the (rather uninteresting) case when $\alpha = 1$, where the governing equations become symmetric:

$$ \begin{array}{l}
\dot{N_1} = N_1 \cdot (r - \frac{r}{k} N_1 - \frac{r}{k} N_2)  \\
\dot{N_2} = N_2 \cdot (r - \frac{r}{k} N_1 - \frac{r}{k} N_2) \\
\end{array} $$

This scenario corresponds to the previous picture when the funnel trap has collapsed on itself. 

By considering the sum $N_s = N_1 + N_2$, we obtain from the above system a decoupled logistic equation in $N_s$, which can be easily solved analytically. Using this result, then, gives a separable equation in $\dot{N_s}$ where $N_s = N_1 - N_2$. Hence, the symmetric case is amenable to a full analytic solution for all $t$. However, the ultimate behavior of this special case is rather apparent without much work: trajectories that begin in the first quadrant will land somewhere along the superposition of the two nullclines at $N_1 + N_2 = k$. Clearly, linear analysis at any of such solutions will give a ``comb'' stability.

\subsection{Introduction of the predator: The trap fails!}

Postponing the theoretical development until the following section, we can simply visualize the dynamics of the full $(N_1,N_2,P)$ system in terms of the $N_1,N_2$ projections. For the diagram below, I have fixed $k$, i.e. the capacity of the system, at 300 units. The initial state of the system is: $(N_1,N_2,P) = (80, 80, p)$ where $p=P_0$ is the parameter being varied.

[Need to indicate that $p$ represents a different advantage that could potentially bring $N_1$ and $N_2$ into coexistence.]

\section{Full $N_1,N_2,P$ dynamics}

Analysis of the three-dimensional system turns out to be algebraically cumberson. Still, the standard techniques of systems analysis can be applied rather procedurally.

\subsection{Nullcline analysis}

By inspection, the dynamical system,

$$ \begin{array}{l}
\dot{N_1} = N_1 \cdot (r - \frac{r}{k} N_1 - \frac{r}{k} N_2 - bP)  \\
\dot{N_2} = N_2 \cdot (r - \frac{r}{k} \alpha N_1 - \frac{r}{k} N_2 - (b - \epsilon)P) \\
\dot{P} = P \cdot (cbN_1 + c(b-\epsilon)N_2 - d)
\end{array} $$

yields the following nullplanes:

[Enter nullplanes here]

If we select a 3-tuple of nullplanes, one from each pair, we obtain a system of equations that can be solved for the fixed points of the full system. As it turns out, there are only 6 unique solutions from such combinations. The first three are of little interest to us now, as we've encountered them already in the previous section when studying the evolution of $N_1$ and $N_2$ in the absence of $P$. They are (again):

\begin{itemize}
\item $(N_1,N_2,P)^{*} = (0,0,0)$: Mass extinction.
\item $(N_1,N_2,P)^{*} = (k,0,0)$: $N_1$ drives $N_2$ into extinction.
\item $(N_1,N_2,P)^{*} = (0,k,0)$: $N_2$ drives $N_1$ into extinction.
\end{itemize}

The following two are only slightly more interesting. They represent the interaction of \emph{one} of the prey species with the predator:

[Insert equations here]

Finally, the fixed point solution that will occupy the remainder of the assignment is the one arising from the three nontrivial nullplanes of the system, i.e. the fixed point that solves the following equation:

\begin{equation}
J^{*} = \left(
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,1} & a_{2,3} \\
a_{3,1} & a_{3,1} & a_{3,3} \\
\end{array}
\right)
\label{}
\end{equation}

For example, for the standard set of parameters and $k = 400$, this fixed point turns out to be $(N_1,N_2,P)^{*} \approx (189.47,105.26,26.31)$ and appears to be stable. [Reference to figure here]. The reason for our particular interest in this solution is because this fixed point represents the coexistence solution that was outlawed in the $(N_1,N_2)$ case by Gause's principle. In the full dynamics of $(N_1,N_2,P)$ we find that there \emph{is} in fact a possible coexistence between the prey. In retrospect, this is not so surprising (and especially considering the specific wording of Gause's principle in the introduction) since $\alpha$ represents $N_1$'s advantage, while $\epsilon$ is an advantage that belongs to $N_2$. Obviously, the latter has no effect in the absence of the predator.

However when the predator is present, $N_1$ and $N_2$ both have a different advantage over each other. It seems plausible, then, this scenario represents a more `fair' competition between the two prey species with, possibly, much richer solutions.

[Enter figure here]

Further justification for our special attention to this particular fixed point comes from the fact that more complicated attractors that arise for higher values of k develop from the instabilities at this location. This proposition will be quantified in the upcoming sections.

\subsection{How big is big enough?}

We now answer the question that opened this assigment. Given that there is now the possibility of mutual coexistence, does the size of the `town' (i.e. $k$) determine whether the two species can coexist? 

We must first establish the range of $k$ for which the fixed point under consideration gives a possible solution to our problem. In the context of populations, all three quantities must be nonnegative. While a symbolic treatment of the solution of [equation reference] turns out to be much too cumbersome, the fixed point is in principle given by the equation:

\begin{equation}
J^{*} = \left(
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,1} & a_{2,3} \\
a_{3,1} & a_{3,1} & a_{3,3} \\
\end{array}
\right)
\label{}
\end{equation}

We then evaluate the curve traced by [equation reference] below.

[Figure]

Evidently, for small values of $k$, $N_2^{*}$ is negative and this fixed point does not give a realistic solution to our population problem. Plotting $N_2$ against $k$ reveals a surprisingly simple relationship between $N_2^{*}$ and $k$. [Figure reference] (This finding motivated the general solution to this problem that follows). Solving this linear relationship in our particular parameter set reveals $k_{crit} \approx 311.5$. The town must be approximately 312 `units' large for the both of us!

[Figure]

\subsubsection{General case}

Motivated by the apparent linearity of $N_2^{*} = f(k)$ in the above, I have considered [equation reference] symbolically. This yields:

\begin{equation}
N_2^{*} = f(k) = \frac{\epsilon b}{\epsilon^2-b (\alpha-1) (b-\epsilon)} k - \frac{b(\alpha-1) + \epsilon}{dc\left[\epsilon^2-b(\alpha - 1)(b - \epsilon)\right]}
\label{}
\end{equation} 

For now, let us take the existence of a nonnegative $N_2$ as the sole criterion for the existence of a valid nontrival fixed point. (Note that we have yet to discuss stability at all!) We'll also assume that our parameter set is such that the intercept of $f$ is negative, i.e. that the constant term in [equation reference] is negative. Since the numerator (and $d$ and $c$) is always positive given our parameter ranges, this requires:

\begin{equation}
M = \epsilon^2 - b(\alpha-1)(b-\epsilon) > 0 
\label{}
\end{equation} 

However, the above condition is not generally true. For the particular parameter values of this investigation, $M = 7.6 \times 10^{-5} > 0$ and hence $N_2$ becomes nonnegative for some $k > k_{crit}$. Suppose instead that $\epsilon = 0.001$; in that case, we obtain $M = -4.4 \times 10^{-5} < 0$. This possibility, along with [equation reference] suggests that for a different set of parameters we may find coexistence for \emph{smaller} values of $k < k_{crit}$, as $N_2$ begins nonnegative for small $k$ and becomes negative for larger $k$. While this possibility is interesting, it will not be considered further. Instead, we will simply assume that in our systems, $N_2 > 0$ for $k > k_{crit}$.

We then solve for $k_{crit}$:

\begin{equation}
k_{crit} = \frac{ab-(b-\epsilon)}{bcd\epsilon}
\label{}
\end{equation}

The above equation gives the critical value of $k$ above which the system produces a valid nontrivial fixed point. Unfortunately, it is not so useful due to the additional assumptions that we have made in the analysis. The finding of the parameters $c$ and $d$ in the expression for $k_{crit}$ is unexpected, however, indicating a lower $k_{crit}$ as the parameters controling the flow of the predator population is increased. 

\subsubsection{Interpretation}

In this subsection I'd like to consier the curve that is parameterized through $k$. One approach to this question would be to very carefully tread through [equation reference] which, despite its algebraic complexity, is linear and hence fully solveable. This is the approach taken in the previous subsection. Having completed such an exercise, one could simply point to the resulting function $N_2^{*} = f(k)$ and report: ``Aha! The reason for the particular behavior of so-and-so is due to the so-and-so functional form of $f$!'' 

[More stuff goes here]

\subsection{Stability}

In the previous sections, attention was paid to the movement of the nontrivial fixed point into the first quadrant of phase space. We now ask about its linear stability. To answer this question, consider again the system:

$$ \begin{array}{l}
\dot{N_1} = r N_1 - \frac{r}{k} N_1^2 - \frac{r}{k} N_1 N_2 - b N_1 P  \\
\dot{N_2} = r N_2 - \frac{r}{k} \alpha N_1 N_2 - \frac{r}{k} N_2^2 - (b - \epsilon)N_2 P \\
\dot{P} = cbN_1 P + c(b-\epsilon)N_2 P - d P
\end{array} $$

and its Jacobian matrix:

\begin{equation}
J^{*} = \left(
\begin{array}{ccc}
r - 2 \frac{r}{k}N_1^{*} - \frac{r}{k}N_2^{*}-bP^{*}& \frac{-r}{k}N_1^{*} & -b N_1^{*} \\
\frac{-r}{k} \alpha N_2^{*} & r-2\frac{r}{k}N_2^{*}-\frac{r}{k}\alpha N_1^{*}-(b-\epsilon) P^{*} & -(b-\epsilon)N_2^{*} \\
cbP^{*} & c(b-\epsilon)P^{*} & cbN_1^{*} + c(b-\epsilon)N_2^{*}-d \\
\end{array}
\right)
\label{}
\end{equation}

However, as $(N_1,N_2,P)^{*}$ solves [equations reference], we simplify the above with:

\begin{itemize}
\item $\frac{r}{k}N_1^{*}+\frac{r}{k}N_2^{*}+bP^{*} = r$
\item $\frac{r}{k}\alpha N_1^{*}+\frac{r}{k}N_2^{*}+(b-\epsilon )P^{*} = r$
\item $cbN_1^{*} + c(b-\epsilon)N_2^{*} = d$
\end{itemize}

and obtain:

\begin{equation}
J^{*} = \left(
\begin{array}{ccc}
-\displaystyle\frac{r}{k}N_1^{*}& \displaystyle\frac{-r}{k}N_1^{*} & -b N_1^{*} \\
\displaystyle\frac{-r}{k} \alpha N_2^{*} & -\displaystyle\frac{r}{k}N_2^{*} & -(b-\epsilon)N_2^{*} \\
cbP^{*} & c(b-\epsilon)P^{*} & 0 \\
\end{array}
\right)
\label{}
\end{equation}

As usual, the characteristic equation is given by 

\begin{equation}
det(J^{*}-\lambda I_3) = 0
\label{}
\end{equation}

Since $(N_1,N_2,P)^{*}$ is parameterized through $k$, in principle the [equation reference] represents a cubic-function whose coefficients are parametrized by $k$. Again, a symbolic treatment is far too cumbersome, so we resort to a numerical approach as described in the following section.

\subsubsection{Approach}

As always, the quantities of interest are the real parts of the eigenvalues, which are determined by [equation reference]. We will proceed as follows:

\begin{enumerate}
\item Consider $det(J^{*}-\lambda I_3) = f(\lambda)$ over the real line. Since the coefficients of the cubic polynomial are positive, we must have at least one real root. Obtain this root, $\lambda_R$, by interpolation.
\item Numerically factor out the term $(\lambda - \lambda_R)$. For the resulting quadratic, interpolate the real parts of the complex roots by calculating the vertex of the parabola. 

(Let $\lambda_{C^{\pm}}$ represent the complex roots. By the quadratic formula, we have $Re\left[\lambda_C\right] = \left[\lambda_{C^{\pm}}\right] = -\frac{b}{2a}$ where $a$ and $b$ are the coefficients associated with the quadratic and linear terms, respectively. From basic calculus, we have that the vertex of a parabola is given by $-\frac{b}{2a}$. Hence:

\begin{equation}
Re\left[\lambda_C\right] = \lambda_{V}
\label{}
\end{equation}

where $\lambda_{V}$ is the vertex of the quadratic.

\end{enumerate}

A typical plot resulting from the above procedure is pictured in [figure reference]. The computational details and the scripts used in the calculation are reported in the appendix of this assignment. Here, we merely state some of the results:

\begin{tabular}[t]{|l||r|r|}\hline
k 		&  $\lambda_R$ 	    & $Re(\lambda_C)$   \\\hline\hline
300 	&	 0.03051893166903 & -0.32887555750354 \\\hline
310 	&  0.00319373625057 & -0.32228363414204 \\\hline
311.1 &  0.00001201063378 & -0.32168591081887 \\\hline
311.11111... & -0.00002013337241 & -0.32178722300964 \\\hline
312   & -0.00250920098777 & -0.32028260248924 \\\hline
350   & -0.16946825476467 & -0.26019842172357 \\\hline
375   & -0.42103166966124 & -0.14692910730032 \\\hline
400   & -0.55943779280922 & -0.08869650680063 \\\hline
425   & -0.64098233317196 & -0.05759456829995 \\\hline
450   & -0.69811479991791 & -0.03763337625298 \\\hline
475   & -0.74126341079928 & -0.02374935563374 \\\hline
500   & -0.77530885078003 & -0.01365102428031 \\\hline
525   & -0.80293799603971 & -0.00610204281474 \\\hline
550   & -0.82584020472873 & -0.00035595010773 \\\hline
551   & -0.82666937586320 & -0.00015662851995 \\\hline
551.79463 & -0.82732067328578 & -0.00000000095928 \\\hline
552   & -0.82748872440641 &  0.00004029256492 \\\hline
553   & -0.82830539969215 &  0.00023540570821 \\\hline
554   & -0.82911941270112 &  0.00042873020269 \\\hline
555   & -0.82993077435794 &  0.00062028508698 \\\hline
560   & -0.83391442890901 &  0.00154948535598 \\\hline
575   & -0.84507461106758 &  0.00406453131344 \\\hline
\end{tabular} 

From the previous section, we have $k_{crit} = \frac{2800}{9} \approx 311.11111...$. Evaluation of $\lambda_R$ and $Re(\lambda_C)$ at $k_{crit}$ gives negative values, which indicates that the nontrivial fixed point considered in the previous sections is \emph{stable} as soon as it appears in the first quadrant. As we continue to increase $k$, we find that $\lambda_R$ remains negative, while $Re(\lambda_C)$ becomes positive for $k^{*} \approx 551.79463$. 

\framebox[6.5in][c]{\parbox[t]{5.5in}{For our parameter set, we have a linearly stable fixed point representing coexistence of $N_1$ \emph{and} $N_2$ while the capacity of the system ranges in: $311.11111... < k < 551.79463$}}

\section{More complicated dynamics}

\subsection{Limit cycle}

As $k \rightarrow k^{*}$ the rate of approach to the fixed point slows. As the capacity of the system is increased even further, we develop a small limit cycle about the fixed point. In [reference to figure] we show the time series in $N_1$ for parameter values of 540, 550, and 560. 

[Enter figure here]

And below, we see the development of the limit cycle in phase space ($k = 560$ pictured):

[Enter figure here]

\subsection{Period doubling}

Given the small limit cycle about the fixed point the interesting question now is to ask how does the cycle's behavior change for further increases in $k$? Given the growing instability [reference to table] of the fixed point for higher k, it is natural to believe that the size of the circuit will expand. However, another familiar behavior of such cycles is the `period doubling' in which the trajectory `completes multiple turns` before closing on itself. This effect is best characterized by a Poincare section of the attractor.

\subsection{Chaos}

\subsubsection{Lyapunov exponent}

\subsubsection{Fractal dimension}

\subsection{Even larger k}

\section{Application: Habitat destruction}

\end{document}