\SweaveOpts{eval=TRUE, prefix=FALSE, width=4, height=4, eps=FALSE, include=FALSE}
\chapter{Density-dependent Growth}
\label{ch:DDG} % Always give a unique label
% use \chaptermark{}
% to alter or adjust the chapter heading in the running head

\label{sec:1}
% Always give a unique label
% and use \ref{<label>} for cross-references
% and \cite{<label>} for bibliographic references
% use \sectionmark{}
% to alter or adjust the section heading in the running head

Let's go back to our Song Sparrow (\index{Melospiza@\emph{Melospiza melodia}}\emph{Melospiza melodia}) data from Chapter 1 on density-independent growth --- now we look at \emph{all} the data.

<<echo=false, results=hide>>=
setwd("~/Documents/Projects/Book/draft/Chap03")
rm(list=ls())
badfonts = FALSE
.detach.stuff = function() {
 s1 = grep("package|Autoloads",search())
 nattach = length(search())
 xattach = search()[c(-1,-s1)]
 for (i in xattach)
   eval(substitute(detach(i),list(i=i)))
}
.detach.stuff()
## ps.options(family="NimbusSan")
ps.options(colormodel="cmyk")
options(width=75, SweaveHooks=list(fig=function() par(mar=c(5,4,.5,1.5))) ) 

## palette(gray((0:8)/8))
library(primer)
trellis.par.set(canonical.theme(color=FALSE))

<<echo=false, results=hide>>=
data(sparrows)
attach(sparrows)
obs.Rs <- Count[-1]/Count[-length(Count)]
<<MelospizaDDG1, fig=TRUE,echo=false, results=hide>>=
plot(Year,Count, type="b")
<<MelospizaDDG2, fig=TRUE,echo=false, results=hide>>=
plot(obs.Rs ~ Count[-length(Count)], 
     ylab=expression("N"[t+1]/"N"[t]),
     xlab="Sparrow Count", type="n") 
text(Count[-length(Count)], obs.Rs,
       substr(as.character(Year), 3,4 ) )
abline(mod <- lm(obs.Rs ~ Count[-length(Count)] ), 
       lty=2 )
#newd <- 20:120
#cts <- Count[-length(Count)]
#mod <- loess(obs.Rs ~ cts )
#lines(newd, predict( mod, newdata=data.frame(cts=newd) ))
@

\begin{figure}[ht]
  \centering
  \subfloat[Song Sparrow Counts]{\includegraphics[width=.47\linewidth]{MelospizaDDG1.pdf}   \label{fig:MeloDDG1}}
  \subfloat[$N_{t+1}/N_t$ $vs$. Counts]{\includegraphics[width=.47\linewidth]{MelospizaDDG2.pdf}   \label{fig:MeloDDG2}}
  \caption{Song Sparrow \emph{Melospiza melodia} counts from 1966--2003 and the relation between observed counts and annual growth rate determined from those counts, fit with ordinary least squares regression. See Chapter 1 for data source.}

\end{figure}

In Chapter 1, we modeled Song Sparrow growth without reference to the abundance, or counts, of sparrows, using a small subset of the data. When we look at all the available data (Fig. \ref{fig:MeloDDG1}), it seems clear that the annual growth rate ($R_t=N_{t+1}/N_t$) depends on the density of the bird population (Fig. \ref{fig:MeloDDG2}). The larger the population, the lower $R$ becomes, until above about 70 counted birds, $R < 1$, indicating population decline. What might limit the population growth of these sparrows? Space available for territories in their successional-scrub type habitat? Availability of seeds, fruits and invertebrates? We don't necessarily know precisely what limits it, but if it is something related to their own abundance, perhaps we can study its \emph{density-dependence}.

Density-dependent population growth is the case where the per capita population growth rate depends statistically on the density of the population. Ecologists consider that negative \index{density-dependence}density-dependence is typically a characteristic of a population undergoing \emph{intraspecific competition}. That is, individuals of the same species are competing for shared resources. So, \ldots how would we represent this algebraically?\footnote{root: Arabic. \emph{Al-jabr}}


\section{Discrete Density-dependent Growth}
\subsection{Motivation}
To begin, we recall our model of geometric growth for a single time
step,

\begin{equation}
N_{t+1}=\lambda N_{t}.
\end{equation}
We can pull this apart into two sections by recalling that we can decompose $\lambda$ into two parts, $\lambda=1+r_{d}$, where
$r_{d}$ is the \emph{discrete growth factor}. This allows
us to state,
\begin{equation}
N_{t+1}=\lambda N_{t}=N_{t}\left(1+r_{d}\right)=N_{t}+r_{d}N_{t}.
\end{equation}
Here we see the $N_{t+1}$ is equal to the previous
year's population, $N_t$, plus a proportional change, $r_dN_t$. To add density dependence (e.g., Fig. \ref{fig:MeloDDG2}), we can build
\emph{density-dependence} into that proportional change. 

Density dependence means that the population will grow or shrink at a rate that depends on its size. We can imagine that each individual in the population exerts some tiny negative effect on $r_dN_t$, so that the realized per capita growth increment shrinks as $N_t$ grows.

Let us specify that the negative effect of each individual in the population is the same, regardless of how many individuals there are. That is, we could represent this negative effect with a constant, perhaps $\alpha$,  where $\alpha N_t$ is the total negative effect of all individuals in the population. On average, each individual exerts the same negative impact, with a magnitude $\alpha$, whether there is one individual or 1000 individuals. A particularly convenient way of implementing this is to keep using $r_d$ but to multiply it by a scaling factor. This scaling factor should equal 1 when the population size is zero so that $r_d \times \left(\mathrm{scaling~factor}\right) = r_d$.
 The scaling factor should be zero when the population is so big that \index{per capita growth increment}per capita growth increment\index{rd@$r_d$|see{per capita growth increment}} is zero, $r_d \times \left(\mathrm{scaling~factor}\right) = 0$. One such scaling factor looks like this,
\begin{equation}
  \label{eq:1}
  \mathrm{Per~Capita~Increment}=r_d\left(1-\alpha N_t\right)
\end{equation}
where $\alpha$ is the per capita negative effect of an individual on the per capita growth increment. As $N_t$ shrinks toward zero, this expression grows toward $r_d$; as $N_t$ grows, this expression shrinks toward zero. 

At precisely what value of $N$ will per capita growth shrink to zero? We can solve this by noting that $r_d$ is a constant and won't change; all that matters is inside the parentheses. We set the per capita increment equal to zero, and solve for $N_t$.
\begin{align}
  0 &= r_d\left(1 - \alpha N_t\right) \label{eq:2add}\\
  0 &= r_d - r_d\alpha N_t\\
  N_t &= \frac{1}{\alpha}.\label{eq:2add3}
\end{align}
When the population reaches the $N=1/\alpha$, the growth increment will be zero, and the population will stop growing. 

In a sense, the algebraic rearrangement (eq. \ref{eq:2add}) is at the core of theoretical ecology. We began with a set of assumptions (constant $r_d$ and $\alpha$), assumed a relation between them ($r_d\left[1-\alpha N_t\right]$), and examined one of the consequences: the population will stop changing when it reaches a density of $1/\alpha$, \emph{ceretus paribus}.\footnote{all else being equal} Welcome to theoretical ecology.

Now instead of density-independent per capita growth, $r_d$, we have density-dependent per capita growth $r_d\left(1 - \alpha N_t\right)$, and our population growth equation becomes,

\begin{equation}
  \label{eq:3dd}
  N_{t+1} = N_t + r_d\left(1-\alpha N_t\right)N_t
\end{equation}
This describes \emph{discrete \index{logistic growth!discrete}logistic growth}. A common alternative representation uses $1/\alpha$ symbolized as $K$,
\begin{equation}
  \label{eq:4dd}
  N_{t+1} = N_t + r_d\left(1-\frac{N_t}{K}\right)N_t
\end{equation}
where $K$ represents the carrying capacity. The \index{carrying capacity}carrying capacity is the population size at which the per capita growth increment has fallen to zero. The context dictates whether we prefer to represent the scaling factor as the per capita effect, $\alpha$, or the population carrying capacity, $K$.

\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{Writing a Function For Discrete Logistic Growth}
An \R~\emph{function} will simplify our explorations. It will return a vector of $N$, given $\alpha$, $r_{d}$, $N_0$, and $t$. The function \emph{arguments} can have defaults (e.g., \texttt{t=10}).
<<>>=
dlogistic <- function(alpha=0.01, rd=1, N0=2, t=15) {
  ## make a holder for initial N and all subsequent N.
  N <- c(N0, numeric(t))
  ## Use a for loop to create the population dynamic. We 
  ## CANNOT use lapply or sapply, because each N depends on
  ## the previous N.
  for(i in 1:t) N[i+1] <- {N[i] + rd * N[i] * (1 - alpha * N[i])}
  return(N) }
@ 
The function first makes a new vector containing $N_0$, uses a for-loop to implement eq. \ref{eq:3dd} for each time step, and then returns the vector $N$.
}
\end{boxedminipage} \medskip

With discrete logistic growth, if we start with a small population ($N<<K$), we will see the population rise and gradually approach $K$ or $1/\alpha$ (Fig. \ref{fig:dlg}). We refer to $K$ as an \emph{\index{attractor}attractor} because $N$ moves in a deterministic fashion toward $K$. We explore the meanings of \emph{attractor} and related terms throughout the book.

  \begin{figure}[ht]
    \centering
    \includegraphics[width=.5\linewidth]{DLG1.pdf}
    \caption{Discrete logistic growth with $r_d=1$, $\alpha = 0.01$. }
    \label{fig:dlg}
  \end{figure}

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Graphing Population Size}
We can use the function created above, \texttt{dlogistic}, with default settings, to generate a population projection.
<<>>=
Nts <- dlogistic()
@ 
Now we plot the projection, and put in a dotted line for $1/\alpha$ or $K$.
<<DLG1, fig=TRUE>>=
t <- 15
a <- 0.01
plot(0:t, Nts)
abline(h=1/a, lty = 3)
@
}\end{boxedminipage} \medskip

\subsection{Relations between growth rates and density}
We already know a lot about one description of density dependent growth, the discrete logistic model. In particular, we know that with a constant per capita negative effect, $\alpha$,  the population size at which growth falls to zero is $K$. Let us explore further how the \emph{per capita} growth increment, and the \emph{population} growth increment, vary with population size.

Casual examination of Fig. \ref{fig:dlg} suggests that the total population growth increment ($\Delta N_t = N_{t+1}-N_t$) starts out small when both $t$ and $N$ are small, accelerates as $N$ grows, and then over time, slows down and reaches an asymptote of $K$. Is this changing rate of population growth a function of time, or a function of density? Let us first consider the growth increment as a function of $N$.

First consider the relation between the population growth increment and population size (Fig. \ref{pgi}). We see it increase as $N$ grows, and then decrease as $N$ approaches $K$. The pattern is fairly symmetric. That is, it increases and decreases at about the same rates. 

Next consider the per capita growth increment ($\Delta N_t / N_t$; Fig. \ref{ppgi}). There is a direct linear relation between the per capita growth increment and the size of the population --- this is \emph{linear density dependence}. This linear dependence on $N$ comes from our assumption that the per capita negative effect is a constant, $\alpha$.

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{(Per Capita) Population Growth Increment $vs$. $N$ (Fig. \ref{fig:DGI1})}
Using the previous projection, we now capture both the total and the per capita growth increment per unit time, from $t$ to $t+1$. We graph these versus $N_t$, population size at $t$.
<<>>=
total.incr <- Nts[1:t+1] - Nts[1:t] 
per.capita.incr <- total.incr/Nts[1:t]
<<PGI, fig=true, include=false>>=
plot(Nts[1:t], total.incr)
<<PPGI, fig=true, include=false>>=
plot(Nts[1:t], per.capita.incr)
@ 
}\end{boxedminipage} \medskip

\begin{figure}[ht]
  \centering
  \subfloat[Population Growth Increment]{\includegraphics[width=.47\linewidth]{PGI.pdf} \label{pgi}}
  \subfloat[Per capita Growth Increment]{\includegraphics[width=.47\linewidth]{PPGI.pdf} \label{ppgi}}
  \caption{Relations between the total and per capita discrete growth increments and population size.}
  \label{fig:DGI1}
\end{figure}

Let's use a simple analytical approach to understand Figs. \ref{pgi} and \ref{ppgi} a little better. Take a look at eq. \ref{eq:3dd}. First let's rearrange the expression so that we can set the \emph{increment} of change equal to zero. How would you do that? Take a minute and try it, before you look at the answer below. 

Setting the increment of change equal to zero, we rearrange with the population growth increment on the left, and then set it to zero.
\begin{align}
N_{t+1}-N_{t} &= r_{d}N_{t}\left(1-\alpha N_{t}\right) \label{eq:pgi}\\
N_{t+1}-N_{t} &= r_{d}N_{t} - r_d\alpha N_{t}^2 \label{eq:quad}\\
0 &= r_{d}N_{t} - r_d\alpha N_{t}^2 \label{eq:quad0}
\end{align}
What do we notice? One observation we could make is that we have a \index{quadratic equation}quadratic equation,\footnote{$ax^2+bx+c=0$} where the intercept is zero. This tells us that perhaps Fig. \ref{pgi} is symmetric because it is a quadratic expression in terms of $N$. 

