% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
% 
% 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.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\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.
%
% 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 only 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.
%
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{url}
%
\renewcommand{\baselinestretch}{1.3}
%
% 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 and delete the footnote. 
% See above for instructions on toggling from Teacher Version to Student Version.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION %\footnote{DELETE THIS FOOTNOTE BEFORE SUBMISSION.  In the \TeX \, source code there is a toggle/option to produce a Teacher Version or a Student Version by ``carving'' out with  $\backslash$teacher\{\} that material which should only appear in the Teacher version. See the source code.} 
\else STUDENT VERSION \fi\\ Plants Versus Herbivores}




\markboth{Plants Versus Herbivores}{Plants Versus Herbivores}

\author{Mary Vanderschoot\\ 
Wheaton College\\
Wheaton IL 60187 USA\\
mary.vanderschoot@wheaton.edu
}

\makeStitlePDFLaTex


\begin{abstract}

In a recent study \cite{van} of plants and herbivores on an island in the North Sea, ecologists made a surprising observation:  Instead of more vegetation resulting in more grazers, more vegetation resulted in fewer grazers.   Consequently, the ecologists hypothesized that, as the vegetation grew more dense, it became increasingly difficult for herbivores to reach and/or to digest the vegetation.  Eventually, the vegetation might grow so dense as to be unsuitable for the herbivores.  In this activity, students will apply a variety of techniques for analyzing nonlinear systems (e.g., nullclines, linearization, and technologies for drawing phase portraits) to study plant-herbivore models.  In the final section, students will determine the long-term predictions of two mathematical models used in the North Sea study.  
\end{abstract}


\section*{SCENARIO DESCRIPTION}

In a recent study \cite{van} of a plant-herbivore ecosystem on an island in the North Sea, ecologists made a surprising observation:  Instead of more vegetation resulting in more grazers, more vegetation resulted in fewer grazers.   The ecologists hypothesized that, as the vegetation grew more dense, it became increasingly difficult for the herbivores to reach and/or to digest the vegetation.  Eventually, the vegetation might grow so dense that it became unsuitable for the herbivores altogether.   In this modeling scenario, you will apply techniques for studying nonlinear systems to analyze plant-herbivore models, including the two models used in the North Sea study.


\section*{Submission Version}

\section{Background}

