\SweaveOpts{ eval=TRUE, prefix=FALSE, width=4.5, height=4.5, eps=FALSE, include=FALSE} \chapter[An Introduction to Food Webs]{An Introduction to Food Webs, and Lessons from Lotka--Volterra Models} <>= setwd("~/Documents/Projects/Book/draft/Chap07") 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(lattice) trellis.par.set(canonical.theme(color=FALSE)) #source("DMBpplane.R") require(primer); @ A food web is a real or a model of a \emph{set of feeding relations among species or functional groups}. This chapter has two foci, (i) a very brief introduction to multi-species webs as networks, and (ii) a re-examination of old lessons regarding the effect of food chain length on a web's dynamical properties. \begin{figure}[ht] \centering \begin{minipage}[c]{.35\linewidth} \centering \includegraphics[width=3cm]{web1} \end{minipage} \begin{minipage}[c]{.35\linewidth} \centering \begin{tabular}{c| c c c c } & \multicolumn{4}{c}{\emph{~As predators~}}\\ \emph{~As prey~} & ~R~ & ~B~ & ~A~ & ~P~\\ \hline R&0&$-$&$-$&$-$\\ B&$+$&0&0&$-$\\ A&$+$&0&0&0\\ P&$+$&$+$&0&0\\ \end{tabular}\\ \end{minipage} \caption{Two representations of a food web with 4 nodes and 8 directed links. The link label ``v'' indicates the proportion of $P$'s diet comprised of $R$.} \label{fig:web} \end{figure} \section{Food Web Characteristics} \index{food web characteristics} We need language to describe the components of a food web, such as \emph{\index{links}links} and \emph{\index{nodes}nodes}, and we also need language to describe the properties of a web as a whole. Generally speaking, networks such as food webs have \emph{emergent properties} \cite{May:2006lr}, such as the number of nodes in the network. Emergent properties are typically considered to be nothing more than characteristics which do not exist at simpler levels of organization. For instance, one emergent property of a population is its density, because population density cannot be measured in an individual; that is why density is an \index{emergent property}emergent property.\footnote{If we assume no supernatural or spectral interference, then we can also assume that density arises mechanistically from complicated interactions among individuals and their environments, rather than via magic. Other scientists will disagree and say that properties like density that appear to be simple additive aggregates do not qualify for the lofty title of ``emergent property.''} While all food web models are based on simple pairwise interactions, the resulting emergent properties of multispecies webs quickly become complex due to indirect interactions and coupled oscillations \cite{Berlow:2004yq,Vandermeer:2006fj}. In addition, any extrinsic factor (e.g., seasonality) that might influence species interactions may also influence the emergent properties of food webs. A few important \index{network}network descriptors and emergent properties include, \begin{description} \item[\textbf{Node}] A point of connection among links, a.k.a. \index{trophospecies}trophospecies; each node in the web may be any set of organisms that share sufficiently similar feeding relations; in Fig. \ref{fig:web}, $P$ may be a single population of one species, or it may be a suite of species that all feed on both $B$ and $R$. \item[\textbf{Link}] A feeding relation; a connection between nodes or trophospecies; may be directed (one way) or undirected. A directed link is indicated by an arrow, and is the effect ($+,\,-$) of one species on another. An undirected link is indicated by a line (no arrow head), and is merely a connection, usually with positive and negative effects assumed, but not quantified. \item [\textbf{\index{connectance}Connectance}] The proportion of possible links realized. Connectance may be based on either directed, $C_{D}$, or undirected, $C_{U}$, links. For Fig. \ref{fig:web} these would be \begin{gather*} C_{D}=\frac{L}{S^{2}}=\frac{8}{16}=0.5 \\ C_{U}=\frac{L}{\left(S^2-S\right)/2}=\frac{4}{6}=0.67 \end{gather*} where $S$ is the number of species or nodes. \item [\textbf{\index{degree distribution}Degree distribution}, $\mathrm{Pr}(i)$] The probability that a randomly chosen node will have degree $i$, that is, be connected to $i$ other nodes \cite{May:2006lr}. In Fig. \ref{fig:web}, $A$ is of degree 1 (i.e., is connected to one other species). $P$ and $B$ are of degree 2, and $R$ is of degree 3. If we divide through by the number of nodes (4, in Fig. \ref{fig:web}), then the \emph{degree distribution} consists of the probabilities $\mathrm{Pr}\left(i\right)=\{0.25,\,0.5,\,0.25\}$. As webs increase in size, we can describe this distribution as we would a statistical distribution. For instance, for a web with randomly placed connections, the degree distribution is the binomial distribution \cite{Cohen1990d}. \item[\textbf{\index{characteristic path length}Characteristic \index{path length|see{characteristic path length}}path length}] Sometimes defined as the average shortest path between any two nodes \cite{Dunne2002}. For instance, for Fig. \ref{fig:web}, the shortest path between $P$ and $R$ is 1 link, and between $A$ and $P$ is 2 links. The average of all pairwise shortest paths is $(1+1+1+1+2+2)/6=1.\bar{3}$. It is also sometimes defined as the average of \emph{all paths} between each pair of nodes. \item[\textbf{Compartmentation}, $C_I$] The degree to which subsets of species are highly connected or independent of other species.\smallskip To calculate \index{compartmentation}compartmentation in a food web, first assume each species interacts with itself. Next, calculate the proportion of shared interactions, $p_{ij}$ for each pair of species, by comparing the lists of species with which each species in a pair interacts. The numerator of $p_{ij}$ is the number of species with which both of the pair interact. The denominator is the total number of different species with which either species interacts. \medskip As an example, let's calculate this for the above food web (Fig. \ref{fig:web}). $A$ interacts with $A$ and $R$, $B$ interacts with $B$, $R$, and $P$. Therefore, $A$ and $B$ both interact with only $R$, whereas, together, $A$ and $B$ interact with $A$, $B$, $R$, and $P$. The proportion, $p_{ij}$, therefore is $1/4 = 0.25$. We do this for each species pair. \smallskip Next we sum the proportions, and divide the sum by the maximum possible number of undirected links, $C_U$. \begin{table}[h] \begin{minipage}[c]{.48\linewidth} \centering \begin{tabular}{l| c c c } Species & A & B & P\\ \hline R&2/4&3/4&3/4\\ A& &1/4&1/4\\ B& & &3/3 \end{tabular} \end{minipage} \hfill \begin{minipage}[c]{.48\linewidth} \begin{align*} \label{eq:ci} C_{I}&=\frac{\sum_{i=1}^{S-1}\sum_{j=i+1}^{S}p_{ij}}{\left(S^2-S\right)/2}\\ &=3.5/6\\ &=0.58 \end{align*} \end{minipage} \end{table} To reiterate: For any pair of species, $i$ and $j$ ($i\neq j$), $p_{ij}$ is the proportion of shared interactions, calculated from the number of species that interact with \emph{both} species $i$ and $j$, divided by the number of species that interact with \emph{either} species $i$ or species $j$. As above, $S$ is the number of species or nodes in the web. \item [\textbf{Trophic~Level}] \index{trophic position}Trophic \emph{position} may simply be categorized as basal, intermediate or top trophic positions. Basal positions are those in which the trophospecies feed on no other species. The top positions are those in which the trophospecies are fed upon by nothing. One can also calculate a quantitative measure of trophic \emph{level}. This is important in part because omnivory, something rampant in real food webs, complicates identification of a trophic level. We can calculate \index{trophic level}trophic level for the top predator, $P$ (Fig. \ref{fig:web}), and let us assume that $P$ gets two-thirds of what it needs from $B$, and gets one-third from $A$. $B$ itself is on the second trophic level, so given that, the trophic level of $P$ is calculated as \begin{equation*} T_{i}=1+\sum_{j=1}^{S}T_{j}p_{ij}=1+\left(2\left(0.67\right)+1\left(0.33\right)\right)=2.67 \end{equation*} where $T_{i}$ is the trophic level of species $i$, $T_{j}$ is the trophic level of prey species $j$, and $p_{ij}$ is the proportion of the diet of predator $i$ consisting of prey $j$. \item [\textbf{\index{omnivory}Omnivory}] Feeding on more than one \emph{trophic} level ($v>0$, Fig. \ref{fig:web}); it is \emph{not} merely feeding on different species or resources. \item [\textbf{\index{intraguild predation}Intraguild predation}] A type of omnivory in which predation occurs between consumers that share a resource; in Fig. \ref{fig:web} $P$ and $B$ share prey $R$. When $P$ gets most of its energy from $B$, we typically refer to that as omnivory ($v<0.5$); when $P$ gets most of its energy from $R$, we typically refer to that as intraguild predation ($v>0.5$). \end{description} This list of food web descriptors is a fine start but by no means exhaustive. \section{Food chain length --- an emergent property} There are many interesting questions about food webs that we could address; let us address one that has a long history, and as yet, no complete answer: What determines the length of a \index{food chain length}food chain? Some have argued that chance plays a large role \cite{Cohen1990d, Hubbell2001, Williams2000a}, and others have shown that area \cite{MacArthur1967, Rosenzweig1995} or ecosystem size \cite{Post:2002bs} may play roles. The explanation that we will focus on here is \emph{dynamical stability}. Communities with more species had been hypothesized to be less stable, and therefore less likely to persist and be observed in nature. Stuart Pimm and John Lawton \cite{Pimm1977} extended this work by testing whether \emph{food chain length} could be limited by the instability of long food chains \cite{Post:2002bs}. \subsection{Multi-species Lotka--Volterra notation} A many-species \index{Lotka--Volterra!food web}Lotka--Volterra model can be represented in a very compact form, \begin{equation} \label{eq:compact} \frac{\D X_{i}}{\D t}=X_{i}\left(b_{i}+\sum_{j=1}^{S}a_{ij}X_{j}\right) \end{equation} where $S$ is the number of species in the web, $b_i$ is the intrinsic rate of increase of species $i$ (i.e. $r_i$), and $a_{ij}$ is a per capita effect of species $j$ on species $i$. When $i=j$, $a_{ij}$ refers to an \emph{intra}specific effect, which is typically negative. Recall that in our earlier chapters on competition, we used $\alpha_{ii}$ to represent intraspecific per capita effects. Here for notational convenience, we leave $i$ and $j$ in the equation, realizing that $i=j$ for intraspecific interactions. Further, we let $a_{ij}$ be any sign, either positive or negative, and sum the effects. If we let $X=N$, $b=r$, and $a=r \alpha$, then the following are equivalent: \begin{gather*} \dot{N_1} =r_1N_1\left(1-\alpha_{11}N_1 - \alpha_{12}N_2\right)\\ \dot{X} =X_1\left(b_1+a_{11}X_1+a_{12}X_2 \right) \end{gather*} The notation in eq. \ref{eq:compact} is at once both simple and flexible. When $a_{ij}$ is negative, it may represent competition or the effect of a predator, $j$, on its prey, $i$. When $a_{ij}$ is positive, it may represent mutualism or the effect of prey~$j$ on a predator~$i$. \subsection{Background} In the early and mid-1970's, Robert May and others demonstrated that important predictions could be made with relatively simple Lotka--Volterra models \cite{May1973ab}, and this work still comprises an important compendium of lessons for ecologists today \cite{May:2006lr}. May used simple Lotka--Volterra models to show that increasing the number of species in a food web tended to make the food web less stable \cite{May1973ab, May1972}. In species-rich webs, species were more likely to become extinct. Further, he showed that the more connections there were among species in the web (higher connectance), and the stronger those connections (higher interaction strength), the \emph{less} stable the web. At the time, this ran counter to a prevailing sentiment that more diverse ecosystems were more stable, and led to heated discussion. \begin{figure}[ht] \centering \subfloat[]{\includegraphics[width=1.55cm]{PimmandLawtonA.pdf}\label{A}}\qquad \subfloat[]{\includegraphics[width=5cm]{PimmandLawtonE.pdf}\label{E}}\qquad \subfloat[]{\includegraphics[width=3cm]{PimmandLawtonB.pdf}\label{B}} \caption[Three of the six webs investigated by Pimm and Lawton (1977).]{\subref{A}, \subref{E}, and \subref{B} correspond to Pimm and Lawton (1977) Figs. 1A, E, and B. Note that \subref{E} differs from \subref{A} in that \subref{E} has only two trophic levels instead of four. Note also that \subref{B} differs from \subref{A} only in that species 4 has an additional omnivorous link. All basal species exhibit negative density dependence.} \label{fig:PLABE} \end{figure} May used specific quantitative definitions of all of his ideas. He defined connectance as the proportion of interactions in a web, given the total number of all possible directed interactions (i.e., directed connectance). Thus a linear food chain with four species (Fig. \ref{A}), and intraspecific competition in the basal (bottom) species would have a connectance of $4/16=0.25$. May's definition of \index{interaction strength, quantified}interaction strength was the square root of the average of all $a_{ij}^2$ ($i \neq j$), \begin{equation} \label{eq:IS} I =\sqrt{ \frac{\sum_{i=1}^S\sum_{j=1, i \neq j}^S a_{ij}}{S^2-S}}. \end{equation} Squaring the $a_{ij}$ focuses on magnitudes, putting negative and positive values on equal footing. An important component of May's work explored the properties of {\em randomly} connected food webs. At first glance this might seem ludicrous, but upon consideration, we might wonder where else one could start. Often, simpler (in this case, random) might be better. The conclusions from the \index{random connection models}random connection models act as null hypotheses for how food webs might be structured; deviations from May's conclusions might be explained by deviations from his assumptions. Since this work, many ecologists have studied the particular ways in which webs in nature appear non-random. One conclusion May derived was a threshold between stability and instability for random webs, defined by the relation \begin{equation} \label{eq:1} I\left( S C_D \right)^{1/2} = 1 \end{equation} where $I$ is the average \index{interaction strength, average}interaction strength, $S$ is the number of species, and $C_D$ is directed connectance. If $I\left( S C \right)^{1/2} > 1$, the system tended to be unstable (Fig. \ref{fig:may1}). Thus, if we increase the number of species, we need to decrease the average interaction strength if we want them to persist. The larger and more tightly connected (larger $I$, $S$, and $C_D$) the web, the more likely it was to come crashing down. Therefore, if longer food chains were longer by virtue of having more species, they would be less stable because of the extra species, if for no other reason. <>= S <- 40; C=.3 curve(C*x^ -2, .1, 1, ylab=expression("Number of Species ("*italic(S)*")"), xlab = expression("Interaction Strength ("*italic(I)*")") ) text(.6, 20, "Unstable Webs", srt=45) text(.2, 3.75, "Stable Webs" , srt=-45, cex=.8) @ \begin{figure}[ht] \centering \includegraphics[width=.5\textwidth]{MayCriteria1} \caption{\small Relation between the average interaction strength and the number of species able to coexist (here directed connectance is $C_D = 0.3$). The line represents the maximum number of species that are predicted to be able to coexist at equilibrium. Fewer species could coexist, but, on average, more species cannot coexist at equilibrium.} \label{fig:may1} \end{figure} @ Pimm and Lawton felt that it seemed reasonable that long chains might be less stable also because predator-prey dynamics appear inherently unstable, and a connected series of unstable relations seemed less likely to persist than shorter chains. They tested whether food chain length \emph{per se}, and not the number of species, influenced the stability of food chains. Another question they addressed concerned omnivory. At the time, surveys of naturally occurring food webs had indicated that omnivory was rare \cite{Pimm1978}. Pimm and Lawton tested whether omnivory stabilized or destabilized food webs \cite{Pimm1977,Pimm1978}. Like May, Pimm and Lawton \cite{Pimm1977} used Lotka--Volterra models to investigate their ideas. They designed six different food web configurations that varied food chain length, but held the number of species constant (Fig. \ref{fig:PLABE}). For each food web configuration, they varied randomly interaction strength and tested whether an otherwise randomly structured food web was stable. Their food webs included negative density dependence only in the basal species. Pimm and Lawton concluded that (i) shorter chains were more stable than longer chains, and (ii) omnivory destabilized food webs (Fig. \ref{fig:RTHist}). While these conclusions have stood the test of time, Pimm and Lawton failed to highlight another aspect of their data --- that omnivory shortened return times for those webs that were qualitatively stable (Fig. \ref{fig:histsAB}). Thus, omnivory could make more stable those webs that it didn't destroy. Subsequent work has elaborated on this, showing that weak omnivory is very likely to stabilize food webs \cite{McCann1998, McCann1997}. \section{Implementing Pimm and Lawton's Methods} Here we use \R~code to illustrate how one might replicate, and begin to extend, the work of Pimm and Lawton \cite{Pimm1977}. In previous chapters, we began with explicit time derivatives, found partial derivatives and solved them at their equilibria. Rather than do all this, Pimm and Lawton bypassed these steps and went straight to the evaluated Jacobian matrix. They inserted random estimates for the elements of the Jacobian into each non-zero element in the food web matrix. These estimates were constrained within (presumably) reasonable limits, given large less abundant predators and small more abundant prey. Their methods followed this approach. \begin{compactenum} \item Specify a food web interaction matrix.\footnote{In the original publication, webs E and D seem to be represented incorrectly.} \item Include negative density dependence for basal species only. \item Set upper and lower bounds for the effects of predators on prey ($0$ to $-10$) and prey on predators($0$ to $+0.1$); these are the Jacobian elements. \item Generate a large number of random Jacobian matrices and perform linear stability analysis. \item Determine qualitative stability (test $\lambda_1<0$), and return time for each random matrix. Use these to examine key features of the distributions of return times (e.g., average return time). \item Compare the stability and return times among different food web configurations that varied systematically in food chain length and the presence of omnivory, but held the number of species constant. \end{compactenum} It is worth discussing briefly the \index{Jacobian elements}Jacobian elements. May \cite{May1973ab} defined interaction strength as the Jacobian element of a matrix, which represents the \emph{total effect of one individual on the population growth rate of another species}. Think about how you calculate the Jacobian --- as the partial derivative of one species' growth rate with respect to the size of the other population. It is the instantaneous change in the \emph{population} growth rate per \emph{unit change} in the population of another, at the equilibrium. The units chosen for the Jacobian elements thus mean that individual predators have relatively much larger effects on the population growth rates of prey than \emph{vice versa}. Let's build a function that does what Pimm and Lawton did. There are an infinite number of ways to do this, but this will suffice. First, we'll create a matrix that represents \emph{qualitatively} the simplest longest food chain (Fig. \ref{A}) where each species feeds only on one other species and where no prey are fed upon by more than one consumer. Initially, we will use the values of greatest magnitude used by Pimm and Lawton. <>= Aq = matrix(c( -1, -10, 0, 0, 0.1, 0, -10, 0, 0, 0.1, 0, -10, 0, 0, 0.1, 0), nrow=4, byrow=TRUE) @ Note that this matrix indicates a negative effect of the basal species on itself, large negative effects ($-10$) of each consumer on its prey, and small positive effects of each prey on its consumer. For subsequent calculations, it is convenient to to find out from the matrix itself how big the matrix is, that is, how many species, $S$, are in the web. <<>>= S <- nrow(Aq) @ Next, we create a random realization of this matrix by multiplying each element times a unique random number between zero and 1. For this matrix, that requires $4^2$ unique numbers. <<>>= M <- Aq * runif(S^2) @ Next we perform eigenanalysis on it, retaining the eigenvalues. <<>>= eM <- eigen(M)[["values"]] @ Pimm and Lawton tested whether the dominant eigenvalue was greater than 0 (unstable) and if less than zero, they calculated return time. We will simply record the dominant eigenvalue (the maximum of the real parts of the eigenvalues). <<>>= deM <- max( Re(eM) ) @ Given the stabilizing effect of the intraspecific negative density dependence, we will hang on to that as well. <<>>= intraNDD <- sqrt(sum(diag(M)^2)/S) @ Given lessons from May's work \cite{May1972}, we might also want to calculate the average interaction strength, not including the intraspecific interactions. Here we set the diagonal interactions equal to zero, square the remaining elements, find their average, and take the square root. <<>>= diag(M) <- 0 IS <- sqrt( sum(M^2)/(S*(S-1)) ) @ Recall that weak omnivory is supposed to stabilize food webs \cite{McCann1997}. For webs that include omnivory, we will calculate the interaction strength of omnivory in the same way we do for other interactions, as the square root of the average of the squared $a_{ij}$ (eq. \ref{eq:IS}). We can wrap all this up in a function where we specify the $i,\,j$ of one of the directed omnivorous links.\footnote{It does not matter which we specify, either the $ij$ or the $ji$.} <>= pimmlawton <- function(mat, N=1, omni.i=NA, omni.j=NA, omega=NULL){ S <- nrow(mat) if(is.na(omni.i)) { out <- matrix(NA, nrow=N, ncol=4) colnames(out) <- c("DomEig", "Im", "IntraDD", "I") for(n in 1:N) out[n,] <- { M = runif(S^2) * mat ## Do eigenanalysis eigs <- eigen(M)[["values"]] mx <- which.max(Re(eigs)) deM = Re(eigs)[mx] deMi = Im(eigs)[mx] intraNDD <- sqrt(sum(diag(M)^2)/S) diag(M) <- 0 IS = sqrt( sum(M^2)/(S*(S-1)) ) c(deM, deMi, intraNDD, IS) } } else { out <- matrix(NA, nrow=N, ncol=5) colnames(out) <- c("DomEig", "Im", "IntraDD", "I", "I.omni") for(n in 1:N) out[n,] <- { M = runif(S^2) * mat ## Adjust for McCann type omnivory if(!is.null(omega)) {M[omni.i,omni.j] <- omega*M[omni.i+1,omni.j] M[omni.i+1,omni.j] <- (1-omega)*M[omni.i+1,omni.j] M[omni.j,omni.i] <- omega*M[omni.j,omni.i+1] M[omni.j,omni.i+1] <- (1-omega)*M[omni.j,omni.i+1]} ## Do eigenanalysis eigs <- eigen(M)[["values"]] mx <- which.max(Re(eigs)) deM = Re(eigs)[mx] deMi = Im(eigs)[mx] intraNDD <- sqrt(sum(diag(M)^2)/S) diag(M) <- 0 IS = sqrt( sum(M^2)/(S*(S-1)) ) omnivory <- sqrt(mean(c(M[omni.i,omni.j],M[omni.j,omni.i])^2)) c(deM, deMi,intraNDD, IS, omnivory) } } return(as.data.frame(out)) } <<>>= args(pimmlawton) @ Now we can check this function for a single simulation for our first web, <<>>= set.seed(1) pimmlawton(Aq) @ Now let's do it 2000 times, as Pimm and Lawton did. Each row will be an independent randomization, and the columns will be the dominant eigenvalue, the intraspecific density dependence, and the average interaction strength. <<>>= out.A <- pimmlawton(Aq, N=2000) @ We might like to look at basic summary statistics of the information we collected --- what are their minima and maxima and mean? <<>>= summary(out.A) @ We see that out of 2000 random food chains, the largest dominant eigenvalue is still less than zero ($\lambda_1<0$). What does that mean? It means that all of the chains are qualitatively stable, and that the return times are greater than zero ($-1/\lambda_1>0$).\footnote{Recall that a negative return time indicates that any ``perturbation'' at the equilibrium would have been closer to zero at some time in the past, i.e., that the perturbation is growing.} May's work showed that stability is related to interaction strength. let's examine how the dominant eigenvalue is related to interaction strength.\footnote{Recall that if we think of stability analysis as the analysis of a small perturbation at the equilibrium, then the dominant eigenvalue is the growth rate of that perturbation.} <>= pairs(out.A) @ \begin{figure}[ht] \centering \includegraphics[width=.9\linewidth]{pairsA} \caption{Perturbations at the equilibrium tend to dissipate more rapidly (more negative dominant eigenvalues) with greater intraspecific negative density dependence (\textsf{IntraDD}) and greater interspecifiic interaction strength (\textsf{I}). This graph also demonstrates the independence of \textsf{IntraDD} and \textsf{I} in these simulations.} \label{fig:pairsA} \end{figure} The results of our simulation (Fig. \ref{fig:pairsA}) show that the dominant eigenvalue can become more negative with greater intraspecific negative density dependence (\textsf{IntraDD}) and greater intersepcifiic interaction strength (\textsf{I}). Recall what this means --- the dominant eigenvalue is akin to a perturbation growth rate at the equilibrium and is the negative inverse of return time. Therefore, stability can increase and return time decrease with increasing interaction strengths. Note also (Fig. \ref{fig:pairsA}) that many eigenvalues seem very close to zero --- what does this mean for return times? The inverse of a very small number is a very big number, so it appears that many return times will be very large, and therefore effectively zero. Let's calculate return times and examine a summary. <<>>= RT.A <- -1/out.A[["DomEig"]] summary(RT.A) @ We find that the maximum \index{return time}return time is a very large number, and even the median is fairly large (\Sexpr{round( median(RT.A) )}). In an every changing world, is there any meaningful difference between a return time of 1000 generations \emph{vs.} neutral stability? Pimm and Lawton addressed this by picking an arbitrarily large number (150) and recording the percentage of return times greater than that. This percentage will tell us the percentage of webs are not ``effectively'' stable. <<>>= sum( RT.A > 150 ) / 2000 @ Now let's extract the return times that are less than or equal to 150 and make a histogram with the right number of divisions or bins to allow it to look like the one in the original \cite{Pimm1977}. <>= A.fast <-RT.A[RT.A < 150] A.slow <-RT.A[RT.A >= 150] # Opens a new plotting device. On Macs use quartz() or x11() histA <- hist(A.fast, breaks=seq(0,150, by=5), main=NULL) @ This histogram (Fig. \ref{hA}) provides us with a picture of the stability for a food chain like that in Fig. \ref{A}. Next, we will compare this to other webs. @ \begin{figure}[ht] \centering \subfloat[Four level chain]{\includegraphics[width=.32\textwidth]{histA.pdf} \label{hA}} \subfloat[Two level chain]{\includegraphics[width=.32\textwidth]{histE.pdf} \label{hE}} \subfloat[Four levels + omnivory]{\includegraphics[width=.32\textwidth]{histB.pdf} \label{hB}} \caption{Histograms for three of the six food chains (A, E, and B) used by Pimm and Lawton.} \label{fig:RTHist} \end{figure} \section{Shortening the Chain} Now let's repeat all this (more quickly) for a shorter chain, but with the same number of species (Fig. \ref{E}). So, we first make the web function. <<>>= Eq = matrix(c( -1, 0, 0, -10, 0, -1, 0, -10, 0, 0, -1, -10, 0.1, 0.1, 0.1, 0), nrow=4, byrow=TRUE) @ Next we run the 2000 simulations, and check a quick summary. <<>>= out.E <- pimmlawton(Eq, N=2000) summary(out.E) @ The summary shows that, again, that all webs are stable ($\lambda_1<0$). A histogram of return times also shows very short return times (Fig. \label{hE}). Plots of $\lambda_1$ \emph{vs.} interaction strengths show that with this short chain, and three basal species that the role of intraspecfic density dependence becomes even more important, and the predator-prey interactions less important in governing $\lambda_1$. <>= layout(matrix(1:2,nr=1)) plot(DomEig ~ IntraDD, data=out.E) plot(DomEig ~ I, data=out.E) @ \begin{figure}[ht] \centering \includegraphics[width=.95\linewidth]{E1} \caption{For a food chain with two levels, and three basal species, perturbation growth rate ($\lambda_1$) declines with increasing intraspecific negative density dependence (\textsf{IntraDD}) and is unrelated to predator-prey interaction strengths.} \label{fig:E1} \end{figure} Note that with the shorter food chain, a greater proportion of the $\lambda_1$ are more negative (farther away from zero) than in the four level food chain. Clearly then, shortening the web stabilizes it, in spite of still having the \emph{number} of species. Let us again categorize these as having long and short \index{return time}return times, and graph the distribution of the short ones. <>= RT.E <- -1/out.E[["DomEig"]] E.fast <-RT.E[RT.E < 150] E.slow <-RT.E[RT.E >= 150] # Opens a new plotting device. On Macs use quartz() or x11() histE <- hist(E.fast, breaks=seq(0,150, by=5), main=NULL) @ \section{Adding Omnivory} Real webs also have \index{omnivory}omnivory --- feeding on more than one trophic level. A nagging question, then and now, concerns the effect of omnivory on food web dynamics. Pimm and Lawton compared food web dynamics with and without omnivory. Let's now create the web (Fig. \ref{B}) that they used to compare directly with their linear food chain (Fig. \ref{A}). <<>>= Bq = matrix(c( -1, -10, 0, 0, 0.1, 0, -10, -10, 0, 0.1, 0, -10, 0, 0.1, 0.1, 0), nrow=4, byrow=TRUE) @ Next we run the 2000 simulations, and check a quick summary. <<>>= out.B <- pimmlawton(Bq, N=2000, omni.i=2, omni.j=4) summary(out.B) @ With omnivory, we now see that most webs have $\lambda_1>0$, and thus are \emph{unstable}. This was one of the main points made by Pimm and Lawton. Let's look at the data. <>= pairs(out.B) @ It means that most of the randomly constructed webs were not stable point equilibria. To be complete, let's graph what Pimm and Lawton did. <>= RT.B <- -1/out.B[["DomEig"]] B.fast <- RT.B[RT.B < 150 & RT.B>0 ] out.B.fast <- out.B[RT.B < 150 & RT.B>0 ,] out.B.stab <- out.B[RT.B>0 ,] histB <- hist(B.fast, breaks=seq(0,150, by=5), main=NULL) @ \subsection{Comparing Chain A versus B} Now let's compare the properties of the two chains, without, and with, omnivory, chains \textbf{A} and \textbf{B} (Figs. \ref{A}, \ref{B}). Because these are stochastic simulations, it means we have \emph{distributions} of results. For instance, we have a distribution of return times for chain \textbf{A} and a distribution for return times for chain \textbf{B}. That is, we can plot histograms for each of them. Pimm and Lawton compared their webs in common sense ways. They compared simple summaries, including \begin{itemize} \item the proportion of random webs that were stable (positive return times), \item the proportion of stable random webs with return times greater than 150. \end{itemize} Now let's try graphical displays. Rather than simply chopping off the long return times, we use base 10 logarithms of return times because the distributions are so right-skewed. We create a histogram\footnote{Note that now we use probabilities for the $y$-axis, rather than counts. The probability associated with any particular return time is the product of the height of the column and the width of the column (or bin).} of the return times for chain \textbf{A}, and nonparametric density functions for both chain \textbf{A} and \textbf{B}.\footnote{These density smoothers do a good job describing empirical distributions of continuous data, often better than histograms, which have to create discrete categories or ``bins'' for continuous data.} <>= hist(log(RT.A,10), probability=T, ylim=c(0,1), main=NULL, xlab=expression(log[10]("Return Time")) ) lines(density(log(RT.A,10)) ) lines(density(log(RT.B[RT.B>0],10)), lty=2, lwd=2 ) legend("topright", c("Chain A", "Chain B"), lty=1:2, lwd=1:2, bty='n') @ \begin{figure}[ht] \centering \includegraphics[width=.5\textwidth]{histsAB} \caption{Comparing the distributions of return times for chain \textbf{A} and \textbf{B}. "Density" is probability density. The distribution of return times for chain \textbf{A} is the solid black line, and the distribution of return times for chain \textbf{B} is the dashed red line.} \label{fig:histsAB} \end{figure} By overlaying the density function of web B on top of web A return times (Fig. \ref{fig:histsAB}), we make an interesting observation. The omnivorous webs with positive return times (those plotted) actually tended to have shorter return times than the linear chain. Pimm and Lawton noted this, but did not emphasize it. Rather, they sensibly focused on the more obvious result, that over 90\% of the omnivorous webs had negative return times, indicating an absence of a stable point equilibrium. \section{Re-evaluating Take-Home Messages} The primary messages made by Pimm and Lawton\cite{Pimm1977} were that \begin{itemize} \item shorter webs are more stable than long chains, \item omnivory destabilized webs. \end{itemize} These conclusions were a major part of the lively debate surrounding these issues. It was consistent with the apparent observation of the time, that empirical food webs revealed little omnivory \cite{Pimm1977,Pimm1978}, and that food chains in nature seemed to be much shorter than could occur, if \index{primary productivity}primary productivity (green plants and algae) was channeled upwards into a linear chain. Let's consider their assumptions. First, Pimm and Lawton made the argument, as many others have (including us), that systems with stable point equilibria are more likely to persist than systems with oscillations, such as stable limit cycles. That is, we presume a strong correlation between the tendency to oscillate, and the tendency for species to become extinct (i.e., the system to collapse). It is easy to show that a system can be pushed from a stable equilibrium into oscillations which eventually become so big as to drive an element of the system to extinction. This is a very reasonable assumption, but one which is not challenged enough. Other measures of system stability could be used, such as the minimum that occurs in a long series of fluctuations \cite{Huxel1998,McCann1997}. Second, Pimm and Lawton ignored the effects of self-regulated basal species. By exhibiting negative density dependence, the basal species stabilized the web. When Pimm and Lawton made a shorter web, they also added more self-regulated populations. Thus, they necessarily confounded chain length with the number of species with negative density dependence. Which change caused the observed differences among webs? We don't know. Third, they assumed that the effect of web topology (i.e., short \emph{vs.} long chain) was best evaluated with the \emph{average} properties of the topology, rather than the maximum properties of the topology. By these criteria, webs without omnivory were clearly better. On average, webs without omnivory were more often stable than chains with omnivory, even if some of their return times tended to be quite long. Therefore, one might argue that if a web assembles in nature, it is more likely to persist (i.e., be stable) if it lacks omnivory. However, let us consider this preceding argument further. The world is a messy place, with constant insults and disturbances, and resources and environmental conditions fluctuating constantly. In addition, there is a constant rain of propagules dropping into communities, and species abundances are changing all the time. In a sense then, communities are being constantly perturbed. The only webs that can persist in the face of this onslaught are the \emph{most} stable ones, that is the ones with the shortest return times. We just showed that Pimm and Lawton's own analyses showed that \emph{the most stable webs tended to be those with \index{omnivory, stabilizing}omnivory}. Subsequent work supports this claim that omnivory is rampant in nature \cite{Polis1991}, and this is supported by theory that shows weak interactions, including omnivory, stabilize food webs \cite{McCann1998, McCann1997}. Pimm and Lawton made very important contributions to this lengthy debate, and we are finally figuring out how to interpret their results. \section{Summary} Over 35 years ago, May started using simple, highly artificial dynamical descriptions of communities, using mathematical approaches that had been well established, if somewhat controversial, more than 50 years earlier by Alfred Lotka, Vito Volterra, and others \cite{Kingsland:1985kx}. Such simple abstractions are still useful today \cite{Bolker2003a}. May's results, and those of Pimm and Lawton remain logical deductions that have resonance throughout community ecology. May, and Pimm and Lawton showed that under very simple assumptions, adding complexity usually destabilizes food webs. We have found that in practice, it is very difficult to build or restore structurally complex, speciose ecosystems \cite{Bakker:1999qa, Byers:2006bs}. Further, we all now realize that omnivory is quite widespread \cite{Thompson:2007uq}, and additional theory indicates that omnivory can actually stabilize food webs by speeding a return to equilibrium and bounding systems farther from minima \cite{McCann1998, Vandermeer:2006fr}. Simple Lotka--Volterra webs will likely reveal more interesting generalizations in the years ahead. % 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{prob1} % The problem\footnote{Footnote} is described here. The % problem is described here. The problem is described here. % \end{prob} % \begin{prob} % \label{prob1} % \textbf{Problem Heading}\\ % (a) The first part of the problem is described here.\\ % (b) The second part of the problem is described here. % \end{prob} \begin{prob} \textbf{General questions}\\ (a) For each web, write out all four species' differential equations, using row, column subscripts for each parameter. Label species 1 $X_1$ or $N_1$, and the others accordingly. \\ (b) State the type of predator functional response used in these models, explain how you know, and explain the effect that this type of response typically has on dynamics. \\ (c) Describe the role played by intraspecific negative density dependence in these models --- which species have it and what is it likely to do?\\ (d) Explain whether the results of this chapter support the contention that longer food chains are less stable.\\ (e) Explain whether the results of this chapter support the contention that omnivory destabilizes food chains. \end{prob} \begin{prob} \label{prob4} \textbf{More models}\\ (a) Rewrite the above code to replicate the rest of Pimm and Lawton's results.\\ (b) Replicate the results of Pimm and Lawton's companion paper \cite{Pimm1978}.\\ (c) Test the effects of intraspecific negative density dependence. Vary the average \emph{magnitude} of negative density dependence.\\ (d) Design a new simulation experiment of your own.\\ \end{prob}