What would satisfy this quadratic expression (eq. \ref{eq:quad0}), that is, cause the growth increment to equal zero? Well, if $r_{d}$ or $N_{t}$ equal zero, those would yield potentially interesting solutions. Assuming neither $r_{d}$ nor $N_{t}$ equal zero, we can divide each side by these, and we are left with the solution we found in eq. \ref{eq:2add3}, that the growth increment will be zero when $N_{t}=\frac{1}{\alpha}=K$.\label{sec:per-capita-popul}

Now let us examine the per capita growth increment (Fig. \ref{ppgi}). If we start with the population growth increment  eq. \ref{eq:pgi}, all we need is to divide through by $N_t$ to get
\begin{equation}
\frac{N_{t+1}-N_{t}}{N_t} = r_{d} - r_d\alpha N_{t}. \label{eq:ppgi}\\
\end{equation}
With $r_d$ and $\alpha$ being constants, and $N$ varying, what is this expression? It is the expression for a straight line,\footnote{$y=mx+b$.} just like we observe (Fig. \ref{ppgi}). When $N_t=0$, the per capita increment equals $r_d$, and when $N_t=1/\alpha$, the per capita increment is zero. This is precisely where we started when we began the motivation for discrete logistic growth. 

\subsection{Effect of initial population size on growth dynamics}
What will be the effect of differences in initial population size? 
We could approach such a question in at least two ways \cite{Case2000}. 
For some of us, the simplest way is to play games with numbers, trial and error, plugging in a variety of initial population sizes, and see what happens to the dynamics. If we do this systematically, we might refer to this as a simulation approach. For very complicated models, this may be the \emph{only} approach.  Another approach, which is often used in concert with the simulation approach, is the analytical approach. We used this above, when we set the growth equation equal to zero and solved for $N$. In general, this analytical approach can sometimes give us a definitive qualitative explanation for \emph{why} something happens. This has been used as a justification for using simple models that actually have analytical solutions --- they can provide answers \cite{May2001}.

For an analytical approach, first consider the endpoint solutions to the discrete logistc model eq. \ref{eq:quad}. The population will stop changing when $N_t=K$. Note that it does not appear to matter what the initial population size was. The only thing that matters is $\alpha$. Recall also that the population would not grow if for any reason $N_t=0$ --- the population will be stuck at zero. Based on these analyses, it appears that the only thing that matters is whether the initial population size is zero, or something greater than zero. If the latter, then initial population size appears to have no effect on the eventual population size.

It always pays to check our analytical answer with a brute force numerical approach, so we will use a little simple simulation to see if we are right. In this case, we can vary systematically the initial population size, and see what happens (Fig. \ref{fig:init}). What our approach shows us is that regardless of the initial conditions (except zero), $N$ converges on $K$ --- $K$ is an attractor. We also might notice that sometimes when $N_0>K$, it crashes below $K$ before converging on $K$ in the end (more on that later). Last, because there is a qualitative shift in the behavior of the population when $N=0$, we might want to investigate what happens when $N$ gets very very close to zero. However, in this situation, the analytical solution is so straightforward that it seems convincing that as long as $N>0$, it will grow toward $K$.

\begin{figure}[ht]
  \centering
  \subfloat[Variation in $N_0$]{\includegraphics[width=.48\linewidth]{DLGinit.pdf} \label{fig:init}}
  \subfloat[Variation in $\alpha$]{\includegraphics[width=.48\linewidth]{DLGalpha.pdf}\label{fig:a}}
  \caption[Confirmation of analytical solutions for discrete logistic growth]{\subref{fig:init} Dynamics due to different initial $N$ (zero was also specifically included, and $\alpha=0.01$).  \subref{fig:a} Dynamics due to different $\alpha$.   All $N$, except $N=0$, converge on $K=1/\alpha$, regardless of the particular value of $\alpha$ ($r_{d}=1$).}
\end{figure}


\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Numerical Evaluation of Initial Conditions  (Fig. \ref{fig:init})}
Here we draw randomly 30 $N_0$ from a uniform distribution between zero and $1.2K$. We also include zero specifically. We then use \texttt{sapply} to run \texttt{dlogistic} for each $N_0$, using defaults for the other arguments.
<<>>=
N0s <- c(0, runif(30) * 1.1 * 1/a)
N <- sapply(N0s, function(n) dlogistic(N0=n) )
<<DLGinit, fig=TRUE>>=
matplot(0:t,N, type="l", lty=1, lwd=.75, col=1)
text(t, 1/a, expression(italic("K") == 1/alpha), adj=c(1,0))
@ 
A serious simulation might include a much larger number of $N_0$. 
}\end{boxedminipage}
\medskip

\subsection{Effects of $\alpha$}

Our conclusions thus far have been based on specific values of $\alpha$ and $r_{d}$. Have we been premature? Just to be on the safe side, we should probably vary these also.

What will happen if $\alpha$ varies? This seems easy. First, when $N_t$ is zero, the population growth increment eq. \ref{eq:pgi} is zero, regardless of the magnitude of $\alpha$. However, when $N_t > 0$, $N$ will increase until it reaches $1/ \alpha$ ($K$; Fig. \ref{fig:a}).  The outcome seems pretty clear --- by decreasing the negative effect of individuals on each other (i.e. decrease $\alpha$) then the final $N$ increases, and $\alpha$ determines the final $N$.

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Numerical Evaluation of $\alpha$  (Fig. \ref{fig:a})}
Here we draw 30 random $K$ from a uniform distribution from 50 to 1000, and convert these to $\alpha$. We use \texttt{sapply} to run \texttt{dlogistic} for each $\alpha$.
<<>>=
a.s <- 1/runif(30, min=50, max=1000)
N <- sapply(a.s, function(a) dlogistic(alpha=a, t=15) )
@ 
We next plot all populations, and use some fancy code to add some informative text in the right locations.
<<DLGalpha, fig=TRUE>>=
matplot(0:t, N, type="l", ylim=c(0,1000), lty=1, lwd=.75, col=1)
text(8, 1/min(a.s), bquote(italic(alpha) == .(round(min(a.s),3))), 
     adj=c(1,0.5))
text(10, 1/max(a.s), bquote(italic(alpha) == .(round(max(a.s),3))), 
     adj=c(0,1.2))
@ 
Note that we use the minimum and maximum of \texttt{a.s} to both position the text, and provide the values of the smallest and largest $\alpha$.
}\end{boxedminipage}
\medskip


\subsection{Effects of $r_{d}$} 
What will variation in $r_{d}$ do? \index{logistic growth!effect of $r_d$}Probably nothing unexpected, if our exploration of geometric growth is any guide. Our analytical approach indicates that it should have no effect on $K$ (sec. \ref{sec:per-capita-popul}). Nonetheless, let us be thorough and explore the effects of $r_{d}$ by varying it systematically, and examining the result.

Yikes --- what is going on (Fig. \ref{fig:DDGr1})?. Perhaps it is a good thing we decided to be thorough. These wild dynamics are real --- let's go back and look more carefully at $r_{d}$. 

\begin{figure}[ht]
  \centering
  \includegraphics[width=.67\linewidth]{DLGrd.pdf}
  \caption{The variety of population dynamics resulting from different values of $r_d$ for the discrete logistic growth model ($r_d = 1,1.2, \ldots, 3$, $\alpha=0.01$). See Fig. \ref{fig:DDGr2} for a more informative view.}
  \label{fig:DDGr1}
\end{figure}

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Simple Numerical Evaluation of $r_d$ (Fig. \ref{fig:DDGr1})}
Here we vary $r_d$ by creating a short systematic sequence $r_d = 1.3,1.6, \ldots, 2.8$. We set $t=50$, and use \texttt{dlogistic} to create a trajectory for each of the six $r_d$.
<<>>=
rd.v <- seq(1.3, 2.8, by=.3)
t <- 15
Ns <- data.frame(sapply(rd.v, function(r) dlogistic(rd=r, t=t) ))
<<DLGrd, fig=TRUE>>=
matplot(0:t, Ns, type="l", col=1)
@ 
Note that many populations do not seem to settle down at $K$.
}\end{boxedminipage}
\medskip

If we examine each projection separately, we see a cool pattern is emerging (Fig. \ref{fig:DDGr2}). At the lowest $r_{d}$, the population grows gradually toward its carrying capacity, $K$, and stays there. Once $r_{d}=1.6-1.9$ it  overshoots $K$ just a bit, creating oscillations; these oscillations, however, dampen back down to $K$. When $r_{d}=2.2$, however, the populations seem to bounce back and forth between two values. When $r_d=2.5$, $N$ bounces around, but now it bounces around between four different values.  When $r_{d}=2.8$, however, it seems to bounce back and forth around $K$, but at values that vary every time. This model is just about as simple as a model can be, and includes no random numbers. What is going on?

\begin{figure}[ht]
  \centering
  \includegraphics[width=.9\linewidth]{DDGr2.pdf}
  \caption{A more informative view of the effects of variation in $r_{d}$ on population dynamics.}
  \label{fig:DDGr2}
\end{figure}

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Presentation of Limit Cycles (Fig. \ref{fig:DDGr2})}
First we make a data frame with the six $r_d$ values in the first column, and the respective populations in rows, using \texttt{t()} to transpose \texttt{Ns}. This puts the data in \emph{wide} format, with a different time step in each column. (This might, for instance, be how you record data in a field experiment with repeated measurements through time).
<<>>=
tmp <- data.frame(rd=as.factor(rd.v), t(Ns))
@
Next, we reshape the data to \emph{long} format, were all $N$ are in the second column, and each is associated with a time step and its $r_d$ value (cols. 3 and 4).
<<results=hide>>=
Ns2 <- reshape(tmp, varying=list(2:ncol(tmp)), idvar="rd", v.names="N", 
               direction="long")
str(Ns2)
@ 
\dots (output omitted) We plot each trajectory separately using \texttt{xyplot} in a different graphics package, \texttt{lattice}. Known as a \emph{conditioning} plot, \texttt{xyplot} graphs  $y~vs.~x$ conditional on $g$ (\texttt{y ~ x | g}).
<<DDGr2, fig=TRUE, width=7, height=4>>=
library(lattice)
print( xyplot(N ~ time|rd, data=Ns2, type="l", layout=c(3,2,1), col=1 ) )
@
}\end{boxedminipage} \medskip


What is going on is the emergence of \emph{\index{stable limit cycles}stable limit cycles}, and \emph{\index{chaos}chaos}.\footnote{Not the evil spy agency featured in the 1960's US television show, \emph{Get Smart}.} At low $r_d$, we have simple asymptotic approach to $K$. As $r_d$ increases, we see the population overshoot the carrying capacity and exhibit \emph{\index{damped oscillations}damped oscillations}. When $2 < r_d < 2.449$, the population is attracted to two-point limit cycles. In this case, \emph{these two points are stable attractors}. Regardless where the population starts out, it is attracted to the same two points, for a given $r_d$. As $r_d$ increases further, the number of points increases to a four-point limit cycle (e.g., at $r_d=2.5$), then an eight-point cycle, a 16-point limit cycle, and so on. These points are stable attractors. As $r_d$ increases further , however, stable limit cycles shift into \emph{chaos} ($r_d>2.57$). Chaos is \emph{a non-repeating, deterministic fluctuating trajectory, that is bounded, and sensitive to initial conditions}.

Robert May \cite{May:1974fk} shocked the ecological community when he first demonstrated stable limit cycles and chaos using this model. His groundbreaking work, done on a hand calculator, showed how very complicated, seemingly random dynamics emerge as a result of very simple deterministic rules. Among other things, it made population biologists wonder whether prediction was possible at all. In general, however, chaos seems to require very special circumstances, including very high population growth.

Is there a biological interpretation of these fluctuations? Consider some simple environment, in which small vegetation-eating animals with high reproductive rates eat almost all the vegetation in one year. The following year, the vegetation will not have recovered, but the animal population will still be very high. Thus the high growth rate causes a disconnect between the actual population size, and the negative effects of those individuals comprising the population. \emph{The negative effects of the actions of individuals (e.g., resource consumption) are felt by the offspring of those individuals, rather than the individuals themselves}. We won't belabor the point here, but it is certainly possible to extend this delayed density dependence to a wide variety of populations. The discrete logistic model has a built in delay, or \emph{\index{time lag}time lag}, of one time step, because the growth increment makes a single leap of one time step. This delay is missing from the analogous continuous time model because the growth increment covers an infinity small time step, thanks to the miracles of calculus.\footnote{A time lag can be explicitly built in to a continuous time model, with only a small effort.}

\subsubsection{Bifurcations}
Up until now, we have examined $N$ as a function of time. We have
graphed it for different $\alpha$ and $N_0$, but time was always on
the X-axis. Now we are going to examine $N$ as a function of $r_d$, so
$r_d$ is on the X-axis. Specifically, we will plot the stable limits
or attractors $vs$. $r_d$ (Fig. \ref{fig:bifurcation}). 
 What does it mean? For $r_{d} < 2$, there is only a single $N$. This is what we mean by a stable point equilibrium, or point attractor --- as long as $r_d$ is small, $N$ always converges to a particular point.\footnote{It need not be the same $N$ for each $r_d$, although in this case it is.} When $2<r_{d}<2.45$, then all of a sudden there are two different $N$; that is, there is a two-point stable limit cycle. Note that when $r_d\approx 2$ these oscilliations between the two point attractors around $K$ are small, but as we increase $r_d$, those two points are farther apart. The point at which the limit cycle emerges, at $r_d=2$, is called a \emph{\index{bifurcation}bifurcation}; it is a splitting of the single attractor into two attractors. At $r_{d}\approx 2.45$, there is another bifurcation, and each the two stable attractors split into two, resulting in a total of  four unique $N$. At $r_{d}\approx 2.53$, there are eight $N$. All of these points are \emph{periodic} attractors\index{attractor!periodic} because $N$ is drawn to these particular points at regular intervals. As $r_d$ increases the number of attractors will continue to double, growing geometrically. Eventually, we reach a point when there becomes an infinite number of unique points, \emph{that are determined by $r_{d}$}.\footnote{They are also determined by the initial $N$, but we will get to that later.} This completely deterministic, non-repeating pattern in $N$ is a property of \emph{chaos}. Chaos is not a random phenomenon; rather it is the result of deterministic mechanisms generating non-repeating patterns.

\begin{figure}[ht]
  \centering
  \includegraphics[width=.67\linewidth]{Bifurcation}
  \caption{Illustration of the long term dynamics of discrete logistic population growth. When a small change in a continuous parameter results in a change in the number of attractors (e.g. a single point equilibrium to a stable 2-point limit cycle), we call this a bifurcation.}
  \label{fig:bifurcation}
\end{figure}

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Bifurcation Plot: Attractors as a Function of $r_d$ (Fig. \ref{fig:bifurcation})}
Here we perform more comprehensive simulations, and plot the point and periodic attractors \emph{vs}. $r_d$. First we pick some constraints for the simulation: the number of different $r_d$, the sequence of $r_d$ values, and the number of time steps.
<<eval=true>>=
num.rd <- 201
rd.s <- seq(1,3, length=num.rd)
t <- 400
<<echo=false, results=hide, eval=FALSE>>=
num.rd <- 1001
rd.s <- seq(1,3, length=num.rd)
t <- 1000
@ 
Next we use \texttt{sapply} for the simulations (If \texttt{num.rd} and \texttt{t} are large, this could take a while).
<<>>=
tmp <- sapply(rd.s, function(r) dlogistic(rd=r, N0=99, t=t) ) 
@ 
Next we convert the output to a data frame and stack up the $N$ in one column. We also rename each of the stacked columns, and add new columns for the respective $r_d$ and time steps.
<<>>=
tmp.s <- stack( as.data.frame(tmp) )
names(tmp.s) <- c("N", "Old.Column.ID")
tmp.s$rd <- rep(rd.s, each=t+1)
tmp.s$time <- rep(0:t, num.rd)
@  
We save just the later dynamics in order to focus on the $N$ after they have converged on the periodic attractors. Here we select the last 50\% of the time steps. (Your figure will look a little different than Fig. \ref{fig:bifurcation} because I used more $r_d$ and time steps.)
<<Bifurcation,fig=TRUE>>=
N.bif <- subset(tmp.s, time > 0.5 * t )
plot(N ~ rd, data = N.bif, pch=".", xlab=quote("r"["d"]))
@ 
}\end{boxedminipage} \medskip