In the 1970s \cite{Feng} ecologists began to apply the predator-prey model used in differential equations to model plant (`prey') and herbivore (`predator') interactions.   Let $P$ represent the plant biomass (usually measured in $kg/m^2$) and let $H$ represent the herbivore density (usually measured in number of herbivores or number of droppings per unit of area).   We first consider a simpler case of the standard plant-herbivore model, in which the herbivores are granivores.  Granivores, such as kangaroo rats, feed primarily on dropped leaves or dropped seeds;  they do not subtract biomass from the plant itself.   The plant-granivore model is based on the following assumptions:

\begin{itemize}
\item The plant biomass exhibits logistic growth with intrinsic growth rate $r$ and carrying capacity $K$.  
\item The seed loss is proportional to the plant biomass (i.e. seed loss is $fP$ for some $f>0$.) 
\item The herbivore density grows at a rate $\gamma$ proportional to the seed loss.  
\item The herbivore density decreases at a rate $m$ proportional to its current value. 
\end{itemize}

\noindent These assumptions give the following plant-granivore model:
\begin{eqnarray}
\frac{dP}{dt} & = &rP(1-P/K)-f P\, , \label{eq:1a} \\
\frac{dH}{dt} & = &\gamma fP - mH. \label{eq:1b}
\end{eqnarray} 

\begin{enumerate}

\item Understanding the model:
\begin{enumerate}
\item Explain why it makes sense to use logistic growth for plant biomass.  (What are some of the limitations to plant growth?)
\item In using $fP$ for seed loss, what assumptions are being made about the plant?
\item In your own words, what does $\gamma$ represent for the herbivores?  What does $m$ represent?
\end{enumerate}

\item Consider the differential equation (\ref{eq:1a}) for the plant biomass:
\begin{enumerate}
\item  Explain why (\ref{eq:1a}) does not depend on $H$.  What is the biological meaning?
\item  Find the equilibrium points and sketch the phase line.
\item  Using the phase line, what will be the plant density in the long-run?
\end{enumerate}

\item Suppose the plant-granivore model, (\ref{eq:1a}) and (\ref{eq:1b}), is applied to a population of granivorous kangaroo rats in a desert where  $r=0.25, K=10, f=0.05, \gamma=2,$ and $m=0.5$.
\begin{enumerate}
\item Sketch the nullclines and determine the equilibrium points of the system.
\item Use the Jacobian to classify the equilibrium points. 
\item Use all of this information to sketch the phase portrait.  Based on the phase portrait, what does this model predict about the plant and herbivore densities in the long-run?
\end{enumerate}

\end{enumerate}

In most plant-herbivore ecosystems, the herbivores consume more than just dropped leaves and/or dropped seeds.  The rate at which an herbivore consumes plant biomass has been found to depend primarily on the amount of available vegetation.  (Intuitively, we expect more grazing when there is more grass growth.)  Thus, a simplifying assumption in many plant-herbivore models is that the herbivore consumption rate can be represented by a function of $P$.  

Let $c(P)$ be the herbivore per-capita consumption rate.  The function $c(P)$ is called the \emph{functional response}.   The total herbivore consumption rate is found by multiplying $c(P)$  by $H$.   The plant-herbivore model is obtained by replacing the seed loss in the plant-granivore model with the total herbivore consumption rate: 
\begin{eqnarray}
\frac{dP}{dt} & = &rP(1-P/K) -c(P) H\, , \label{eq2:a} \\
\frac{dH}{dt} & = &n(P) H\, , \label{eq2:b}
\label{Standard}
\end{eqnarray} 
 
\noindent  where $n(P)=(\gamma c(P) - m)$.  The function $n(P)$ is called the \emph{numerical response}, which is the herbivore per-capita growth rate.  

The remaining exercises in this activity require use of a technology for plotting phase portraits.  The phase portraits in this scenario were created using pplane, a free MATLAB application written and developed by John C. Polking \cite{Polking}.  The most recent version of pplane (pplane8) can be found on the MathWorks File Exchange.  (See \cite{harvey}.)   After downloading pplane8.m, place it in a special folder, open the .m file, and then press run.   A setup screen will appear for entering the  system of differential equations and the parameters.  In addition to generating the vector field, pplane8 can be used to sketch trajectories and nullclines as well as to find and classify equilibrium points as saddles, spiral sinks, etc.  (See  \cite{PplaneDirections} for additional directions on downloading and using pplane8.)
 
\begin{enumerate}
\item[4.] Suppose that (\ref{eq2:a}) and (\ref{eq2:b}), with $c(P)= \frac{aP}{b+P}$, is used to model a population of snowshoe hares in a boreal forrest where $K=18$, $r=1$, $a=1$, $\gamma = 0.4$, $b=10$, and $m=0.1$.

\begin{enumerate} 

\item Sketch the graph of $c(P)$.  What is the maximum consumption rate (also called the saturation constant)? What is the biological meaning of  $b$? 

\item Sketch the graph of $n(P)$.  Explain why $n(P)$ is negative for small values of $P$.  What is the biological meaning of this negativity?

\item Note that $(0,0)$ and $(18,0)$ are equilibrium points.  Use technology to find the third equilibrium point.

\item Use technology to sketch the phase portrait.  Observe that solution curves approach a periodic orbit.  Is the direction of this orbit clockwise or counterclockwise?   Based on the phase portrait, what does this model predict about the plant and herbivore densities in the long-run?
\end{enumerate}
\end{enumerate}

Note:  The functional response in Exercise 4  is said to be \emph{saturating} because it satisfies the assumption that the foraging ability of the herbivore increases as more plant biomass becomes available, but it levels off as the herbivore reaches satiation.   

\section{More Grass, Fewer Grazers} 

The herbivores in \cite{van} consisted of geese, hares, and rabbits, which grazed on a standing crop of vegetation in an island salt marsh.  From field observations (using herbivore droppings to measure herbivore density), the ecologists found that the maximum foraging pressure was not at a maximum level of the crop, but at an intermediate level of the crop.  They hypothesized that as the vegetation became taller and/or more dense, it became more difficult for smaller herbivores to digest and/or to reach.    In this section you will analyze two models from this study.   The parameters for the study are $r=1$, $a=1$, $\gamma = 0.4$, $b=10$, and $m=0.1$.

The following system of differential equations was used by van der Koppel \emph{et al.} to model reduced digestion efficiency:   
\begin{eqnarray}
\frac{dP}{dt} & = &rP(1-P/K) -  \frac{aP}{b+P} H\, , \label{eq3:a} \\
\frac{dH}{dt} & = &\left(\gamma  \frac{aP}{b+P} e^{-.065 P} - m\right) H. \label{eq3:b}
\label{Standard}
\end{eqnarray} 

\noindent Note that this is the standard plant-herbivore model given by (\ref{eq2:a}) and (\ref{eq2:b}), with $c(P)= \frac{aP}{b+P}$ and 
$n(P)  = (\gamma  \frac{aP}{b+P}  e^{-.065 P} - m)$.    


\begin{enumerate}
\item[5.] Sketch the graph of $n(P)$.  Observe that $n(P)$ is negative for small and large values of $P$.  Explain what this means biologically.  How does this differ from the numerical response in Exercise 4?  (The parameter values are the same as those in Exercise 4, only the numerical response is different.) 
\end{enumerate}

As in Exercise 4, we can use technology to generate phase portraits for various values of $K$.    Figure \ref{Kis18} shows the phase portrait (including nullclines) for  $K=18$.  Observe that the $P$ and $H$ nullclines intersect at four points (the equilibrium points).   Here $(0,0)$ is a source, $(18,0)$ is a sink,  $(5.6, 10.7)$ is a spiral sink, and $(12, 7.3)$ is a saddle.  From Figure \ref{Kis18} we observe that  solution curves with initial conditions on one side of the saddle's stable separatrix approach the spiral sink, while those on the other side approach $(18,0)$.  
\begin{figure}[htbp]
\hspace*{-.4in}
  \begin{minipage}[b]{0.6\linewidth}
  \centering
    \includegraphics[width=7 cm]{K=18.png}
 
   \caption{Phase portrait for K=18.}
    \label{Kis18}
  \end{minipage}
  \hspace{0.5cm}
  \begin{minipage}[b]{0.4\linewidth}
   \centering
    \includegraphics[width=\textwidth]{phaserep.png}
    \caption{Phase representation for K=18.}
    \label{phaserep}
  \end{minipage}
\end{figure}


The sketch shown in Figure \ref{phaserep} (called a phase representation) provides a summary of all the possible long-term outcomes.  This sketch includes the nullclines (the $P$ nullclines are dashed and the $H$ nullclines are solid) and the equilibrium points.  Equilibrium points are indicated by filled circles for sinks or open circles for sources or saddles.  If the initial plant and herbivore populations, $(P(0),H(0))$, is in the interior of the shaded region  then $(P(t), H(t))$ approaches the sink $(5.6, 10.7)$ as $t \to \infty$.   If $(P(0),H(0))$ is not in the shaded region then the solution curve $(P(t), H(t))$ will approach $(18,0)$, which means that the vegetation becomes too dense for the herbivores to graze.  Ecologists call the shaded region the \emph{domain of attraction} of the plant-herbivore equilibrium and the unshaded region the domain of attraction of the plant equilibrium.  


\begin{enumerate}
\item[6.] Use technology to sketch phase representations for (\ref{eq3:a}) and (\ref{eq3:b}) (as in Figure \ref{phaserep}) for $K = 5, 10, 14, 22$ and $24$.  (In some cases you will find that solutions in the shaded region approach a periodic orbit in the region.)  At what value of $K$ are the herbivores no longer able to keep the vegetation in check?  That is, at this level of $K$, the vegetation will become so dense that it is no longer suitable for grazing, regardless of the initial number of herbivores.
\end{enumerate}

The second model in \cite{van} was also designed to incorporate the assumption of reduced consumption efficiency:  
\begin{eqnarray}
\frac{dP}{dt} & = &rP(1-P/K) -  \frac{aP}{b+P}  e^{-0.065 P} H\, , \label{eq5:a} \\
\frac{dH}{dt} & = &\left(\gamma \frac{aP}{b+P} e^{-.065 P} - m\right) H. \label{eq5:b}
\label{Standard}
\end{eqnarray} 

\begin{enumerate}
\item[7.]  Sketch the functional response for this model and estimate the value of $P$ for which $c(P)$ is maximum.  Is this functional response saturating?  Explain what this means biologically in terms of herbivore consumption efficiency.  


\item[8.] Use technology to sketch phase representations as in Figure \ref{phaserep} for $K = 5, 10, 18,$ and $22$.   How do the results compare with those in Exercise 6 (the reduced digestion efficiency model)?
\end{enumerate}

\begin{thebibliography}{98}


\bibitem{Feng} Feng Z., DeAngelis, D. 2018. {\it Mathematical Models of Plant-Herbivore Interactions}. Boca Raton FL USA: CRC Press.

 \bibitem{harvey} Harvey, Hugh.  2021. Pplane. \url{https://www.mathworks.com/matlabcentral/fileexchange/61636-pplane}.  MATLAB Central File Exchange. Retrieved 8 March 2021. 
 
\bibitem{PplaneDirections} Professor Vladimir Dobrushkin's Pplane Directions.  \url{http://www.cfm.brown.edu/people/dobrush/am33/Matlab/ch1/pplane.html}. Accessed 8 March 2021.

\bibitem{Polking}  Polking,  J.  2004. {\it Ordinary Differential Equations Using MATLAB}. Upper Saddle River NJ USA: Pearson/Prentice Hall.

\bibitem{van} van der Koppel, J., Huisman, J., van der Wal, R., and H.~Olff. 1996. Patterns of herbivory along a productivity gradient: an empirical and theoretical investigation. {\it Ecology} 77: 736-745.



\end{thebibliography}




\teacher{
\section*{NOTES FOR TEACHERS}
This scenario is intended to be used as a group project after students have seen some of the techniques for analyzing nonlinear systems.   It assumes students have experience sketching phase portraits using nullclines, linearization, and technology.  The phase portraits in this scenario were created using pplane8.  Professor Vladimir Dobrushkin's online MATLAB tutorial for differential equations \cite{PplaneDirections} is a student-friendly resource with directions for downloading and using pplane8.   

      
\begin{enumerate}

\item Understanding the Model:

\begin{enumerate}
\item Logistic growth is often used to describe plant growth, which is  limited by availability of resources such as sunlight, water and nutrients.
\item More plant biomass results in more seed loss.  
\item The constant $\gamma$ is the conversion rate from plant biomass to herbivore density. The constant $m$ is the herbivore mortality rate.

\end{enumerate}

\item Analyzing the differential equation (\ref{eq:1a}):

\begin{enumerate}

\item The differential equation for $P$ does not depend on $H$ because the granivores only feed on dropped leaves and/or dropped seeds.

\item The equilibrium points are $P=0$ (source) and $P=\frac{K(r-f)}{r}$ (sink).

\item In the long-run, the plant density approaches  $\frac{K(r-f)}{r}$.

\end{enumerate}

\item Kangaroo rats:
\begin{enumerate}
\item The differential equation for the plants is $\frac{dP}{dt}=0.25P(1-P/10)-.05P$.  The equilibrium solutions are $P=0$ and $P=8$.  In the long-run, the plant density approaches $8$.  The $P-$nullclines are $P=0$ and $P=8$, and the $H-$nullcline is  $H=0.2P$.  The nullclines intersect at the two equilibrium points: $(0,0)$ and $(8,1.6)$.  (See Figure \ref{kangaroo}.)
\begin{figure}[h]
\centering
\includegraphics[scale=0.75]{kangaroo.png}
\caption{A sketch of the nullclines.}
\label{kangaroo}
\end{figure}



\item The Jacobian at $(0,0)$ is $\left(\begin{array}{rr}0.2&0\\0.1&-0.5\end{array}\right)$ with eigenvalues $\lambda_1=0.2$ and $\lambda_2=-0.5$.  So $(0,0)$ is a saddle.  The Jacobian at $(8,1.6)$ is $\left(\begin{array}{rr}-0.2&0\\0.1&-0.5\end{array}\right)$ with eigenvalues $\lambda_1=-0.2$ and $\lambda_2=-0.5$.  So $(8,1.6)$ is a sink. 
\item Students should be able to sketch the phase portrait using parts (b) and (c) without technology.  In the long-run this model predicts that the kangaroo rat density will approach 1.6 and the plant density will approach 8.
\end{enumerate}

\item[4.] Snowshoe Hares:
\begin{enumerate}
\item The graph of $c(P)$ increases with $\lim_{P \to \infty} c(P) = a$.  So the value of $a$ is the maximum consumption rate of the herbivore.  Since $c(b)=a/2$, $b$ is called the \emph{half saturation constant}.
\item At small values of plant density the herbivore cannot consume enough nutrients to survive.
\item The third equilibrium point is $(10/3, 880/81)$.
\item Using technology we see (Figure \ref{hares}) that $(10/3, 880/81)$ is a spiral source and that all solutions approach a periodic orbit with counterclockwise direction.  The plant and herbivore populations oscillate.  The herbivore population decreases when plant density is small and this decrease in herbivore density results in an increase in plant density.  An increase in  plant density results in an increase in herbivores which then results in a lower plant density and the cycle repeats.
\end{enumerate}
\begin{figure}
\centering
\includegraphics[scale=1.5]{hares.png}
\caption{Phase portrait for snowshoe hares created using pplane8.}
\label{hares}
\end{figure}

\item[6.] Reduced Digestion Efficiency:  Figure \ref{phaserep2} shows the phase representations for $K=5, 10, 22,$ and $24$.   In the case $K=5$, there are no interior equilibrium solutions and all solutions approach the plant equilibrium $(5,0)$;  the plant density is too low to support the herbivore population.  For $K=10$, there is one interior equilibrium solution, a sink.  Its domain of attraction is the first quadrant, so the model predicts that the plant density reaches a level that can support the herbivores.  The case of $K=14$ is similar to that of $K=18$, but with a larger domain of attraction.  For $K=22$, there are two interior equilibrium solutions, one saddle and one spiral source.  Solutions in the shaded area approach a periodic orbit.  At $K=24$ the plant vegetation is so dense that it is unsuitable for herbivore grazing.  
\begin{figure}
\centering
 \includegraphics[scale=1.4]{phaserep2.png}
 \caption{Phase Representations for Reduced Digestion Efficiency \cite{van}.}
 \label{phaserep2}
 \end{figure}
 
 \item[8.] Reduced Digestion and Reduced Consumption Efficiency:  The predictions for $K=5$ and $K=10$ are similar to those in Exercise 4.  The case of $K=14$ is similar to the behavior of $K=22$ in Exercise 6.  At  $K=18$ the plant vegetation is so dense that it is unsuitable for herbivore grazing.


\end{enumerate}

}

\end{document}



