% SIMIODETemplate-TeacherORStudent.tex and simiode.cls for LaTeX2e
% Version 29 July 2015.
% 
% This class (simiode.cls) is a modification of that used by 
% the journal PRIMUS - Problems, Resources, and Issues in 
% Mathematics Undergraduate Studies. 
% 
% See http://www.tandfonline.com/loi/upri .
% 
% The class was developed by  Marc Frantz with suggestions of the 
% PRIMUS Editorial Board, especially John Thoo.
%
% This template is to be used with the  class file "simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
%
\documentclass{simiode}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Basically the Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review (but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
% This option was created  by Eric Sullivan, Carroll College, Helena MT.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{hyperref}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding double spacing %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manuscripts should be submitted doubled spaced.  This is 
% accomplished via the \baselinestretch command below.
% However, you may obtain an idea of how your article will 
% look in final form, so that you can better plan your 
% presentation.  To see your manuscript formatted very 
% similarly to the on-line published format of  SIMIODE, 
% comment the following line out by inserting a percent sign % (%) at its beginning.

\renewcommand{\baselinestretch}{1.3}

\newcommand{\pwrt}[2]{$\frac{\partial #1}{\partial #2}$} % Macro for partials

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding macros %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Top matter is above and 
% the beginning of the document below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\newpage
% You only need to change the NAME OF SCENARIO in the next line.
\title{\ifnum\TeacherEdition=1 TEACHER VERSION \else STUDENT VERSION \fi\\Algae Population Self-Replenishment}


\markboth{Algae Population Self-Replenishment}{Algae Population Self-Replenishment}



\author{Richard Corban Harwood\\

George Fox University\\
Newberg  OR 97132USA\\
rharwood@georgefox.edu
}

\makeStitlePDFLaTex


\begin{abstract}
This modeling scenario investigates the massive algal blooms that struck Lake Chapala, Mexico, in 1994. After reading a summary of news articles on the incident, students create an ODE system model from a verbal description of the factors, visualize this system using an executable Java applet (PPLANE) to predict overall behavior, and then analyze the nonlinear system using the Jacobian matrix, eigenvalues, phase plane, and feasibility conditions on parameters to fully describe the system behavior. Students are expected to be familiar with systems of differential equations, equilibria, jacobian matrices, and eigenvalues. Students will learn modeling from qualitative descriptions, nondimensionalization, applying feasibility conditions to parameters, and how to use technology to interactively analyze a system of differential equations.
\end{abstract}


\section*{SCENARIO DESCRIPTION}

\subsection*{Background Reading}  %%%%%%%%%%%%%%%%%%%
Lake Chapala, Mexico, hosts a diverse throng of fish and migratory birds as well as numerous fishing and farming communities \cite{Tony}. The Lerma Rivers system, the ecological and hydrological core of central Mexico, drains into Lake Chapala before draining out to the Pacific Ocean. As such, problems in Lake Chapala are indicators of the health of the entire ecosystem.

%Add explanation of other factors (from Tony's article)
The most visible ailment of Lake Chapala is dehydration which has several side effects including algal blooms. This dehydration is due to a combination of evaporation, irrigation of it's headwaters to area farms, drainage by the Santiago River to the Pacific Ocean, and piping of potable water to Guadalajara's 4 million people. The lake has come close to drying up completely a couple of times in the last century and currently sits half-full with an average depth of 4.5 meters ($\approx 15$ ft). In addition, the lake's once pristine waters are now muddy, with transparency rarely exceeding 30 centimeters ($\approx 1$ ft). As a combined result of the dehydration and increased farming in the area, Lake Chapala has high levels of phosphorus and other sediments. Finally, the lake has been frequented by algal blooms, rapid increases in algae population, which have been irritations to the tourists and communities in the surrounding area \cite{Tony}.

Of the many factors that threaten Lake Chapala, we focus our attention on the algae, specifically referring to the free-floating phytoplankton. These algae are tiny drifting plants that grow and reproduce using photosynthesis.  The first evidence of algal blooms in Lake Chapala was in 1994. Since then algal blooms have freqented the lake and the growth of seaweed, another form of algae, has hampered the use of motorized fisher boats \cite{Tony}. Since the initial appearance of algae in 1994 in Lake Chapala, the blooms, though non-toxic, have at times tainted the smell and taste of potable water for the citizens of Guadalajara \cite{Tony}. Currently no toxic strains of algae have been found in the lake, however, questions still linger concerning how the current strains of algae arrived and how the danger of toxic strains may be avoided \cite{Lopez}. 

A Citizens Report to the Commission for Environmental Cooperation and further investigation purport that the Mexican government is not complying to protect the water quality and ecology of the Lerma Rivers system and is lax in their enforcement of the environmental laws concerning the region \cite{CEC}. The Mexican goverment has assured environmental law enforcement with its compliance. Specifically, the National Water Program of Mexico has defended its actions in restoring and preserving the water quality for the Lerma Rivers system \cite{CEC:response}. Algae is just one factor in this environmental contention, but is highly visible due to its effect on the taste and smell of drinking water in Guadalajara. Thus it has been at the forefront in media reportings of the crisis of Lake Chapala \cite{Tony}.

Considering the government efforts to eliminate the algae and the media reportings of algal blooms of increasing freqency and magnitude \cite{Lopez}, the resilience of the algae in Lake Chapala cannot be neglected. Let us consider the possibility that algal blooms are being eliminated near the surface, but dormant algae deeper in the lake is unaffected. Churning in the water due to fishing motorboats brings this dormant population to the surface where they are activated by access to sunlight. The increase in agricultural chemical residues, heavy metals, and dissolved solids provide a rich food base for the algae to rapidly grow \cite{Tony}. 

Reading questions:
\begin{enumerate}
\item What are some of the problems that algal blooms cause for the communities around Lake Chapala?
\item How does water depth impact a free-floating algae?
\end{enumerate}

\subsection*{Modeling Assumptions}%%%%%%%%%%%%%%%%%%%
You will develop a sunlight-dependent model to describe the resilient ability of an algae population to bounce back after near-annihilation.

The modeling process starts with applying common simplifying assumptions, such as spatial uniformity. Assuming that the density of algae is uniform horizontally, we reduce population change to time and depth. Assuming a threshold depth, after which photosynthesis is no longer possible, we can discretize the depth variability into two regions or compartments with spatially uniform population densities. Thus, we are left with an upper population density of active algae $u(t)$ which are photosynthesizing and a lower dormant population $v(t)$ with units $g/m^3$ or other volumetric density units, splitting the lake into two domains represented by points where these population densities reside but still vary in time. Time, $t$, is usually measured in days to average out the daily cycle of sun intensity and focus on the scale of days to months for population growth. Volumetric population densities are used instead of population counts so that the model can be applied to any volume of water tested. 

%We start by splitting the water body, a vertical cross-section of the water, into two regions dependent upon the amount of sunlight received: upper, where photosynthesis is possible, and lower, where it is not. Assuming that the density of algae is uniform horizontally and can be discretized vertically into two regions, we can represent the algae populations as time-varying densities of active $u(t)$ and dormant $v(t)$ algae in the upper and lower region with units $g/m^3$ or other volumetric density units. Time, $t$, is usually measured in days to average out the daily cycle of sun intensity and focus on growth on the scale of days to months. Volumetric population densities are used instead of population counts so that the model can be applied to any volume of water tested. 

Since algae depend on sunlight for photosynthesis, only the $u$ population is capable of growth \cite{Encyclopedia}. Due to circulation in the water (based on river inflow and watercraft propulsion), these densities mix with one another. As the algae circulate with the water, some are lost by either sinking to the bottom, dying, or being eaten. Note, this loss term provides only a first order approximation of the effect of predation. Expanding the system to include a predator term may mitigate significant errors when comparing to subsequent data. In essence, $v$ is a perishable storage population that $u$ circulates into and partially receives back from. This gives the following system:
\begin{eqnarray*}
\frac{du}{dt}&=&\textrm{Growth -  Circulation Down +  Circulation Up},\label{verbal_model}\\ 
\frac{dv}{dt}&=&\textrm{Circulation Down - Circulation Up - Loss}.\nonumber
\end{eqnarray*}

\subsection*{Using PPLANE Software}%%%%%%%%%%%%%%%%%
A powerful tool for visualizing and interactively analyzing systems of differential equations is PPLANE, an executable Java applet  exported from MATLAB source code which can be run without having MATLAB, but it does require \href{http://www.java.com/en/download/help/download_options.xml}{Java Runtime Environment} installed \cite{Polking1}. A short tutorial with example image follow below. For a more extensive tutorial, see \cite{Polking}. PPLANE 2005.10 was updated in 2015 by Dr. Joel Castellanos at the University of New Mexico from code originally developed in 1994 by Dr. John Polking at Rice University. The file pplane.jar is available to download at \href{http://math.rice.edu/~dfield/dfpp.html}{http://math.rice.edu/{$\sim$}dfield/dfpp.html}. Documentation and original MATLAB source code is also availabe at the site.



\begin{figure}[t]
 \centerline{\includegraphics[width=4.7in]{6-20-1-PPLANEdisplayex2}} 
\vspace*{-7.8 true cm}
\begin{center}
\parbox{9 true cm}{\caption{\label{fig:6-20-1-PPLANEdisplayex2} Phase plane representation of circulation between two compartments with $a=\frac12,b=\frac13$, with several plotted solution curves and a graph of an example function.}}
\end{center}



\end{figure}

Figure \ref{fig:6-20-1-PPLANEdisplayex2} shows an example display of PPLANE with a system similar to model (\ref{verbal_model}), but focusing solely on circulation between the two compartments:
\begin{eqnarray*}
\frac{du}{dt}&=&-au+bv,\nonumber\\ 
\frac{dv}{dt}&=&au-bv.
\end{eqnarray*}

Here is a quick tutorial for first time users of PPLANE. Download and open the PPLANE.jar file. The copyright page gives a brief summary of what the PPLANE software can do; click ``Ok". In the PPLANE equation window, enter your system by typing the dependent variables $u,v$ into the top left boxes and their derivative functions in the boxes to their right. Note that multiplied variables can be written sequentially with or without an `*' between them, unlike the syntax of the underlying MATLAB code. Parameter values ($a=1/2$ and $b=1/3$ for example) need to be defined in the boxes in the lower left and the window size dimensions (min $x = 0$, max $x = 1$, min $y = 0$, max$ y = 1$) go in the lower right. Note, the `$x$' and `$y$' here refer to the first and second dependent variables. Clicking on ``Graph Phase Plane" will produce the phase plane where you can click multiple times on the screen to see curves solved numerically using the clicked position as the initial position $u(0), v(0)$. To see the individual functions of time, select the ``Graph" header, ``both $x-t$ \& $y-t$", and then click a point on a solution curve you have already plotted. To clean up the graph of the individual function to a more reasonable interval, you can click the vertical bar following the cursor twice for lower and upper bound and click ``Crop" to reduce the graph.


\section*{ACTIVITIES}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\begin{enumerate}
\item {\bf Set Up Model:} Fill in the verbal model with your own functions of $u,v$ for each slope function $u'=f_1(u,v)$ and $v'=f_2(u,v)$. Start simple and explain your selections to someone else for feedback. Nondimensionalize your model by scaling and combining parameters to reduce the number of parameters in the model. 

Hint: One simple method is to scale everything by a common factor and make substitutions, such as turning $\frac{du}{dt}=bu(1-au)+bv$ into $\frac{dx}{d\tau}=x(1-x)+y$, where $x=au, y=av, \tau = bt$. Note, it may be helpful to rename all variables and parameters in a common way. Note, in this example, all parameters have been built into the new variables and thus will impact the solution behavior indirectly based upon bounds on $x$ and $y$.

\item {\bf Predict Behavior:} Type your model into PPLANE (see tutorial and Figure \ref{6-20-1-PPLANEdisplayex2}) and choose parameter values that seem feasible to you based upon conditions on the original parameters (a good initial guess). Check the plot and describe the behavior. Are solutions converging (stable or cyclic) or diverging (unstable)? Are solutions moving in exponential or straight paths (real eigenvalues) or spiraling (nonreal eigenvalues)? Alter the parameter values. The solution behaviors should shift, but does the type of behavior (different answers the the above questions) change?

\item {\bf Linearize Model:} Using your nondimensionalized version of your model, find equilibirum points by setting your slope functions equal to zero and solving for values of the variables satisfying both equations. Determine the feasibility conditions which need to be placed upon your new variables and the original variables (reversing your substitutions) to have nonnegative equilibrium points. These represent conditions for biologically feasible equilibrium populations.

\item {\bf Analyze Jacobian:} Then form the Jacobian matrix at each point and analyze the solution behavior based upon the trace and determinate of each Jacobian matrix.

Hint: Using $\frac{dx}{d\tau}=f_1(x,y), \frac{dy}{d\tau}=f_2(x,y)$ as an example nondimensionalized system, the Jacobian matrix would be 
$$
J(x,y)=\left(
\begin{tabular}{cc}
\pwrt{f_1}{x} & \pwrt{f_1}{y} \\
\pwrt{f_2}{x} & \pwrt{f_2}{y} \\
\end{tabular}
\right)
.$$

\item {\bf Analyze Parameters:} Summarize feasibility (conditions for nonnegative equilibria) and stability conditions (conditions which make the nonnegative equilibria a stable node) on parameters. Use PPLANE to explore this parameter space and describe the amount of increasing/decreasing for each parameter upon the scaled and original variables. Demonstrate an example of each type of behavior for feasible equilibria.
\end{enumerate}

%\begin{thebibliography}{98}

%\bibitem{Harwood2008}V. S. Manoranjan, M. A. O. Gomez, and  R. C. Harwood, ``Modeling Algae Self-Replenishment,'' {\it Journal of Interdisciplinary Mathematics}, Vol. 11, pp. 681-694, 2008.
%\bibitem{Polking1} John Polking and Joel Castellanos. {\it PPLANE}. Software available at  \href{http://math.rice.edu/~dfield/dfpp.html}{http://math.rice.edu/{$\sim$}dfield/dfpp.html}. Accessed on 17 July 2015.

%\end{thebibliography}


\teacher{
\section*{NOTES FOR TEACHERS}

\section*{SOLUTION}
The following is an example solution from \cite{Harwood2008} which leads from the given assumptions and works through the outlined steps. See the discussion afterwards for ways to encourage and handle diversity in your students solutions.

\subsection*{Set Up Model}
Assuming a sufficient amount of sunlight and nutrients, as is the case at Lake Chapala, we can start with a constant proportion of growth rate to population, befitting exponential growth. Yet, since this model is sunlight-dependent, the growth of algae in the upper region will be heirarchically inhibited as top layers of algae cast shadows on those below. One way to model this restricted growth is Growth $=\frac{\alpha u}{1+\beta u}$, where $\alpha$ is the growth proportionality constant and $\beta$ is the constant restriction rate. Note that this is similar to the logistic growth model by geometric series if $|-\beta u|<1$ since
$$ \frac{\alpha u}{1-(-\beta u)} = \alpha u \sum_{k=0}^\infty (-\beta u)^k \approx \alpha u (1 - \beta u).
$$
Though dehydration of the lake changes the volume of the water in each region, we do assume that they remain constant for the duration of the simulation (since the volume varies by year whereas the algae vary more by day). Thus, ignoring loss, the circulation between the two regions, though not equal, will be a constant percentage of the current algae density in the region of origin: Circulation Down $=\gamma u$ and Circulation Up $=\delta v$. Since the loss is also a percentage of the lower population, it can be combined with the upward circulation as Circulation Up + Loss $=\lambda v$. The model system is then
\begin{eqnarray}
\frac{du}{dt} &=& \frac{\alpha u}{1+\beta u}-\gamma u+\delta v,\nonumber\\
\frac{du}{dt} &=& \gamma u-\lambda v
\label{original}
\end{eqnarray}
where all parameters are constant, nonnegative, $\delta\leq \lambda$ due to some loss in the lower region, and $\gamma, \delta, \lambda \leq 1$, since they are percentage rates.

Seeking a simpler form of this model, we look for ways to group parameters and apply the substitution transformations:
\begin{eqnarray}
x=\beta u,\hspace{.5in} y=\frac{\beta \lambda}{\gamma} v,\hspace{.5in} \tau =\gamma t
\end{eqnarray}
The system, now with only three parameters from the original five, becomes
\begin{eqnarray}
\frac{dx}{d\tau} =& \frac{A x}{1+x}-x+B y &= f_1(x,y),\nonumber\\
\frac{dy}{d\tau} =& C x-C y &= f_2(x,y),
\label{scaled}
\end{eqnarray}
where $A=\alpha/\gamma$, $B=\delta/\lambda$, and $C=\lambda/\gamma$. Note, by construction, $A,B,C$ are all nonnegative, and to circumvent trivial solutions, we will assume they are all positive. Also, using the chain rule, $\frac{du}{dt}=\frac{du}{d\tau}\frac{d\tau}{dt}$ we can transform the independent time variable to compensate for the extra $\frac{1}{\gamma}$ on the derivative terms.

\subsection*{Predict Behavior}

\begin{figure}[h!]
 \centerline{\includegraphics[width=4.0in]{6-20-2-pplanealgae}} 
\begin{center}
\vspace*{-6.8 true cm}
\parbox{11 true cm}{\caption{\label{fig:6-20-2-pplanealgae} Example phase plane of the model(\ref{scaled}) with $A=\frac34, B=\frac{1}{2}, C=\frac14$, resulting from original parameters $\alpha=\frac35, \beta=\frac12, \gamma=\frac45, \lambda=\frac15, \delta=\frac{1}{10}$. Several solutions (blue curves) are shown converging on the coexistence equilibrium point (red point), along with the two  (yellow and red curves).}}
\end{center}
\end{figure}

To get a feel for the behavior of the system (\ref{scaled}), we typed the slope functions in to PPLANE, chose example parameter values for $\alpha, \beta, \gamma, \delta, \lambda$ between 0 and 1 to satisfy the original conditions and computed values for $A,B,C$. See Figure~\ref{fig:6-20-1-PPLANEdisplayex2} for placement of slope functions and parameter values in the PPLANE applet. Figure~\ref{fig:6-20-2-pplanealgae} enumerates these example parameter values as well as demonstrating an example phase plane, several plotted solutions, and a visualization of the  which we will be using to analyze the model. The  can be plotted by selecting ``Plot " under the ``Solution" tab. 

From this snap shot of the system behavior we see solutions in the first quadrant that are converging to an equilibrium point and are moving in exponential (nonspiral) paths. From altering the parameters a bit, we continue to see similar behavior and predict that the system has a real stable sink in the first quadrant which represents a point of coexistence which initial positive populations will converge to. Next, we analyze the model to find conditions on when this prediction is true and find any other feasible behaviors.

\subsection*{Linearize Model}
We will first analyze the system (\ref{scaled}) by considering the associated  and Jacobian matrix. With these, we will obtain the equilibrium points and investigate the criteria for their stability. The  $\frac{dx}{d\tau}=0$ and $\frac{dy}{d\tau}=0$ separate positive and negative slope of the variables $x$ and $y$, result in $y=\frac{x^2+(1-A)x}{B(1+x)}$ and $y=x$ as demonstrated in Figure~\ref{fig:6-20-2-pplanealgae}. 

The intersection points of these two ,
\begin{eqnarray*}
(x_0,y_0) &=& (0,0)\\ 
(x_*,y_*) &=& \left(\frac{A+B-1}{1-B},\frac{A+B-1}{1-B}\right), 
\end{eqnarray*}
are the equilibrium points of the system: paired values at which both variables remain constant. Notice the parameter $C=\lambda/\gamma$ does not affect the  or the equilibrium points. This, however, does not eliminate the effect of $\lambda$ or $\gamma$, since $B$ and $A$ depend on them, respectively. To be biologically feasible, equilibria of algae population densities must be nonnegative. To keep the system from digressing into a one trivial equilibrium point, i.e. $(0,0)$, the second equilibrium point must be strictly positive, forcing 
$$
\frac{A+B-1}{1-B}>0.
$$
To maintain this inequality, either the numerator and denominator are both positive or both negative. Due to the assumption of loss in the lower region, $\delta\leq \lambda$, $B=\frac{\delta}{\lambda}\leq 1$. Then, $1-B>0$ where the equality is lost to avoid dividing by zero, and consequentially $A+B-1>0$ as well. Since $A,B,C$ are defined as ratios of nonzero parameters, they are also nonzero. These feasibility criteria of the model are summarized as follows and are summarized in terms of the original parameters in Table \ref{table:feas}.
%Arguing by contradiction, suppose both are negative, $A+B-1>0$ and $1-B<0$. Then $B>1-A$ and $B>1 impossible since it would assume a negative denominator, then $B>1$, which forces $A+B>1$, making it impossible to have a negative numerator. Thus, both are positive, leading us to the feasibility criteria of the model:
\begin{table}[h]
\centering
\begin{tabular}{|lcc|}\hline
Model setup: &$A,B,C>0$ & $\alpha, \beta,\gamma, \delta, \lambda>0$\\ \hline
Circulation: && $\gamma, \delta, \lambda \leq 1$\\ \hline
Loss: &$B<1$ & $\delta<\lambda$\\ \hline
Coexistence: &$A+B<1$ & $\frac{\alpha}{\gamma}+\frac{\delta}{\lambda}< 1 $\\ \hline
\end{tabular}
\caption{Feasibility Conditions on Parameters}
\label{table:feas}
\end{table}

%\begin{eqnarray}
%A,B,C &>& 0 \nonumber \\
%A+B &>& 1 \nonumber \\
%B &<& 1.\label{feas}
%\end{eqnarray}

\subsection*{Analyze Jacobian}
The Jacobian matrix is the matrix of all first partial derivatives of the slope functions, $\frac{dx}{d\tau}=f_1$ and $\frac{dy}{d\tau}=f_2$, from the system (\ref{scaled}), 
$$
J=\left(
\begin{tabular}{cc}
\pwrt{f_1}{x} & \pwrt{f_1}{y} \\
\pwrt{f_2}{x} & \pwrt{f_2}{y} \\
\end{tabular}
\right) 
%=
%\left(
%\begin{tabular}{cc}
%\frac{A}{(1+x)^2-1} & B \\
%C & -C \\
%\end{tabular}
%\right)
.$$
The eigenvalues of the Jacobian matrix determine the behavior of the system. Using the determinant and trace of the Jacobian matrix with their combined discriminant, we can analyze the stability of the equilibrium points by using characteristics of the eigenvalues. This simplifies the task of equilibria characterization due to the complexity of the eigenvalues under three unknown parameters. The determinant, $det()$, and trace, $tr()$, are the product and sum of the eigenvalues of the matrix, respectively. In two dimensions, unstable saddle points are caused by eigenvalues with real parts of opposite sign, making the determinant of a two-dimensional system negative. Stable systems require eigenvalues to have negative real parts, making the trace negative. The discriminant, $\Delta=tr^2(J) - 4 \cdot det(J)$, discriminates between real and complex eigenvalues; nonnegative for the former and negative for the latter \cite{Kreyszig}.

At the origin the Jacobian matrix is
$$
J(0,0)=\left(
\begin{tabular}{cc}
A-1 & B \\
C & -C 
\end{tabular}
\right)
.$$
The determinant of this matrix is
\begin{eqnarray*}
det(J) &=& (A-1)(-C)-BC \\ 
&=& -C(A+B-1),
\end{eqnarray*}
which is negative by the feasibility criteria in Table \ref{table:feas} as $A+B-1>0$ and $A,B,C>0$. Thus, the origin, zero populations in each region, is an unconditional saddle characterized by real eigenvalues of opposite sign.


Evaluated at the second equilibrium point, $(x_*,y_*)=(\frac{A+B-1}{1-B},\frac{A+B-1}{1-B})$,
the Jacobian matrix is
$$
J(x_*,y_*) = \left(
\begin{tabular}{cc}
$\frac{(1-B)^2}{A}-1$ & B \\ 
C & -C 
\end{tabular}
\right)
,$$
with determinant
$$
det(J) = \frac{BC(1-B)}{A}
,$$
which is positive by the feasibility criteria in Table \ref{table:feas} as $1-B>0$ and $A,B,C>0$. So, $(x_*,y_*)$ cannot be a saddle point. The trace of the Jacobian is 
\begin{eqnarray*}
tr(J) &=& \frac{(1-B)^2-A(1+C)}{A}\\
     &<& \frac{(1-B)^2-(1-B)(1+C)}{A} \\
     &=& -\frac{1-B}{A}(B+C),
\end{eqnarray*}
which must be negative as it is bounded above by a negative value. Note that the feasibility criteria in Table \ref{table:feas} imply $-A<-(1-B)$, $1-B>0$, and $A,B,C>0$. Since the determinant is positive and the trace is negative, $(x_*,y_*)$ is a stable sink, either real sink or a spiral sink. To further describe the stability, we look at the discriminant:
\begin{eqnarray*}
\Delta &=& \frac{(1-B)^2-A(1+C)}{A}^2-4\left(\frac{BC(1-B)}{A}\right)\\  
&=& \frac{1}{A^2}[A^2((C+1)^2-4C+4BC)\\
&&+ 2A(C-1)(1-B)^2+(1-B)^4]\\ 
%&=& \frac{A^2((1+C)^2+4C(B-1))+2A(C-1)(1-B)^2+(1-B)^4}{A^2}\\ 
&=& \left(\frac{A(C-1)+(1-B)^2}{A}\right)^2+4BC.
\end{eqnarray*}
Through simplification and completing the square, the discriminant is shown to be positive since $B,C>0$. With a positive discriminant, no cyclic behavior is possible. Therefore, $(x_*,y_*)$ is an unconditionally stable node characterized by two real, negative eigenvalues. Since the other equilibrium points are unstable and lie on the $x,y$ axes and the origin, they create separatrices along the $x,y$ axes which feed any initial condition in the first quadrant to $(x_*,y_*)$. Biologically speaking, any initial nonzero populations of algae in both regions will converge to nonzero equilibrium populations corresponding to the scaled populations $(x_*,y_*)$.

\subsection*{Analyze Parameters}

\begin{figure}[h]
 \centerline{\includegraphics[width=4.0in]{6-20-3-solutionsalgae1b}}
\begin{center}
\vspace*{-6.8 true cm}
\parbox{11 true cm}{\caption{\label{fig:6-20-3-solutionsalgae1b} Numerical solutions to the model (\ref{scaled}) with $A=\frac34, B=\frac{1}{2}, C=\frac14)$, and initial condition $(x_0,y_0)=(1, 0.1)$ or in original parameters $\alpha=\frac35, \beta=\frac12, \gamma=\frac45, \lambda=\frac15, \delta=\frac{1}{10}$, and initial population densities $(u_0,v_0)=(2, 0.8)$.}}
\end{center}

\end{figure}

By direct inspection of (\ref{scaled}), we see that $C$ dictates a limit of convergence in the second equation, rewritten as $x=y+\frac{1}{C} \frac{dy}{d\tau}$. As $C$ increases, $y$ converges to $x$. Considering original variables $u$ and $v$, recall that $y=\beta \frac{\lambda}{\gamma}v = \beta C v$. So as $C$ increases, $v$ diminishes since $y$ converges to $x$. See Table \ref{peffecttable} for a summary of these parameter effects on the population dynamics.

Figures~\ref{fig:6-20-3-solutionsalgae1b} and \ref{fig:6-20-4-solutionsalgae1} are graphs of example solutions taken from the phase plane in Figure~\ref{fig:6-20-2-pplanealgae}. Figure \ref{fig:6-20-3-solutionsalgae1b} demonstrates an influx of algae in the upper region being circulated to a balanced equilibrium between the two regions. Figure~\ref{fig:6-20-4-solutionsalgae1} demonstrates algae self-replenishment where a saturated lower region circulates to replenish a decimated upper population. This example better matches the explanation offered for why algae can be killed off on the surface only to reappear without a consistent frequency.

\begin{figure}%[t]
 \centerline{\includegraphics[width=4.0in]{6-20-4-solutionsalgae1}}
\vspace*{-6.8 true cm}
\begin{center}
\parbox{12 true cm}{\caption{\label{fig:6-20-4-solutionsalgae1} Numerical solutions to the model (\ref{scaled}) with same parameter values as Figure~\ref{fig:6-20-3-solutionsalgae1b} but with an initial condition of $(x_0,y_0)=(0.1, 1)$ representing initial population densities of $(u_0,v_0)=(0.2, 8)$. This solution demonstrates through scaled variables how a decimated upper population of algae ($x$) can be revived through water circulation of the lower storage population ($y$).}}
\end{center}

\end{figure}


% Tables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table}[h!]
\centering
\begin{tabular}{|lp{4.5cm}|}\hline
Parameter Change & \multicolumn{1}{c|}{Effect} \\ \hline
1. $C\rightarrow \infty$ & 1. $y\rightarrow x$ or $\frac{\lambda}{\gamma} v\rightarrow u$. Upper and lower population densities become proportional. Upward circulation is unevenly favored.\\ \hline
2. $A+B\rightarrow 1$ & 2. Population densities increase more steeply. Either resulting from downward circulation greatly exceeding growth rate with a minimal loss, or downward circulation approaches growth rate and loss is great. \\ \hline
3. $B\rightarrow 1$ & 3. Loss is minimized, and as a result the populations are quicker to replenish each other.\\ \hline
4. $C\rightarrow 0$ & 4. Downward circulation is unevenly favored. Loss is minimal. Lower population becomes fixed.\\ \hline
\end{tabular}
\caption{Effect of Parameter Changes}
\label{peffecttable}
\end{table}

\section*{DISCUSSION}
How would our results change if we chose a different form of the model? What is the range of behaviors produced by the original verbal model? To see for yourself, encourage diversity of the models your students produce by only giving them the verbal model and an unrelated example system to analyze, like that given in the student handout. Let them plug it into PPLANE with some randomly chosen parameter values first. Once they have some experience in predicting the type(s) of behaviors, have them go through the analysis to prove the types of behaviors they predicted and find the conditions necessary for it to occur.

From my experience, the main alternatives differ only in the growth of the upper algae: using either a logistic growth term, such as
\begin{eqnarray*}
\frac{du}{dt} &=& \alpha u(1-\beta u)-\gamma u+\delta v,\\
\frac{du}{dt} &=& \gamma u-\lambda v
\end{eqnarray*}
 or a linear growth term, such as 
\begin{eqnarray*}
\frac{du}{dt} &=& \alpha u-\gamma u+\delta v,\\
\frac{du}{dt} &=& \gamma u-\lambda v.
\end{eqnarray*}
The model with logistic growth has similar equilibrium points, stable region, and parameter conditions to the example solution shown below while the linear model either converges to or diverges from the origin, or has an infinite number of equilibrium points similar to the model shown in Figure~\ref{fig:6-20-1-PPLANEdisplayex2}.



} 

\begin{thebibliography}{98}

\bibitem{Tony} Burton, T. 1997. Can Mexico's Largest Lake Be Saved? {\it Ecodecision}. 23:68-71.

\bibitem{CEC} Najera, R.G., Society of Friends of Lake Chapala, {\it et al\/}. 2003. Citizen Submissions on Enforcement Matters: Lake Chapala II. {\it Commission for Environmental Cooperation}.

\bibitem{CEC:response} Garver, G. 2003. Response From Mexican Government: Lake Chapala II. {\it Commission for Environmental Cooperation}.

\bibitem{Kreyszig} Kreyszig, E. 1999. {\it Advanced Engineering Mathematics, 8th ed}. Hoboken, New Jersey: John Wiley and Sons, Inc.

\bibitem{Lopez} Lopez, L. 2000. El Cultivo del Charal en Chapala [the Culture of the Charal in Chapala]. {\it Gaceta Universitaria}. Publication of the University of Guadalajara.

\bibitem{Harwood2008} Manoranjan,V. S., M. A. O. Gomez, and R. C. Harwood. 2008. Modeling Algae Self-Replenishment. {\it Journal of Interdisciplinary Mathematics}, 11: 681-694.

\bibitem{Polking} Polking, J. and D. Arnold. 1999 {\it Ordinary Differential Equations Using Matlab, 2nd ed}. Upper Saddle River, NJ: Prentice Hall. 

\bibitem{Polking1} Polking, J. and J. Castellanos. 2015. {\it PPLANE}. Software available at  \href{https://www.cs.unm.edu/\~joel/dfield/}{https://www.cs.unm.edu/\~joel/dfield/}. 

\bibitem{Encyclopedia} Rittner, D. and T. McCabe. 2004. {\it Encyclopedia of Biology}. New York: Facts on File, Inc.

\end{thebibliography}




\end{document}