There has been a great deal of effort expended trying to determine whether a particular model or real population exhibits true \index{chaos}chaos. In any practical sense, it may be somewhat unimportant whether a population exhibits true chaos, or merely a higher order periodic attractor \cite{Ellner:2005fk}. The key point here is that very simple models, and therefore potentially simple mechanisms, can generate very complex dynamics.

@
\subsubsection{Sensitivity to initial conditions}
Another very important characteristic feature of chaotic populations is that they are very sensitive to initial conditions. Thus emerges the idea that whether a \index{butterflies}butterfly in Sierra Leone flaps its wings twice or thrice may determine whether a hurricane hits the southeastern United States in New Orleans, Louisiana, or in Galveston, Texas.\footnote{Clearly, this suggests that the solution to increased storm severity due to global warming is to kill all butterflies.}

If we generate simulations where we vary initial population size by a single individual, we find that this can have an enormous impact on the similarity of two populations' dynamics, and on our ability to predict future population sizes (Fig. \ref{fig:DDGChaosInitN}). Note how the populations start with similar trajectories, but soon diverge so that they experience different sequences of minima and maxima (Fig. \ref{fig:DDGChaosInitN}). This is part of what was so upsetting to ecologists about May's 1974 paper --- perhaps even the simplest deterministic model could create dynamics so complex that we could not distinguish them from random \cite{May:1974fk}. Over time, however, we came to learn that (i) we could distinguish random dynamics from some chaos-like dynamics, and (ii) the hunt for chaos could be very exciting, if most frequently disappointing \cite{Becks2005, Constantino1995,Kendall:1998ys}.

\begin{figure}[ht]
  \centering
  \includegraphics[width=.67\linewidth]{DDGChaosInitN.pdf}
  \caption{Effects of differences in initial population size on the short term and long term dynamics, and their correspondence, of three populations.}
  \label{fig:DDGChaosInitN}
\end{figure}

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
  \paragraph{Sensitivity to Intitial Conditions}
We start with three populations, all very close in initial abundance. We then propogate with a $r_d$ to generate chaos for 100 time steps.
<<>>=
N.init <- c(97,98,99)
t <- 30
Ns <- sapply(N.init, function(n0) dlogistic(rd=2.7, N0=n0, t=t) )
@ 
Now we would like to graph them over the first 12 times, and look at the correlations between $N_1$ and the other two populations.
<<DDGChaosInitN, fig=TRUE>>=
matplot(0:t,Ns, type="l", col=1)
@
}\end{boxedminipage} \medskip

\subsubsection{Boundedness, and other descriptors}
One last issue that we should note is the extent to which our populations are \emph{\index{chaos!boundedness}bounded}. A population may have complex dynamics, but we may be able to characterize a given population by its upper and lower bounds. In spite of the differences created by the initial conditions, the upper and lower bounds of our chaotic populations were very similar (Fig. \ref{fig:DDGChaosInitN}). Note also as $r_d$ increases (Fig. \ref{fig:bifurcation}) the oscillations increase very systematically. 

In general, we can describe many characteristics of populations, even if we cannot predict exactly what the population size will be five years hence. For instance, we can describe the shape of density dependence (linear, nonlinear), and characteristics of $N$, such as the average, the variance, the periodicity, the upper and lower bounds, and the color (i.e. degree of temporal auto-correlation). Indeed, these characteristics may vary among types of organisms, body sizes, or environments.


\section{Continuous Density Dependent Growth}

The classic version of continuous density dependent growth \cite{Kingsland:1985kx} is the continuous logistic growth equation, the continuous version of eq. \ref{eq:4dd} above, \index{K@\emph{K}}
\begin{equation}
\label{clogisticK}
\frac{\D{N}}{\D{t}} = rN\left( \frac{K-N}{K} \right)
\end{equation}
which, as we saw above, is no different than
\begin{equation}
\label{clogistica}
\frac{\D{N}}{\D{t}} = rN\left( 1-\alpha N \right).
\end{equation}
Take a moment to see the parallels between these and eq. \ref{eq:3dd}.

\subsection{Generalizing and resimplifying the logistic model}
\index{logistic growth!generalizing}
By now you might be wondering why we are studying such simplistic models that make such unrealistic assumptions. Let's use the next couple of pages to figure out how we might relax some of these assumptions, and generalize eq. \ref{clogistica}. 

First let's realize that the density independent per capita growth rate, $r$, is really a net rate --- it is the difference between the density independent, instantaneous birth and death rates, $b$ and $d$. If we start with per capita logistic growth, the generalization would look like this.
\begin{align}
\label{clogisticderive1}
\frac{\D{N}}{N\D{t}} &= r\left(1-\alpha N\right)\\
 &= \left(b - d\right)\left(1-\alpha N\right)\\
 &= b -  d  - b\alpha N  + d\alpha N
\end{align}
where $b-d=r$. 

Note that the positive, density-independent effect on growth rate of $b$ is counterbalanced a bit by a negative density-dependent effect of $b$. As $N$ increases, increasing births tend to rein in growth rate a bit, because more births mean a larger negative density-dependent effect. Similarly, mortality, $d$ includes a small \emph{positive} density-dependent effect  that helps enhance growth rate because it reduces the negative effect of $\alpha N$ --- death frees up resources.

Density dependence is the effect of density on growth rate, and thus far we have let that be $1-\alpha N$. Let's represent this as the function $F\left(N\right) = 1 -\alpha N$. We can now note that $\alpha$ is also a net effect --- $N$ will have separate effects on the birth and death rates, and $\alpha$ could be just the sum of more arbitrary constants (e.g., $\alpha N = (x+y)N$);  $1-\alpha N$ is merely the simplest form of $F\left(N\right)$. We can generalize further and let  density affect the birth and death rates separately,
\begin{gather}
  \label{eq:ddd}
\frac{\D{N}}{N\D{t}} = rN\,F\left(N\right)\\
  F\left(N\right) = B\left(N\right) + D\left(N\right)
\end{gather}

We might anticipate that density has no effect on death rates when $N$ is low, because there is no particular limiting factor (Fig. \ref{fig:dd1}). Paradoxically, however, when $N$ becomes large, the ensuing mortality benefits growth rate a bit because death frees up limiting resources (Fig. \ref{fig:dd1}). Such an idea might have the simple form
\begin{equation}
  \label{eq:7}
  D\left(N\right) = gN^2
\end{equation}
where $g$ is a constant. This function starts out at zero (no effect) and increases rapidly  (Fig. \ref{fig:dd1}). This means that at low density, mortality does not influence growth rate, whereas at high density, mortality enhances growth rate.

The density dependence for birth rates could also be more complicated. The \emph{Allee effect} \cite{Stephens:1999mz} arises when a population gets very small and mating success declines because mates can't find each other, leading to a large negative per capita impact of $N$. We can describe this with a quadratic function, 
\begin{equation}
  \label{eq:6}
  B\left(N\right) = - aN^2 + eN - f
\end{equation}
where $a$, $e$, and $f$ are constants. Here, the negative effect of $N$ on birth rates, $B(N)$, will be important at low $N$ because when $N=0$, $B(N)=-f$. The negative effect will also be large at very high $N$ because $-aN^2$ becomes a large negative number at large $N$ (Fig. \ref{fig:dd1}). The point here is that we can alter density dependence to incorporate specific biology.

If we sum the density dependences, we have a function $F\left(N\right)=B(N)+ D(N)$ that is analogous to $1-\alpha N$, but which is a different shape, because it describes a  particular underlying biology (Fig. \ref{fig:dd2}). Note the differences between the simple linear density dependence of the logistic model, and the nonlinear density dependence of our generalization.
\begin{figure}[ht]
  \centering
  \subfloat[Birth and Death Rates]{\includegraphics[width=.47\linewidth]{DD1.pdf}\label{fig:dd1}}
  \subfloat[Net Density Dependence]{\includegraphics[width=.47\linewidth]{DD2.pdf}\label{fig:dd2}}
  \caption{Elaborating on simple logistic growth. \subref{fig:dd1}, a specific example of a generalized density dependence via effects on birth rate and death rate. \subref{fig:dd2} Net density dependence for our generalized model and logisitc linear density dependence.}
\end{figure}

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Density Dependence on Birth and Death Rates (Figs. \ref{fig:dd1}, \ref{fig:dd2})}
To make expressions for eqs. \ref{eq:6}, \ref{eq:7}, use \texttt{expression}.
<<>>=
B.N <- expression(-a*N^2 + e*N - f)
D.N <- expression(g*N^2)
@ 
We then provides constants and evaluate the expressions to plot the density dependence on birth and death rates.
<<>>=
a <-1/1600; e <- 80*a; f <- .2; g <- 0.00004
N <- 0:100
<<DD1, fig=TRUE>>=
plot(N, eval(B.N), type="l", ylab="B(N), D(N)")
lines(N, eval(D.N), lty=2)
abline(h=0, lty=3)
legend("bottom", c("Effect on Birth Rate", "Effect on Death Rate"), lty=1:2, bty='n')
@ 
The sum of these two rates is merely density dependence,
<<results=hide, echo=false>>=
myQ()
<<>>=
plot(N, eval(B.N) + eval(D.N), type='l', ylim=c(-0.4,1), ylab="Net Density Dependence, F(N)")
abline(h=0, lty=3)
@ 
and we can compare this to linear density dependence of the logistic model ($1-\alpha N$). If $\alpha=0.01$, we have
<<>>=
curve(1-0.01*x, 0,100, lty=2, add=T)
legend('topright', c("Generalized", "Logistic"), lty=1:2, bty='n')
<<echo=false, results=hide>>=
dev.print(pdf, "DD2.pdf")
@ 
}\end{boxedminipage} \medskip

The extreme flexibility of eq. \ref{clogisticderive1} has a serious downside---we now have to estimate more parameters. Clearly, if $r=0.5$, then $b$ and $d$ can take on any values, provided their difference is $0.5$. We can make the same argument with regard to linear density dependence as well --- $B(N)$ and $D(N)$ could take an infinite variety of forms, provided their sum equals $1- \alpha N$.

In an attempt to simplify our lives a bit, and in the absence of any additional information, let us impose the simplest function for $B$ and $D$, a linear constant
 \begin{equation}
 \label{clogisticderive2} 
\frac{\D{N}}{N\D{t}} = \left(b -  m_{b}N\right) - \left(d - m_{d}N\right) 
\end{equation} 
where $m_{b}$ and $m_{d}$ are constants that describe the effect of an additional individual on the base density-independent per capita birth and death rates, respectively.  If we really wanted to simplify our lives, we could let $m_b+m_d=m$, $b-d=r$, and simplify to
\begin{equation}
\label{clogisticdiff}
\frac{\D{N}}{N\D{t}} = r - mN
\end{equation}
Thus, the \emph{net} effect of $N$ is linear, even though the effects of $N$ on birth and death rates may not be. For reasons that probably have more to do with history than with mathematics, we create another constant, $\alpha=m/r$, and we find ourselves back to 
\begin{equation}
\label{clogisticderive3}
\frac{\D{N}}{N\D{t}} = r \left( 1 - \alpha N \right).
\end{equation}

Let's review the circular path of elaboration and simplification upon which we have travelled. We started with a historical representation of logistic growth, using $K$ to represent an equilibrium population size toward which $N$ is attracted. We represented this in a tad more mechanistic fashion (I would argue) by rephrasing the limitation factor in terms of a per capita effect of an individual on the growth rate. We then expanded and generalized this expression to see that the \emph{net} effects, $r$ and $\alpha$, could each be thought of as two components related to birth rates and death rates. We further expanded on density dependence, allowing it to be curvilinear. We then collapsed this back to a simple difference ($r-mN$, eq. \ref{clogisticdiff}) and finally recast it in the original expression.

\subsubsection{The integral of logistic growth}
\index{logistic growth!integral}
We might as well add that there is an expression for $N_{t}$, analogous to $N_{t}=N_{0}\exp(rt)$ for exponential growth. If we integrate eq. \ref{clogisticderive3}, we have
\begin{equation}
N_{t}=\frac{N_{0}e^{rt}}{1+\alpha N_{0}\left(e^{rt}-1\right)}.
\label{fig:clogisticNt}
\end{equation}
Note the resemblance to the exponential equation. We have the exponential equation in the numerator, divided by a correction factor which starts at 1 when $t=0$. The correction factor then grows with time, at a rate which is virtually the same as exponential growth but adjusted by $\alpha$ and scaled by $N_{0}$. Thus when $t=0$, and if $N_0$ is small, we have near-exponential growth, decreased merely by a factor of $\alpha N_0$. We will take advantage of this integral below, when we describe some data in terms of the logistic growth model.


\subsection{Equilibria of the continuous logistic growth model}
\index{logistic growth!equilibrium, stability, and dynamics}
As we did with the discrete model, we can determine the \emph{equilibria} of the continuous logistic growth model, that is the values of $N$ for which $\D{N}/\D{t}=0$. Consider again eq. \ref{clogisticderive3} 
\begin{equation}
\label{clogistica2}
\frac{\D{N}}{\D{t}} = rN \left( 1 - \alpha N \right).
\end{equation}
What values of $N$ will cause the growth rate to go to zero, thus causing the population to stop changing? As with the discrete model, we would find that if $N=0,\,1/\alpha$, the growth rate would be zero, and the population stops changing. Therefore the equilibria for continuous logistic growth are $N^*=0,\,1/\alpha$.  

\subsection{Dynamics around the equilibria --- stability}
The \emph{stability} of a system\footnote{A system might be nearly any set of interacting parts, such as a population, or an economy, or the internet.} is a measure of how much it tends to stay the same, in spite of external disturbances or changes in the state of the system. 

The term \emph{stability} has been given many more specific meanings, including resilience, resistant,  reactivity, and permanence. We won't go into these here, but one could consider them different facets of stability \cite{Chen2001}.\footnote{Resistance refers to the tendency to resist change in the first place; resilience refers to the long term rate of return to an equilibrium; reactivity refers to the rate of change immediately following such a disturbance \cite{Neubert1997}, and permanence is the degree to which a system does not simplify by entering one of its boundary equilibria.} We will focus on resilience, the tendency for a population to return an equilibrium, or be \emph{attracted toward} an equilibrium, if it is perturbed from it.

Consider the stability of a marble inside a wok. If the wok doesn't move, then the marble just sits in the lowest spot on the bottom. If the wok is bumped, the marble jiggles and rolls around a bit, but settles back down to the bottom. This is a stable system because there is a point at the bottom of the bowl (the attractor) toward which the marble rolls --- all points inside the wok allow the marble to find the lowest point (the attractor). For this reason, we call the collection of points on the inside surface of the wok the \emph{basin of attraction}. The steeper the sides, the more quickly the marble moves toward the bottom. This rate is referred to as its \emph{resilience}.

To translate the notion of a basin into a population, let the position of marble inside the wok be the value of $N$ at any particular point. Bumping the wok is any disturbance that causes a change in $N$. The slope of the sides of the wok reflects the biology and its mathematics that cause $N$ to move quickly or slowly back toward the bottom. For the logistic model, this rate is determined by $r$ and $\alpha$, and also by $N$ itself. The attractor at the very bottom of the bowl is the population's carrying capacity, $K=1/\alpha$.\footnote{The marble's inertia, which contributes to its momentum, is analogous to a time lag in density dependence that we saw in the discretet model; with a time lag, a population can overshoot the equilibrium.}

When we imagine a marble in a wok, it becomes easy to envision $K$ as an attractor at the bottom of the bowl.  That is why we consider the carrying capacity a stable equilibrium point, or attractor, even if other factors, such as a disturbance, may cause $N$ to deviate from it.

The equilibrium we have focused on has been $K$, but recall that $N=0$ is also an equilibrium. This equilibrium is actually the edge of the wok --- the slightest jiggle, and the marble falls in and moves toward the center of the wok, $K$.  The biological analog of this ``jiggle'' is the additional of one or a very small numbers of individuals to an otherwise extinct population. For example, consider that a sterile petri dish with nutrient agar has an \emph{E. coli} population size of zero. If that \emph{E. coli} $N=0$ gets ``perturbed'' with a single added cell, the population will move quickly toward its carrying capacity $K$. In this sense, $N=0$ is an equilibrium, but it is an \emph{unstable} equilibrium. We also call such an unstable equilibrium a \emph{repeller}.

\subsubsection{Analytical linear stability analysis}
We are interested in the dynamics or stability of $N$ at each of the equilibria, $N^*$. Here we use \emph{analytical linear \index{stability analysis!single species}stability analysis} to show mathematically what we described above with the marble and the wok. While somewhat limited in its scope, linear stability is nonetheless a powerful technique for dynamic model analysis.

In a nutshell, what we do is determine whether the growth rate, which is zero at each equilibrium, becomes positive or negative in response to a small change in $N$. That is, if $N$ takes a tiny step away from the equilbrium, does that tiny step shrink, or grow? If a tiny step results in a shrinking of that step, back toward the equilibrium, that demonstrates stability. If a tiny step results in growth and a bigger step, that demonstrates instability. That is what analytical stability analysis tells us.

Consider a plot of the growth rate, $\D N/\D t$ \emph{vs.} $N$ (Fig. \ref{fig:FN.N}).  The equilibria, $N^*$, are the $N$ ($x$-axis) at which $\D{N}/\D{t}=0$ (points a, d). Note where population growth rate is positive and where it is negative. What will happen to this population if it finds itself at $N=50$? It will grow, moving along the $x$-axis, until population growth rate slows so much that it comes to rest where $N=1/\alpha=K=100$. Thus, $N$ changes until it converges on the equilibrium, $N^*$, where $\D{N}/\D{t}=0$. Alternatively at $N=110$, $\D N/\D t$ is negative, and so $N$ shrinks back down to $K$ (point d). This return toward $K$ means that $K$ is an attractor and a stable equilibrium.

\begin{figure}[ht]
  \centering
  \includegraphics[width=.85\linewidth]{Stability1.pdf}
  \caption{Population growth rate, $\D{N}/\D{t}$, as a function of $N$. Points a--e are labelled for convenience. At values of $N$ associated with points $a$ and $d$, growth rate (here labelled as $F(N)$) equals zero. At values of $N$ associated with $b$ and $c$ growth rate is positive, and for $e$ growth rate is negative. \emph{Note this is growth rate as a function of $N$} (time is not explicitly part of this graph).}
  \label{fig:FN.N}
\end{figure}

Next, consider $N=0$, at point a; it cannot change on its own. However, if ``perturbed'' at all, with the addition of one or more individuals, then this ``perturbation'' will grow, sending $N$ across the $x$-axis, away from $N=0$, toward $N=K$. This stasis at $N=0$, and movement away from $N=0$ with the slightest perturbation means that $N=0$ is a repeller and an unstable equilibrium.

\emph{Linear stability analysis will calcculate the degree of attraction or repulsion.}


\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Growth rate $vs.$ $N$}
 We first define an expression, and constants.
<<results=hide, echo=false>>=
myQ(w=6)
<<>>=
pop.growth.rate <- expression( r * N * (1 - alpha * N) )
r <- 1; alpha <- 0.01; N <- 0:120
@ 
A basic plot.
<<>>=
plot(N, eval(pop.growth.rate), type="l",
     ylab="Population Growth Rate (dN/dt)", xlab="N")
abline(h=0)
legend('topright', "r=1", lty=1)
@ 
Add a few points with labels,
<<>>=
N <- c(0, 10, 50, 100, 115)
points(N, eval(pop.growth.rate), cex=1.5)
text(N, eval(pop.growth.rate), letters[1:5],adj=c(.5,2))
@ 
and arrows for the direction of change in $N$.
<<>>=
arrows(20,2, 80,2, length=.1, lwd=3)
arrows(122,-2,109,-2, length=.1, lwd=3)
<<results=hide, echo=false>>=
dev.print(pdf, "Stability1.pdf")
@ 
}\end{boxedminipage} \medskip

How fast will $N$ return to $N^*=K$, if perturbed? Let's start by imagining that we have two populations with different intrinsic rates of increase ($r=1,\,0.5$), and focus in on the growth rate at $N^*=K$ (Fig. \ref{stab2}). For which one of these populations will $\D{N}/\D{t}$ go to zero more quickly? If each population loses a few individuals and declines by $x$ to $N^*-x$, which population will return to $N^*$ more quickly? Not surprisingly, it is the population with the higher $r$ --- when $r=1$, its population growth rate ($y$-axis) is greater than when $r=0.5$, and therefore will return at a faster rate.

Note also that this same population ($r=1$) has a negative growth rate of greater magnitude at $N^*+x$. Regardless of whether it gains or loses individuals, it will return to $N^*$ more quickly than the population with the lower $r$.

\begin{figure}[ht]
  \centering
  \subfloat[Growth Rate $vs$. $N$]{
  \includegraphics[width=.48\linewidth]{Stability2.pdf}\label{stab2}}
  \subfloat[Exponential Recovery]{%
  \includegraphics[width=.48\linewidth]{Stability3.pdf}\label{stab3}}
  \caption[Linear Stability]{Very close to an equilibrium, a population recovers from a perturbation, moving toward the equilibrium attractor, at rate $e^{-r}$. \subref{stab2}, Slopes of population growth $vs.$ $N$ near the equilibrium; around the equilibrium attractor, small decreases in $N$ lead to positive growth rate, whereas small increases lead to negative growth rate; the population with the steeper slope changes faster. \subref{stab3}, Regardless of the particular $r$ of a population, a population recovers from a perturbation, moving toward the equilibrium attractor, at rate $e^{-r}$. }
\end{figure}

How do we quantify the rate of return at $N^*$ so that we can compare two or more populations? The slope of the curve at $N^*$ quantifies the rate of return, because the slope is the exponential return rate. To calculate the slope of a curve we use calculus, because the slope of a curve is a derivative. In this case, we need the derivative of the growth rate ($\D N/\D t$) with respect to the variable on the $x$-axis, $N$. Such a derivative is called a \emph{\index{partial derivative}partial derivative}. Therefore we need the partial derivative of growth rate, with respect to $N$. This will be the slope of the growth rate with respect to $N$.

To find the partial derivative, let us denote the growth rate as $\dot{N}$ (``N-dot''),
\begin{equation}
  \label{eq:dot}
  \dot{N} = rN\left(1-\alpha N\right).
\end{equation}
$\dot{N}$ is a common symbol for a time derivative such as a growth rate.

Given $\dot{N}$, its partial derivative with respect to $N$ is
\begin{equation}
  \label{eq:pdlog}
  \frac{\partial \dot{N}}{\partial N} = r - 2r\alpha N.
\end{equation}
Further, we are interested in the equilibrium, where $N = 1/\alpha$. If we substitute this into eq. \ref{eq:pdlog} and simplify, this reduces to 
\begin{equation}
  \label{eq:pdlog2}
  \frac{\partial \dot{N}}{\partial N} = -r \, .
\end{equation}
At $N^*=1/\alpha$, the slope is $-r$. \emph{This negative slope at the equilibrium demonstrates the stability of this equilibrium.} 


Near $N^*$, population growth rate is $-rx$.\footnote{Because $\Delta y/\Delta x = - r$, i.e. the change in $y$ equals the change in $x$ times the slope.} That is, the rate of change of $x$ (Figs. \ref{stab2}, \ref{stab3}) is 
\begin{equation}
  \label{eq:9}
  \frac{\D{x}}{\D{t}} = -rx
\end{equation}
This allows us to pretend that the perturbation will diminish exponentially, because the growth rate is constant. We say ``pretend,'' because we realize that this rate applies only in the small region where we can assume the slope is a straight line. We can describe the size of $x$ at any time $t$ by integrating this well-known differential equation with respect to time as
\begin{equation}
  \label{eq:8}
  x_t = x_0e^{-r t}
\end{equation}
where $x_0$ is the size of the initial displacement relative to the equilibrium, and $x_t$ is the size of the displacement at time $t$ relative to the equilibrium.\footnote{We also did this in Chapter 1, for exponential growth.}

If we choose an $x_t$ carefully, and do so for any and all such analyses, we can determine the time required to reach $x_t$ and we call that time the \emph{characteristic return time}. To determine the characteristic return time, we will let $x_t=x_0/e$, thus defining the characteristic return time as the amount of time required to reduce the perturbation by 63\% (i.e. $1/e$). We can solve for $t$ by dividing through by $x_0$ and taking the $\log$ of both sides,
\begin{align}
  \label{eq:6}
  \frac{x_0}{e} &= x_0e^{-r t}\\
  \log\left(e^{-1}\right) &= -r t\\
  t &= -\frac{1}{\left(-r\right)}.
\end{align}
Thus we see return time, $t$, here is a positive number, with the same units as $r$, and depends on the slope of the the growth rate with respect to $N$. It also depends on the assumption of linearity very close to the equilibrium.


<<Stability2, fig=TRUE, results=hide, echo=false>>=
mi <- 99.4; ma <- 100.6
N <- seq(mi,ma, length=30)
plot(N, eval(pop.growth.rate), type="l",
     ylab="Population Growth Rate (dN/dt)", xlab="N")
abline(h=0)
r <- .5
lines(N, eval(pop.growth.rate), lty=2); r <- 1
segments(c(100-.3, 100+.3), c(-.05,-.05), c(100-.3,100+.3), c(.05,.05))

text(99.7,-.15, quote("N*"-"x"))
text(100.3,.15, quote("N*"+"x"))
points(100, 0)
text(100, -.1, "N*")

legend('topright', c("r=1", "r=0.5"), lty=1:2)
<<results=hide, echo=false>>=
myQ( mr = c( 1,  3, .25, 1.25) )
par( mgp = c( 2, .75, 0 ) )
x <- 1
plot(c(-1,8), 1.2*c(-x,x), type='n',
     ylab="N", xlab='', axes=FALSE)
#box()
axis(2, at=c(-x,0,x), labels=c(expression("N*"-"x"),"N*", expression("N*"+"x")) )
abline(h=0)
t <- 1
points(c(0, t), c(1, x*exp(-1*t) ), pch=c(19,1))
curve(1*exp(-1*x), 0, .98, add=T, lty=3)
arrows(0,-x-.2, 7,-x-.2)
mtext("Time", side=1)
text(0, 0, "Perturbation", srt=90, adj=c(-.1,-.1) )  
arrows(0,0,0,.8, length=.1)
text(3, .7, expression("Recovery at rate"=="e"^"-r"))

points(c(0+2, t+2), c(-1, -x*exp(-1*t) ), pch=c(19,1))
x2 <- seq(0,.98, length=30)
y <- 0 - 1 * exp(-1*x2)
lines(x2+2, y, lty=3)
text(2, 0, "Perturbation", srt=90, adj=c(1.02,-.1) )  
arrows(2, 0, 2,-.8, length=.1)
text(5, -.7, expression("Recovery at rate"=="e"^"-r"))

<<results=hide, echo=false>>=
dev.print(pdf, "Stability3.pdf")

@ 
To review, we took the partial derivative of $\D{N}/\D{t}$ with respect to $N$, and then evaluated the partial derivative at $N^*$ to get the rate of change of the perturbation, $x$. A negative value indicates that the perturbation will decline through time, indicating a stable equilibrium or attractor. A positive value would indicate that the perturbation will grow over time, thus indicating an unstable equilibrium, or a repellor.

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
    \paragraph{\index{symbolic differentiation}Symbolic differentiation}
We can use \R's minimal symbolic capabilities to get derivatives. Here we get the partial derivative and evaluate for the two equilibria (Fig. \ref{fig:FN.N}).
<<>>=
dF.dN <- deriv(pop.growth.rate,"N")
N <- c(0, 1/alpha) 
eval(dF.dN)
@
The first value, 1, corresponds to the first value of $N$, which is 0. Because it is positive, this indicates that the perturbation will increase with time, meaning that $N=0$ is a repellor. The second value, $-1$, is negative, and so indicates that the perturbation will decrease with time, meaning that $N=1/\alpha$ is an attractor. 
}\end{boxedminipage} \medskip

\subsection{Dynamics}
How does the continuous logistic growth equation behave? It is very well behaved, compared to the discrete version. It has a very characteristic shape (Fig. \ref{c1}), and regardless of $r$ or $N_0$, it moves converges boringly toward $K=1/\alpha$ (Fig. \ref{c2}). 

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
  \paragraph{Function for an ODE}
 Making a function to use with \R's ODE solver is pretty easy, provided we follow the rules (see Appendix, secs. \label{sec:writing-your-own}, \ref{sec:numer-integr-ordin}). To make the code transparent, I translate the vector parameters and the vector of populations (in this single species model, we have only one population).
<<>>=
clogistic <- function(times, y, parms) {
  n <- y[1];  r <- parms[1];  alpha <- parms[2]
  dN.dt <- r * n * (1 - alpha * n)
  return( list( c(dN.dt) ) )
}
@ 
We create vectors for the parameters and the initial densities for all of the populations in the model. We also specify the time.
<<>>=
prms <- c(r=1, alpha=.01); init.N <- c(1)
t.s <- seq(0.1, 10, by=.1)
@ 
We load the \texttt{deSolve} library, and run \texttt{ode}. The output is a matrix that includes the time steps and the $N$ (Fig. \ref{c1}).
<<>>=
library(deSolve)
out <- ode(y=init.N, times=t.s, clogistic, parms=prms)
<<ContDyn1, fig=TRUE, include=FALSE>>=
plot(out[,1], out[,2], type='l', xlab="Time", ylab="N")
@ 
}\end{boxedminipage} \medskip

@
\begin{figure}[ht]
  \centering
  \subfloat[One Population]{\includegraphics[width=.48\linewidth]{ContDyn1} \label{c1} }
  \subfloat[Many Different Populations]{\includegraphics[width=.48\linewidth]{ContDyn2} \label{c2} }
  \caption{Dynamics of continuous logistic population growth. \subref{c1} The characteristic shape of logistic growth. \subref{c2} Regardless of $N_0$ or $r$, populations converge slowly (small $r$) or quickly (high $r$) on $K=1/\alpha$.}
\end{figure}

So why is the continuous version so boring, while the discrete version is so complex? Remind yourself what is going on with the discrete version. The model could only take steps from one generation to the next. The step into generation $t+1$ was a function of $N$ at $t$, and not a function of $N_{t+1}$. Therefore the rate of change in $N$ was not influenced by a contemporaneous value of $N$. There was a delay between $N$ and the effect of $N$ on the population growth rate. For many organisms, this makes sense because they undergo discrete reproductive events, reproducing seasonally, for instance. In contrast, the continuous logistic population growth is always influenced by a value of $N$ that is updated continuously, instantaneously. That is the nature of simple differential equations. They are instantaneous functions of continuous variables. We can, if we so choose, build in a delay in the density dependence of continuous logistic growth. This is referred to as ``delayed density dependence'' or ``time-lagged logistic growth''

\begin{equation}
\label{clogisticDDD}
\frac{\D{N}}{\D{t}} = rN \left( 1 - \alpha N_{t-\tau} \right)
\end{equation}
where $\tau$ is the degree of the delay, in time units associated with the particular model. Just as with the discrete model, the dynamics of this continuous model can get very complicated, with sufficient lag and sufficient $r$.

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
  \paragraph{Plotting Random Populations ( (Fig. \ref{c2}))}
 We use the above function to create 20 populations with different traits. We start with an empty matrix, and then for each of the populations, we draw random $N_0$ and $r$, run the ODE solver, keeping just the column for $N$. Last we plot the output.
<<>>=
outmat <- matrix(NA, nrow=length(t.s), ncol=20)
for(j in 1:20) outmat[,j] <- { y <- runif(n=1, min=0, max=120)
                          prms <- c(r=runif(1, .01,2), alpha=0.01 )
                           ode(y, times=t.s, clogistic, prms)[,2] }
<<ContDyn2, fig=true, include=false>>=
matplot(t.s, outmat, type='l', col=1, ylab="All Populations")
@ 
}\end{boxedminipage} \medskip

\section{Other Forms of Density-dependence}

There are commonly used types of density-dependence other than logistic. Some, like the Richards model, are more flexible and have more parameters. Others, like the Gompertz model, are used more widely in other fields, such as cancer research for tumor growth, von Bertlanaffy for body size growth, and the Ricker model for fisheries; the Richards model provides even more flexibility, at the cost of more paramteters. Here we explore a simple extension of the logistic model, the \emph{\index{logistic growth!theta-logistic}theta-logistic} model, \index{$\theta$-logistic|see{logistic growth}}which adds a parameter to increase flexibility and generality.
\begin{equation}
  \label{eq:11}
  \frac{\D{N}}{\D{t}} = rN\log \left(1 - \left(\alpha N\right)^\theta \right) 
\end{equation}
Here $\theta$ is strictly positive ($\theta>0$); $\theta=0$ means zero growth, and $\theta < 0$ would mean negative growth below $K$, and unbounded positive growth above $K$. Approximations of eq. \ref{eq:11} that allow $\theta \leq 0$ are possible \cite{Sibly:2005ys}, but the interpretation becomes strained.

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Theta-logistic function}
Here we make a function that we can use with \texttt{ode}, the numerical integration function.
<<>>=
thetalogistic  <- function(times, y, parms) {
  n <- y[1]
  with(as.list(parms), {
    dN.dt <- r * n * (1-(alpha*n)^theta)
    return( list( c(dN.dt) ) ) }  )
}
@
Using \texttt{with()} and \texttt{as.list()} creates an environment in which \R~will look inside \texttt{parms} for named elements. This will work as long as we name the parameters in \texttt{parms}.
}\end{boxedminipage} \medskip

By varying $\theta$, we can change the linear density dependence of the simple logistic model to curvilinear density dependence (Fig. \ref{th1}). This curvilinearity arises because when $\alpha N < 1.0$ (i.e. $0 < N < K$),
\begin{align}
  \label{eq:12}
   \left(\alpha N\right)^\theta &< \alpha N,\quad \theta > 1\\
   \left(\alpha N\right)^\theta &> \alpha N,\quad \theta < 1.
\end{align}
In contrast, when $\alpha N=1.0$, then $\left(\alpha N\right)^\theta < \alpha N$. This means that $\theta$ does not affecct $K$. 

When $\theta > 1$, this weakens density dependence at low $N$, so the population grows faster than logistic, all esle being equal. When $\theta < 1$, this strengthens density dependence at low $N$, causing the population to grow more slowly than logistic, all else being equal.

The effects of $\theta$ on density dependence controls the shape of relation between growth rate $vs$. $N$ (a.k.a. the production function, Fig. \ref{th2}). First, note that for a given $r$, growth rate for $N<K$ increases with $\theta$. Second, note that the position of peak of the production function shifts to the right as $\theta$ increases. That is, as $\theta$ increases, the $N$ at which the maximum growth rate occurs also increases. If we wanted to shift the peak growth rate to a higher $N$ without also increasing the height of the peak, we could decrease $r$ simultaneously. 
 
\begin{figure}[ht]
  \centering
  \subfloat[Density Dependence]{\includegraphics[width=.32\linewidth]{thetaDD.pdf}\label{th1}}
  \subfloat[Growth Rate]{\includegraphics[width=.32\linewidth]{thetaGR.pdf}\label{th2}}
  \subfloat[$N$ Dynamics]{\includegraphics[width=.32\linewidth]{thetaN.pdf}\label{th3}}
  \caption{Theta-logistic growth. In a review of population dynamics, Sibly et al. \cite{Sibly:2005ys} use theta-logistic density dependence, $1-(N/K)^{\theta}$, to show that populations most frequently have a concave-up $\theta < 1$ pattern. }
\end{figure}

We could speculate on biological meaning of $\theta$ and the shape of the denisty dependence. For instance, very high $\theta$ suggests a threshold, wherein the population exhibits little density dependence until very close to $K$. Perhaps this is related to territoriality, or a spatial refuge from predators. Alternatively, low $\theta$ might suggest resource preemption, wherein a small number of individuals can sequester large amounts of resources, but increasing $N$ results in reduced preemption. For organisms with very plastic body sizes (plants, fish), this could mean that at low $N$, average body size is large, but as $N\rightarrow K$, average body size decreases. While $\theta$, \emph{per se}, is devoid of mechanism, discovering the magnitude of $\theta$ for different types of species \cite{Sibly:2005ys} could lead to the generation of new testable hypotheses about what controls populations. 

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Theta-logistic density dependence}
We first graph theta-logistic, for $\theta<1$, $\theta=1$, and $\theta>1$ (Fig. \ref{th1}).
<<ThetaDD, fig=TRUE>>=
r <- .75; alpha <- 0.01; theta <- c(.5, 1, 2); N <- 0:110
theta.out <- sapply(theta, function(th) {1-(alpha*N)^th})
matplot(N, theta.out, type='l', col=1)
abline(h=0)
legend('topright', legend=paste("theta =", c(2, 1, 0.5) ), lty=3:1)
@ 
}\end{boxedminipage} \medskip


\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
  \paragraph{Theta-logistic growth rate}
We plot the growth rate (a.k.a. the production function) for the  theta-logistic model with $\theta<1$, $\theta=1$, and $\theta>1$ (Fig. \ref{th2}).
<<results=hide, echo=false>>=
myQ()
<<>>=
thetaGR.out <- sapply(theta, function(th) {r*N*(1-(alpha*N)^th)})
matplot(N, thetaGR.out, type='l', col=1)
abline(h=0)
legend('bottomleft', legend=paste("theta =", c(2, 1, 0.5) ), lty=3:1)
@ 
We also add an example of growth with low $\theta$, but higher $r$.
<<eval=false, echo=false>>=
r <- 3.4
lines(N, r*N*(1-(alpha*N)^theta[1]))
text(12, 24, "r=3.4")
<<results=hide, echo=false>>=
dev.print(pdf, "ThetaGR.pdf")
@ 
}\end{boxedminipage} \medskip


\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
  \paragraph{Theta-logistic dynamics}
We solve numerically for $N$, and plot the dynamics for $N$ with $\theta<1$, $\theta=1$, and $\theta>1$ (Fig. \ref{th3}).
<<ThetaN, fig=TRUE>>=
prms <- c(r <- 0.75, alpha <- 0.01, theta=1)
thetaN <- sapply(theta, function(th) {prms["theta"] <- th
                                      ode(y=1, t.s, thetalogistic, prms)[,2] } )
matplot(t.s, thetaN, type='l')
legend('topleft', legend=paste("theta =", c(2, 1, 0.5) ), lty=3:1, bty='n')
@ 
}\end{boxedminipage} \medskip

\section{Maximum Sustained Yield}
The classical model of harvesting fisheries and wildlife populations is based on a particular conceptualization of \emph{\index{maximum sustained yield}maximum sustained yield} (\emph{MSY}). Maximum sustained yield is historically defined as the population size where both population growth rate and harvest rate are at a maximum.  For the logistic growth model, this is half the carrying capacity. We can solve this by finding the maximum of the production function. To do this we use calculus (differentiation) to find where the slope of the production function equals zero. Starting with the partial derivative with respect to $N$ (eq. \ref{eq:pdlog}), we solve for zero.
\begin{align}
  \label{eq:msy1}
    0 &= r - 2r\alpha N\\
    N &= \frac{r}{2r\alpha} = \frac{K}{2}
\end{align}
Similarly, we could determine this peak for the $\theta$-logistic model as well. 
\begin{align}
  \label{eq:msyth}
  \frac{\partial \dot{N}}{\partial N} &=r - \left(\theta+1\right)r\left(\alpha N\right)^\theta\\
0 &= \frac{K}{\left(\theta+1\right)^{1/\theta}}
\end{align}
When $\theta=1$, this reduces to $K/2$, as for the logistic model. As we saw above, $\theta < 1$ causes this value to be less than $K/2$.

Even if we assume the logistic model is appropriate for a given population, harvests in practice seldom reduce populations to $K/2$, because often $N$ and $K$ vary with time and are difficult to estimate.  Moreover, there are economic forces that often lead to over-harvesting.  More contemporary models of population \index{harvesting, fisheries}harvesting incorporate these concepts \cite{Roughgarden:1996fr}. We should note that when more detailed information is available on age- or stage-specific survivorship and reproduction, age- or stage-structured models of harvesting are the preferred method (e.g., Chapter 2). For many populations, such as marine fisheries, however, data are often  available only on stock size ($N$) and total catch ($H$).  In some fisheries models, the total catch is fixed and does not vary with $N$.  This leads to highly unstable population dynamics and fortunately is no longer practiced in most fisheries.  Usually the total catch $H$ is modeled as a function of $N$. 

Let us assume that total catch or harvest, $H$, is a simple linear function of $N$; the harvest rate $FN$ is usually described this way.\footnote{We could, but will not, further describe this by $F=qE$, where $q$ is the catchability coefficient for the fishery, and $E$ is fishing effort in hours.} In this case, a maximum sustained yield, \emph{MSY}, can be determined from the logistic growth-harvest model. This basic model is then
\begin{align}
  \label{eq:FN}
dN/dt&=rN(1-\frac{N}{K})-FN
\end{align}
where $FN=H$ harvest rate, and $F$ is the per capita fishing mortality rate. Here we assume a fixed per capita catch mortality $F$ so that total harvest is $H=F\times N$. \emph{MSY} occurs at $K/2$ (Fig. \ref{msy2}). Therefore would like to know the value of $F$ for which a new $N^* = K/2$ rather than $K$. We can determine the value of $F$ for which this is true by finding when $\D{N}/\D{t}=0$ and $N=K/2$.
\begin{align}
  \label{eq:F}
  0 &= r\frac{K}{2}\left(1-\frac{\frac{K}{2}}{K}\right)-F\frac{K}{2}\\
  F &= \frac{r}{2}
\end{align}


\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{\emph{MSY} and harvesting (Fig. \ref{msy1})}
Here we illustrate the interaction harvesting at a rate associated with \emph{MSY} for logistic growth. We set logistic model parameters, and first plot logistic growth without harvest.
<<results=hide, echo=false>>=
myQ()
<<>>=
r <- .5; alpha <- 0.01
N <- 0:105
plot(N, eval(pop.growth.rate), type="l", ylim=c(-3, 15),
     ylab="dN/dt and FN")
abline(h=0)
@ 
We then calculate $F$ based on our solution eq. \ref{eq:5}, and plot the linear harvest function with an intercept of zero, and slope of $F$.
<<>>=
F <- r/2
abline(c(0,F), lty=2)
<<results=hide, echo=false>>=
dev.print(device=pdf, file="MSY.pdf")
@ 
}\end{boxedminipage} \medskip

\medskip \noindent
\begin{boxedminipage}{\linewidth}{\footnotesize
\paragraph{Equilibrium solution for logistic growth with harvesting (Fig. \ref{msy2})}
When we add harvesting at rate $F=r/2$, we get a new equilibrium. Here we illustrate this using the same parameters as above, but now using the entire function with both growth and harvest.
<<MSY2, fig=TRUE>>=
pgr.H <- expression(r * N * (1 - alpha * N)-F*N)
N <- 0:55
plot(N, eval(pgr.H), type="l", ylab="dN/dt (growth - harvesting)")
abline(h=0)
@ 
This merely represents the new equilibrium (where $\D N/\D t$ crosses the $x$-axis) with a harvest of rate $F=r/2$. 
}\end{boxedminipage} \medskip

\begin{figure} \centering
\subfloat[Logistic Growth and Harvest Functions]{\includegraphics[width=.47\linewidth]{MSY.pdf}\label{msy1}}
\subfloat[Logistic Growth \emph{with} Harvesting]{\includegraphics[width=.47\linewidth]{MSY2.pdf}\label{msy2}}
\caption{Illustration of maximum sustained yield under harvesting. \subref{msy1} both the logistic growth model eq. \ref{clogisticK} and the linear harvest model eq. \ref{eq:5}, and \subref{msy2} the combined growth model that includes harvesting as a source of mortality. Equilibria occur when the growth rate, $y$, equals zero, crossing the $x$ axis.}
\end{figure}
If harvesting occurs at a rate of $F<r/2$, the slope of the harvest line is shallower (cf. Fig.~\ref{msy1}), and we refer to this as \emph{under-harvesting}. If the rate is $F>r/2$, then the slope is steeper (cf. Fig.~\ref{msy1}), and we refer to this as \emph{over-harvesting}.  Roughgarden and Smith (1996) show that under-harvesting is probably the most ecologically stable, whereas over-harvesting is frequently thought to be the economically optimal approach.  

Typically the economically optimal approach leads to over-harvesting because the marginal income on a fishery is related to \emph{the interest rate on invested income},  $F= (r + \rho)/2$, where $\rho$ is the interest rate (usually $\rho$ is set to 0.05). We can substitute this into eq. \ref{eq:FN} and solve for $N$ (ignoring $N=0$).
\begin{align}
  \label{eq:ecopt}
  0&=rN(1-\frac{N}{K})-\frac{r + \rho}{2}N\\
   N\frac{r}{K}&=r -\frac{r}{2} -\frac{\rho}{2}\\
  N &=\frac{K}{2}\left(1-\frac{\rho}{r}\right).
\end{align}
Thus, unless $\rho=0$, the typical economic optimum $N_{e}^*$ will be smaller than $N*=K/2$.  This makes intuitive economic sense if we consider income from invested profits. The margin of error, however, is much smaller, putting the fisheries population, and the investment of fishing fleets and associated econmoies at greater risk. Roughgarden and Smith \cite{Roughgarden:1996fr} showed that ecological stability of a fishery is likely to beget economic and cultural stability. They show that the supposed ``economic optimal'' harvesting rate was so low that it increased the extinction risk of the harvested population to unacceptably high levels. The collapse of a fishery is, of course, both economically and culturally devastating.


\section{Fitting Models to Data}
\label{sec:data}
Using theory to guide our data collection, and using data to guide our theory is, in part, how science progresses. In this section, we will use the use the theory of simple continuous logistic growth (linear density dependence) to frame our examination of the effects of nutrients on interactions of an alga that is embedded within a simple food web.

The focus of this section is the practice of describing data in terms of deterministic and stochastic models; this has been greatly facilitated by the \R~language and environment \cite{Bolker:2008rr}.

\subsection{The role of resources in altering population interactions within a simple food web}
Increasing resource availability may increase or, paradoxically, decrease the populations that use those resources \cite{Riebesell:1974qr, Rosenzweig1971,Tilman1982}. Here we explore how adding resources influences the strength of per capita denisty dependence, $\alpha$, of an alga, \index{Closterium@\emph{Closterium acerosum}}\emph{Closterium acerosum}, embedded in a simple food web \cite{Stevens:2006on}. 

The variety of observed responses to resource availability is sometimes mediated by the configuration of the food web in which a population is embedded. To explore the interaction between food web configuration and edibility, Stevens and Steiner \cite{Stevens:2006on} performed a simple experimental resource enrichment using a food web with one predator (\emph{Colpidium striatum}) and two competing consumers which differed in their edibility; bacteria, which were edible and \emph{Closterium acerosum}, which was inedible. Resources were added at two concentrations, a standard level, and a ten-fold dilution. We were interested in whether the inedible competitor \emph{Closterium} could increase in abundance in the face of competition from bacteria. Although \emph{Closterium} is inedible, it may be unable to increase in response to resource additions because it is competing with many bacterial taxa which appear to also be predator resistant. In addition, many bacteria are able to grow quickly in enriched environments. However, keystone predation theory \cite{Leibold1996} suggested that \emph{Closterium} might prosper because the bacteria which coexist with the predator may be poor competitors, and the bacteria with the greatest capacity for utilizing higher resource levels may be kept in check by a predator which also might increase in abundance. Stevens and Steiner \cite{Stevens:2006on} presented time-averaged responses, perhaps related to $K$, because the predictions of keystone predation theory were focused on equilibrium predictions. 

Here, we will explore the effects of the resource supply rate on two parameters of linear density dependence (logistic growth) $r$, and $\alpha$ of \emph{Closterium acerosum}.

In the rest of this section, we will follow a generally cautious and pragmatic step-wise scheme which we adapt for our purposes here. We will
\begin{itemize}
\item import the data and examine its structure.
  \item look at the population dynamics of each replicate microcosm.
    \item consider whether we have outliers due to sampling.
      \item estimate parameters for one of the populations.
        \item fit logistic models to each population.
          \item fit a model to all populations and test the effects of nutrients on logistic growth parameters.
\end{itemize}

<<echo=false, eval=false, results=hide>>=
library(nlme)
library(lattice)
ClostExp <- groupedData(No.per.ml ~ Day|ID, tmp)
@ 
\subsection{Initial data exploration}
First we load two more libraries, the data, and check a summary.
<<>>=
library(nlme)
library(lattice)
data(ClostExp)
summary(ClostExp)
@ 
Next we plot each replicate of the high and low nutrient treatments.
<<echo=false, results=hide>>=
trellis.device(width=8, height=6, color=FALSE)
<<>>=
xyplot(No.per.ml ~ Day|Nutrients, ClostExp, 
       groups=rep, type='b', 
       scales=list(relation="free"),
       auto.key=list(columns=4, lines=T)
       )
<<echo=false, results=hide>>=
dev.print(pdf, "ClostRaw.pdf")
@ 
\begin{figure}[ht]
  \centering
  \includegraphics[width=.67\linewidth]{ClosteriumRaw.pdf}
  \caption{Raw data for \emph{Closterium} time series. Each line connects estimates of the \emph{Closterium} popoulation in a single replicate microcosm. Each estimate is based on a $\sim 0.3$\,mL sample drawn from the microcosm after it has been mixed. Low nutrient concentration is $1/10\times$ the standard (high) concentration.}
  \label{fig:closraw}
\end{figure}
Right away, we might notice several things (Fig. \ref{fig:closraw}). First, the high nutrient populations achieve higher abundances (note values on the $y$-axes). This seems consistent with the ten-fold difference in resource concentration \cite{Stevens:2006on}. Perhaps the second thing we notice is that the low nutrient populations seem more variable, relative to their own population sizes. Third, several observations seem remarkably high --- which days?
<<>>=
subset(ClostExp, Nutrients=="high" & No.per.ml > 1000)
subset(ClostExp, Nutrients=="low" & No.per.ml > 100)
@ 
What might we make of these observations? Are they mistakes, or important population dynamics, or some combination thereof? There may be several reasonable explanations for these data which we revisit at the end of this section. For now, we will assume that these points are no less accurate than the other values, and so will retain them.

Finally, we might note that the populations seem to decline before they increase (days 0--5, Fig. \ref{fig:closraw}).  These declines might well represent a bit of evolutionary change, under the new environment into which these populations have been introduced \cite{Fussmann:2005uq} from their stock cultures; we might consider removing the first two sample points to avoid confound mechanisms of evolution versus simple growth. In general, however, we should be cautious about throwing out data that do not conform to our expectations --- let us leave these data in, for the time being.

Now we estimate the parameters of deterministic logistic growth, first without, and then with, an explicit consideration of time. 

\subsection{A time-implicit approach}
How might we estimate the parameters without explicitly considering the temporal dynamic? Recall the relation between per capita population growth $vs$. population size,
\begin{equation}
  \label{eq:percapita}
  \frac{\D N}{N \D t} = r - r \alpha N.
\end{equation}
Here we see that per capita growth rate is a linear function (i.e., $y=mx+b$) of population density, $N$. Therefore, all we have to do is calculate per capita growth rate from our data and fit a simple linear regression of those growth rates versus $N$. The $y$-intercept (when $N=0$) will be our estimate of $r$, and the slope will be an estimate of $r\alpha$.

We will estimate per capita growth rate, \emph{pgr}, with
\begin{equation}
  \label{eq:5}
  pgr = \frac{\log \left(N_{t+i}/N_t \right) }{i}
\end{equation}
where $i$ is the time interval between the two population densities \cite{Sibly:2005ys}. In this and other contexts, logarithm of two populations is sometimes referred to simply as the \emph{log ratio}.

Let's estimate this for one of the microcosms, high nutrient, replicate c.
<<>>=
Hi.c <- subset(ClostExp, Nutrients=="high" & rep=="c")
@ 
We calculate the differences in density, the time intervals over which these changes occured, and the per capita growth rate.
<<>>=
n <- nrow(Hi.c)
N.change <- Hi.c$No.per.ml[-1]/Hi.c$No.per.ml[-n]
interval <- diff(Hi.c$Day)
pgr <- log(N.change)/interval
@ 
We can then (i) plot these as functions of $N_t$ (as opposed to $N_{t+1}$), (ii) fit a linear regression model through those points, and (iii) plot that line (Fig. \ref{pgr1}). 
<<results=hide, echo=false>>=
  myQ()
<<>>=
Nt <- Hi.c$No.per.ml[-n]
plot(pgr ~ Nt)
mod1 <- lm(pgr ~ Nt)
abline(mod1)
<<results=hide, echo=false>>=
dev.print(pdf, "PGR.pdf")
@ 
\begin{figure}[ht]
  \centering
  \label{fig:pgr1}
  \includegraphics[width=.67\linewidth]{PGR.pdf}
  \caption{Per capita growth rate $vs$. population size at the beginning of the interval (data from one high nutrient replciate).}
\end{figure}
The linear regression model defines a best fit straight line through these points, 
<<>>=
summary(mod1)
@ 
where $y$-intercept is the first coefficient of the linear regression; this is an estimate of $r$ (eq. \ref{eq:percapita}). The slope of the line is the second coefficient and is $r\alpha$ (eq. \ref{eq:percapita}).

Now we can do this with all the replicates. Rearranging the data is a bit complicated, but this will suffice. Essentially, we split up the \texttt{ClostExp} data frame into subsets defined by nutrient and replicate ID, and then apply a function that performs all of the steps we used above for a single microcosm. The output will be a \emph{list} where each component is a new small data frame for each microcosm.
<<>>=
EachPop <- lapply( split(ClostExp, list(ClostExp$Nutrients, ClostExp$rep)),
                       function(X) {
                        n <- nrow(X) 
                        N.change <- (X$No.per.ml[-1]/X$No.per.ml[-n])
                        interval <- diff(X$Day)
                        data.frame(Nutrients=as.factor(X$Nutrients[-n]), 
                              rep = as.factor(X$rep[-n]), 
                              pgr = log(N.change)/interval,
                               Nt = X$No.per.ml[-n])
                      })

@ 
Next we just stack all those individual data frames up into one.
<<>>=
AllPops <- NULL
for(i in 1:length(EachPop) ) AllPops <- rbind(AllPops, EachPop[[i]])
@ 
Finally, we plot all of the microcosms separately.
<<results=hide, echo=false>>=
trellis.device(width=8, height=6, color=FALSE)
<<>>=
xyplot(pgr ~ Nt|rep*Nutrients, AllPops, layout=c(4,2,1), 
       #scales=list(x=list(rot=90, alternating=F, relation="free")),
       scales=list(x=list(rot=90)),
       panel=function(x,y){
         panel.grid(h=-1, v=-4)
         panel.xyplot(x,y, type=c('p','r') )
         })
<<results=hide, echo=false>>=
dev.print(pdf, "pgrall.pdf")
@ 
\begin{figure}[ht]
\centering
\includegraphics[width=\linewidth]{pgrall.pdf}
\caption{Per capita growth rates $vs$. population size, $N$, at the beginning of the interval for which the growth rate is calculated (includes all microcosms). Note that $y$-intercepts are all fairly similar, while the maximum $N$, slopes and $x$-intercepts differ markedly. }
\label{fig:pgrall}
\end{figure}
The biggest difference about these populations appears to be that the carrying capacities differ between the high and low nutrients (Fig. \ref{fig:pgrall}). First note that the $y$-intercepts are all remarkably similar, indicating that their $r$ are all very similar. The big differences between the high and low nutrients are that the slopes of the relations and the $x$-intercepts are both very different. The slope is $r\alpha$, while the $x$-intercept is $K$. Next we fit a statistical model to these data.

\subsubsection{A statistical model of per capita growth}
A \emph{\index{mixed model}mixed model}\footnote{This section is advanced material. You may want to skim over the parts you don't understand, and focus on conclusions, or perhaps the general process.} is a statistical model that includes both random effects (e.g., blocks, or subjects) and fixed effects (e.g., a treatment) \cite{Pinheiro2000}. Let's not worry about the details right now, but we'll just use it to test for differences in slope, where we have measure the same microcosms repeatedly through time.

Previously, we showed that it makes at least a modicum of sense to fit eq. \ref{eq:percapita} to each microcosm (Fig. \ref{fig:pgrall}).  Having done that, we next fit all of these to a single mixed model, in which the microcosm is our \emph{block} or \emph{subject} which may exert some unknown random effect on the populations. First we make sure we have a unique ID label for each microcosm.
<<>>=
AllPops$ID <- with(AllPops, Nutrients:rep)
modSlope <- lme(pgr ~ Nt + Nutrients + Nt:Nutrients, data=AllPops, 
                random=~1|ID )
@ 
The above code specifies that per capita growth rate, \texttt{pgr}, depends in part upon population size, \texttt{Nt} --- this will be the slope of the \texttt{pgr~Nt} relation. The effect of nutrient level, \texttt{Nutrients}, will simply alter the intercept, or height, of the lines. The code also specifies that the slope the relation can depend upon the nutrient level, \texttt{Nt:Nutrients}. Finally, the code also specifies that each microcosm makes its own unique and \emph{random} contribution to the per capita growth rate of its \emph{Closterium} population. 

In a proper statistical analysis, it would be necessary to first check the assumptions of this model, that the noise (the residuals) and the random effects are both normally distributed. A plot of the residuals $vs$. the fitted or predicted values should reveal a scatter of points with relatively constant upper and lower bounds.
<<echo=false, results=hide>>=
trellis.device(width=4, height=4, color=FALSE)
<<>>=
xyplot(resid(modSlope) ~ fitted(modSlope))
<<echo=false, results=hide>>=
dev.print(pdf, "pgrres1.pdf")
@ 
Our plot (Fig. \ref{pgr1}) reveals larger residuals for larger fitted values --- this is not good.
\begin{figure}[ht]
  \centering
  \subfloat[Residuals $vs$. Fitted]{\includegraphics[width=.47\linewidth]{pgrres1.pdf}\label{pgr1}}
  \subfloat[Random Effects Q-Q Plot]{\includegraphics[width=.47\linewidth]{pgrranqq.pdf}\label{pgr2}}\\
  \subfloat[Residuals Q-Q Plot]{\includegraphics[width=\linewidth]{pgrres2.pdf}\label{pgr3}}
  \caption{Diagnostic plots for fitting the per capita growth rate to population size.}
  \label{fig:diagnostic}
\end{figure}
Quantile-quantile plots of the residuals and random effects should reveal points along a straight diagonal line. 
<<echo=false, results=hide>>=
trellis.device(width=8, height=3, color=FALSE)
<<>>=
qqmath( ~resid(modSlope)|ID, data=AllPops,layout=c(8,1,1) )
<<echo=false, results=hide>>=
dev.print(pdf, "pgrres2.pdf")
<<echo=false, results=hide>>=
trellis.device(width=4, height=4, color=FALSE)
<<>>=
qqmath(~ranef(modSlope))
<<echo=false, results=hide>>=
dev.print(pdf, "pgrranqq.pdf")
@ 
Our random effects (Fig. \ref{pgr2}) are very small, indicating only small systematic differences among microcosms that are not accounted for by the treatment. Nonetheless, their scatter of points along the diagonal are fairly consistent with a normal distribution. The residuals (Fig. \ref{pgr3}) for each microcosm seem fairly normally distributed except that several populations have outliers (replicates low:a-c, high:c). Taken as a group (Fig. \ref{pgr1}), we see that the residual variation increases the mean (fitted). This is at least in part due to having more observations at larger values, but is also a typical and common problem --- variables with small absolute values simple cannot vary as much as variables with large absolute value. Therefore, we specify that that model takes that into account by weighting residuals observations differently, depending upon the expected (i.e. fitted or mean) value.
<<>>=
modSlope2 <- update(modSlope, weights=varExp())
@ 
We then compare these two models,
<<>>=
anova(modSlope, modSlope2)
@ 
and find that the slightly more complicated one (\texttt{modSlope2} with 7 degrees of freedom) is significantly better (likelihood ratio = 7.1, $P < 0.05$) and is more parsimonious (smaller \index{AIC}AIC) \cite{Burnham:2002mi}.

Now we provide a preliminary test of the fixed effects (nutrients and population size) with analysis of variance. 
<<>>=
anova(modSlope2)
@ 
The first thing we see in this \emph{analysis of variance} output is that the effect of population size, $N_t$, depends upon nutrient level (\texttt{Nt:Nturients}, $P < 0.0001$). Given this, we should avoid drawing conclusions based on the main effects of nutrients. Next we go on to estimate intercepts and slopes associated with different treatments.
<<>>=
summary(modSlope2)$tTable
@ 
Given \R's default contrasts (linear model parameters), this output indicates that the intercept (our estimate of $r$) for the high nutrient group is $\sim 0.13$, and the slope of the high nutrient microcosms is  $-0.000286$. The intercept for the low nutrient level is calculated as $0.13-0.07$, and we find it is not significantly smaller ($P > 0.05$). On the other hand, the slope of the low nutrient microcosms is $-0.000278 - 0.00573$, which is significantly smaller than the slope of the high nutrient microcosms ($P < 0.05$). This information can guide another take on these data.
<<>>=
cfs <- fixef(modSlope2); cfs
- cfs[2]/cfs[1]
- (cfs[2]+cfs[4])/(cfs[1]+ cfs[3])
@ 
indicating that the per capita effects are much smaller in the high nutrient treatment. (The label \texttt{Nt} is an unimportant artifact, retained from using the coefficients to calculate the values.) Next, we use these estimates to help us fit a times series model.

\subsection{A time-explicit approach}
Here we use the above information to help us fit the time-explicit version of the logistic growth model. In this section, we 
\begin{enumerate}
\item create and examine a nonlinear logistic model, 
  \item fit the model to each microcosm using \texttt{nlsList} to confirm it makes sense, and
    \item  fit a single nonlinear mixed model that tests which logistic growth parameters varies between treatments.
\end{enumerate}
At each step, we examine the output and diagnose the fit; it should give you a very brief introduction to this process. For more thorough treatments, see \cite{Bolker:2008rr, Pinheiro2000}.

First, we code in \R~a model of eq. \ref{fig:clogisticNt}, that is, logistic growth with parameters $\alpha$, $r$, and $N_0$. 
<<>>=
ilogistic <- function(t, alpha, N0, r) {
 N0*exp(r*t) / (1 + alpha*N0*(exp(r*t) - 1)) }
@ 
Let's make sure this makes sense by plotting this curve on top of the raw data, focusing on the raw data and model parameters for the high nutrient microcosms.
<<echo=false, results=hide>>=
myQ()
<<>>=
plot(No.per.ml ~ Day,  ClostExp, subset=Nutrients=="high")
curve(ilogistic(x, alpha=-cfs[2]/cfs[1], N0=6, r=cfs[1] ), 
      1, 60, add=T, lty=2)
<<echo=false, results=hide>>=
dev.print(pdf, "hi.pdf")
@ 
\begin{figure}[ht]
  \centering
  \includegraphics[width=.67\linewidth]{hi.pdf}
  \caption{High nutrient microcosms and the logistic growth model using coefficients from \texttt{modSlope}.}
  \label{fig:hi}
\end{figure}
Our model and the coefficients we generated above seem to be a reasonable approximation of the data (Fig. \ref{fig:hi}). Now let's use \texttt{nlsList} to fit the logistic model to each microcosm separately.
<<fit.list>>=
Cmod.list <- nlsList(No.per.ml ~  ilogistic(Day, alpha, N0, r)|ID,
                  data=ClostExp,
                         start= c(alpha=1/100,  N0=5,  r=.15),
                         control=list(maxiter=1000))
@ 
We note that models cannot be fit for three microcosms. For now, we plot the coefficients for each microcosm for which a model was fit.
<<"q", echo=false, results=hide>>=
quartz(, 6,3)
<<>>=
plot(coef(Cmod.list), layout=c(3,1,1))
<<echo=false, results=hide>>=
dev.print(pdf, 'coefs.pdf')
@ 
\begin{figure}[ht]
  \centering
  \includegraphics[width=.75\linewidth]{coefs.pdf}
  \caption{Coefficients estimated for each microcosm separately. Missing coefficients occurred when model fits could not be easily acheived.}
  \label{fig:coefs}
\end{figure}
We can see that the coefficients (Fig. \ref{fig:coefs}) are somewhat similar to those we estimated above with \texttt{modSlope2}. We note, however, that for one replicate $\alpha<0$, which does not make sense. If we check the residuals,
<<"resids.Comd.list">>=
plot(Cmod.list)
@ 
we see that the residuals increase dramatically with the fitted values (graph not shown). This is consistent with our previous observations. 

In this next step, we create a single nonlinear mixed model using \index{nlme@\texttt{nlme}}\texttt{nlme}. We use \texttt{pdDiag} to specify that each microcosm has its own unique effect on each of the parameters, aside from nutrient concentrations. We use \texttt{weights} to specify that the variance of the residuals is a power function ($a^{2t}$) of the mean, and which fits an extra parameter, $t$.
<<>>=
Cmod.all <- nlme( Cmod.list , random=pdDiag(form=alpha + N0 + r~1), weights=varPower()  )
@ 
We saw above that the nutrients caused a big difference in the slopes of the \emph{pgr} $vs$. $N$ relations. Now we would like to estimate the $alpha$ for each level of nutrients. In estimating these, we can also test whether they are different.

First we capture the fixed effect coefficients of the current model.
<<label="cfs2">>=
cfs2 <- fixef(Cmod.all);cfs2
@ 
Next we update the model, specifying that \texttt{alpha} should depend on nutrient level, that is, estimate one $\alpha$ for each level of nutrient. We specify that the other paramters are the same in both treatments. This will mean that we have one additional coefficent, and we need to provide the right number of starting values, and in the correct order (\texttt{alpha1, alpha2, N0, r}).
<<>>=
Cmod.all2 <- update(Cmod.all, fixed = list( alpha ~ Nutrients, N0 + r ~1),
                start=c(cfs2[1], 0, cfs2[2],  cfs2[3]) )
@ 
It would be good practice to plot diagnostics again, as we did above. Here, however, we push on and check the confidence intervals for our parametersa summary.
<<>>=
summary(Cmod.all2)
@ 
If we look first at the random effects on the parameters, we see that the variance components for \texttt{alpha} or \texttt{r} are vanishingly small, we should therefore remove them from the analysis, and then compare the two models.
<<>>=
Cmod.all3 <- update(Cmod.all2, random=N0~1)
anova(Cmod.all2, Cmod.all3)
@ 
We see that the less complex model is not significantly different ($P > 0.05$) and is more parsimonious (smaller AIC, Akaike's information criterion, \cite{Burnham:2002mi}).

At last, we can examine the confidence intervals to get the estimates  on the intraspecific competition coefficents, $\alpha_{ii}$.
<<>>=
intervals(Cmod.all3)$fixed
@ 
see that the analysis, such as it is, indicates that the confidence intervals for the estimated $\alpha$ do not overlap, indicating that the expected values are very, very different from each other. We could calculate the bounds on $K_{high}$ and $K_{low}$ as well, based on our lower and upper confidence limits on $\alpha_{high},~\alpha_{low}$.
<<>>=
1/intervals(Cmod.all3)$fixed[1:2, c(1,3)]
@ 
This aproximately ten-fold difference in $K$ seems surprisingly consistent with the ten-fold difference in nutrient concentration initially established in the experiment. This ten-fold increase in \emph{Closterium} may be merely an odd coincidence because it implies that all species in the food web should have increased ten-fold, but that was not the case \cite{Stevens:2006on}. On the other hand, it may guide future experiments in keystone predation.

We can plot the fixed effects, and also the added variation due random variation in $N_0$ and $r$ due to the ID of the microcosms.
<<results=hide, echo=false>>=
trellis.device(width=8, height=6, color=FALSE)
<<>>=
plot(augPred(Cmod.all3, level=0:1), #scales=list(y=list(log=10)),
     layout=c(4,2,1), scales=list(y=list(relation="free")))
<<results=hide, echo=false>>=
dev.print(pdf, "dyns.pdf")
@ 
\begin{figure}[ht]
  \centering
  \includegraphics[width=\linewidth]{dyns.pdf}
  \caption{Times series plots augmented with the overall predictions (solid lines for the fixed effects) and predictions for individual microcosms (dashed lines).}
  \label{fig:dyns}
\end{figure}
When we plot our data along the with the model (Fig. \ref{fig:dyns}), we can see the disparity between the predictions and the data --- the source of those big residuals we've been ignoring.

There may be at least a couple issues to consider when we think about the disparities.
\begin{description}
\item[{\bf Sampling error}] We do not have much idea about how accurate or precise our sampling is; this might be remedied in the next experiment by using a mean with multiple samples per time interval.
\item[{\bf Experimental Design}] In this system, 10\% of the medium was replaced each week (to remove waste products and replenish nutrients), and this is likely to introduce variability.
\item[{\bf Dynamics}] Species interactions may drive dynamics, derived from feedback between (i) bacteria recycling nutrients, (ii) bacteria consuming nutrients, (iii) the predator \emph{Colpidium} consuming bacteria, and (iv) \emph{Closterium} consuming nutrients (Chapters 5, 6).
\end{description}
Perhaps the first steps for an improved experiment might include sampling at least twice for each sample period, to at least find out what our sampling error is. It may, or may not, be a problem, but a few trials would let us know. We also might increase the frequency and decrease the magnitude of the medium replacement, and to make sure the timing of the sampling corresponds appropriately with the medium replace. After we take care of these issues, then we can take a harder look at the dynamics.

\subsubsection{Conclusions}
Our results show that lowering nutrient concentration increased the intraspecific interaction strengths, $\alpha_{ii}$, ten-fold. Alternatively, we could express this as a ten-fold change in the carrying capacity $K$ of these populations. The effect of resources on interaction strength is sometimes a heated topic of discussion, and our precise definitions, derived from simple growth models, provides an unambiguous result. Our analysis has also pointed to ways that we might perform better experiments. More generally, by framing our experiment in terms of existing theory, we can interpret our results in a more precise manner, and facilitate better advances in the future \cite{Bolker:2008rr}. This will help us make more constructive contributions to the field.

@ 
\section{Summary}
In this chapter, we motivated density dependence, focusing specifically on discrete and continuous logistic growth which assumes linear density dependence. We explored how algebraic rearrangements can yield useful and interesting insights, and how we might think about generalizing or simplifying logistic growth. We showed how the time lag built into the discrete model results in chaos. We explored ideas of maximum sustained yield, and how we might fit the continuous logistic growth model to real data.

% Problems or Exercises should be sorted chapterwise
\section*{Problems}
\addcontentsline{toc}{section}{Problems}
%
% Use the following environment.
% Don't forget to label each problem;
% the label is needed for the solutions' environment
\begin{prob}\label{prob3:1}
\textbf{Dynamics of an annual plant}\\
(a) Calculate $r_d$ of an annual plant that has a maximum growth rate of $N_{t+1}/N_t=2$ at very, very small population sizes.\\
(b) Calculate the appropriate per capita density dependent effect of an annual plant with a carrying capacity $K$ of 100\,inds$\cdot$m$^{-2}$.\\
(c) Write the appropriate logistic growth equation that incorporates the intrinsic growth of (a) and the density dependence of (b).\\
(d) Graph the 10\,y dynamics ($t=0,\ldots,\,10$) of the annual plant in (a) and (b), starting with $N_0=1$.
\end{prob}

\begin{prob}\label{prob3:2}
\textbf{Dynamics of \emph{E. coli}}\\
(a) Calculate $r$ of \emph{E. coli} that has a doubling time of 30\,min. Express this rate in hours. \\
(b) Calculate the per capita density dependent effect of an \emph{E. coli} culture that grows logistically over a 24 \,h period, and which levels off at a density of 10$^7$\,CFU$\cdot$mL$^{-1}$ (CFU is colony forming units --- for \emph{E. coli} its is equivalent to individuals).\\
(c) Graph the 50\,h dynamics ($t=0,\ldots,\,50$) of the \emph{E. coli} population in (a) and (b), starting with $N_0=1000$.\\
\end{prob}
<<eval=FALSE, echo=FALSE, results=hide>>=
prms <- c(r=log(2)/.5, alpha=10^(-7)); init.N <- c(1000)
t.s <- seq(0.1, 40, by=.1)
library(deSolve)
out <- ode(y=init.N, times=t.s, clogistic, parms=prms)
plot(out[,1], out[,2], type='l', xlab="Time", ylab="N")
@ 

\begin{prob}
\label{prob3:theta}
\textbf{Nonlinear Density Dependence}\\
Sibly et al. \cite{Sibly:2005ys} found that most species have nonlinear and concave-up density dependence. They use the $\theta$-logistic growth model.
(a) Create a theta-logistic continuous growth model for use with the \texttt{ode()} function in the \texttt{deSolve} package.\\
(b) Show that with $\theta=1$, it is identical our function \texttt{clogistic} above.\\
(c) Graph $N$ for $t=0,\ldots ,\,100$ using a few different values of $\theta$ and explain how $\theta$ changes logistic growth.
\end{prob}

\begin{prob}
\label{prob3:Harvest}
\textbf{Harvested Populations}\\
The logistic growth equation and other similar equations have been used and abused in efforts to achieve a \emph{maximum sustained yield} of a harvested population. The immediate goal of maximum sustained yield management practices is to kill only the number of individuals that reduces the population to half of its carrying capacity, assuming that eq. \ref{clogisticK} describes the population growth. Answer the questions below to help you explain why this would be a goal.\\
(a) Find expressions for population growth rates when $N = K/4,\,K/2\,3K/4$ (substitute these values for $N$ in eq. \ref{clogisticK}, and show your work). Which of these results in the highest growth rate? How do these relate to the management of a harvestetd population?\\
(b) Show the derivation of the partial derivative of the continuous logistic growth model, with respect to $N$ (i.e., $\partial \dot{N} /\partial N$). Solve for zero to determine when total population growth rate reaches a maximum. Show your work.\\
(c) What would be the ecological and economic rationale for \emph{not} killing more individuals, and keeping $N>K/2$?\\
(d) What would the consequences be for the population if you \emph{assume} linear density dependence ($1-\alpha/N$), but in fact the population is governed by nonlinear density dependence where $\theta<1$ and $theta>1$ (Figs. \ref{th1}-\ref{th3})?\\
(e) What economic benefit would you gain if you harvested the entire population all at once (and eliminated it from the face of the planet)? What could you do with all that money?\\
(f) How would you incorporate both harvesting and economic considerations into your logistic growth model?
\end{prob}

\begin{prob}
\label{prob3:EnvVar}
\textbf{Environmental Variability}\\
Most environments change continually. Temperature, resource availability, changes in predator or pathogen abundance all influence the carrying capacity of the environment. \\
(a) Use the discrete version of the logistic growth equation to model a population in a variable environment. Do this by creating a discrete logistic growth function that adds (or subtracts) a random amount to $K$ in each time step. Use one of the many functions that can draw random numbers from particular distributions (e.g., \texttt{rpois()}, \texttt{rlnorm()}, \texttt{runif()}). You might start by playing with one of the random number generators:\\ \texttt{Kstart <- 100; time <- 1:20; K <- numeric(20); \\
  for(i in 1:20) K[i] <- Kstart + rnorm(1, m=0, sd=5); \\
  plot(time, K)}.\\
(b) All distributions are characterized by their \emph{moments}. For instance, the Normal distribution is typically defined by its mean, $\mu$, and variance, $\sigma^2$. Focus on just one moment of your selected distribution in (a), and use your simulations to determine quantitatively the relation between $K$ and the resulting $N$ derived from the discrete growth model. For instance, you might vary the standard deviation of the random normal distribution that contributes to $K$, and examine how the standard deviation of $K$, $\sigma_{K}$ relates to mean $N$, $\mu_{N}$.\\
(c) Create a reddened time series for $K$.%
\footnote{Environmental factors such as temperature frequently vary in a gradual fashion, such that if the weather is hot today, it is likely to be hot tomorrow. Such variation is described in terms of auto-correlation or a spectral distribution or color \cite{Halley1996,Petchey2001}.  Spectral decomposition of a times series involves a description of the series in terms of waves of different wavelengths or frequencies. White noise is variation that is described by equal contributions of all wavelengths (hence the term ``white''), and it is a series that is completely random from one time step to the next. Reddened environmental variation is variation that is described by a predominance of longer wavelengths (hence the term ``red''), and it is a series in which values tend to change gradually, where values that are close together in time tend to be more similar, or auto-correlated. Spectral variation is also referred to as $1/f$ noise (``one over f noise'').}%
 (Hint: What are the values of x and y when you do \texttt{x <- sample(c(-1,0,1), 20, replace=TRUE); y <- cumsum(x)} ?). Use this time series to create a reddened population dynamic. (Hint: First create the vector of reddened $K$'s equal in length to your time series. Then create a growth function that can access the vector, e.g. \texttt{DLG.RK <- function(alpha=0.001, rd=1, N0=10, gens=20, K=K.red)}. Then, inside the for loop, use an indexed K, such as \texttt{1-N[t]/K[t]}.
\end{prob}