\SweaveOpts{ eval=TRUE, prefix=FALSE, width=4, height=4, eps=FALSE, include=FALSE}

\chapter{Populations in Space}
\label{cha:metapopulations} % Always give a unique label
% use \chaptermark{}
% to alter or adjust the chapter heading in the running head

% 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

<<CandG,echo=false, results=hide, width=6, height=4.5>>=
setwd("~/Documents/Projects/Book/draft/Chap04")
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))
myQ <- function(w=4,h=4, mr = c(5,4,.5,1.5)) {
  quartz(width=w,height=h)
  par(mar=mr)
}
CandG <- c(32,16,10,9,8,7,8,6,4,5,7,4,2,1,1,1,3,8,12)
NoS <- 1:19
par(mar=c(3,3,1,1), mgp=c(1.5,.5,0))
bp <- barplot(CandG,ylim=c(0,35), 
        ylab="No. of Species", xlab="Number of Sites Occupied")
axis(1, at=bp[seq(1,19, by=2)], labels=seq(1,19, by=2) )

@ 
\begin{figure}[ht]
  \centering
  %\resizebox{12cm}{!}{\includegraphics{Chap04/CollinsGlenn1991fig} }
  \includegraphics[width=.67\linewidth]{CandG.pdf}
  \caption{A frequency distribution of the number of plant species ($y$-axis) that occupy different numbers of  grassland remnants ($x$-axis). Note the \textsf{U}-shaped (bimodal) distribution of the number of sites occupied. Other years were similar \cite{Collins:1991fk}}
\label{fig:CandG}
\end{figure}
Over relatively large spatial scales, it is not unusual to have many species that seem to occur everywhere, and even more species that seem to be found in only one or a few locations. For example, Scott Collins and Susan Glenn \cite{Collins:1991fk} showed that in grasslands, each separated by up to 4\,km, there were more species occupying only one site (Fig. \ref{fig:CandG}, left-most bar) than two or more sites, and also that there are more species occupying all the sites than most intermediate numbers of sites (Fig. \ref{fig:CandG}, right-most bar), resulting in a \textsf{U}-shaped frequency distribution. Illke  Hanski \cite{Hanski1982} coined the rare and common species ``satellite'' and ``core'' species, respectively, and proposed an explanation. Part of the answer seems to come from \emph{the effects of immigration and emigration in a spatial context}. In this chapter we explore mathematical representations of individuals and populations that exist in space, and we investigate the consequences for populations and collections of populations.

\section{Source-sink Dynamics}
In Chapters 1-3, we considered \emph{closed} populations. In contrast, one could imagine a population governed by \emph{b}irths plus \emph{i}mmigration, and \emph{d}eaths plus \emph{e}migration (a \emph{BIDE} model). Ron Pulliam \cite{Pulliam:1988ez} proposed a simple model that includes all four components of \emph{BIDE} which provides a foundation for thinking about connected subpopulations. We refer to the dynamics of these as \emph{source-sink dynamics}. Examples of such linked populations might include many different types of species. For instance, a source-sink model could describe linked populations of a single species might occupy  habitat patches of different quality, where organisms might disperse from patch to patch.
\begin{figure}[ht]
 \centering
\begin{tikzpicture}[>=latex',join=bevel,]
  \pgfsetlinewidth{1bp}
%%
\begin{scope}
  \pgfsetstrokecolor{black}
  \draw (113bp,17bp) node[right,above=3cm] {$\mathrm{Habitat}~1 \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \mathrm{Habitat}~2$};
\end{scope}
  \pgfsetcolor{black}
  % Edge: p1 -> p2
  \draw [->] (70bp,63bp) .. controls (92bp,63bp) and (122bp,63bp)  .. (156bp,63bp);
  \draw (113bp,73bp) node[right,above=-.30cm] {$i_{21} = e_{12}$};
  % Edge: p2 -> p1
  \draw [->] (160bp,48bp) .. controls (153bp,44bp) and (146bp,41bp)  .. (138bp,39bp) .. controls (117bp,34bp) and (110bp,34bp)  .. (88bp,39bp) .. controls (84bp,40bp) and (80bp,42bp)  .. (66bp,48bp);
  \draw (113bp,52bp) node[right,below=.5cm] {$i_{12} = e_{21}$};
  % Node: p2
\begin{scope}
  \pgfsetstrokecolor{black}
  \draw (183bp,63bp) ellipse (27bp and 27bp);
  \draw (183bp,63bp) node {$\lambda_2 < 1$};
\end{scope}
  % Node: p1
\begin{scope}
  \pgfsetstrokecolor{black}
  \draw (43bp,63bp) ellipse (27bp and 27bp);
  \draw (43bp,63bp) node {$\lambda_1 > 1$};
\end{scope}
%
\end{tikzpicture}
% End of code
\caption{The simplest source-sink model.}
\end{figure}


\paragraph{The concept}
The general idea of \index{source-sink}source-sink populations begins with the idea that spatially separated subpopulations occupy distinct patches, and each exhibit their own intrinisic dynamics due to births and deaths; that is, we could characterize a \index{lambda!source-sink}$\lambda$ for each subpopulation. In addition, individuals move from one patch to another; that is, they immigrate and emigrate from one patch (or subpopulation) to another. Therefore, the number of individuals we observe in a particular patch is due, not only to the $\lambda$ in that population, but also to the amount of immigration, $i$, and emigration, $e$.\index{sink|see{source-sink}} 

Subpopulations with more births  than deaths, $\lambda > 1$, and with more emigration  than immigration, $e>i$, are referred to as \emph{source populations}.  Subpopulations with fewer births than deaths, $\lambda < 1$, and with more immigration than emigration, $i>e$, are referred to as \emph{sink populations}.  

When we think about what might \emph{cause} variation in $\lambda$, we typically refer to the \emph{quality} of patches or habitats. Quality might be inferred from $\lambda$, or it might actually be the subject of investigation and independent of $\lambda$ --- typically we think of high quality habitat as having $\lambda > 1$ and poor quality habitat as having $\lambda < 1$. 

\paragraph{The equations}
Pulliam envisioned two linked bird populations where one could track adult reproduction, and adult and juvenile survival and estimate $\lambda$, per capita growth rate separately for each population. For the first population, the number of birds in patch 1 at time $t+1$, $n_{1,t+1}$, is the result of adult survival $P_A$, reproduction $\beta_1$, and survival of the juveniles $P_J$. Thus,
\begin{equation}
  \label{eq:ss1}
  n_{1, t+1} = P_An_{t} + P_J \beta_1 n_{1,t} = \lambda_1 n_1.
\end{equation}
Here $\beta_1 n_{1,t}$ is production of juveniles, and $P_J$ is the survival of those juveniles to time $t+1$. Pulliam described the second poopulation in the same fashion as 
\begin{equation}
  \label{eq:ss2}
  n_{2, t+1} = P_An_{t} + P_J \beta_2 n_{1,t} = \lambda_2 n_2.
\end{equation}

Pulliam then assumed, for simplicity's sake, that the two populations vary only in fecundity ($\beta$), which created differences in $\lambda_1$ and $\lambda_2$. He called population 1 the \emph{source population} ($\lambda_1>1$) and population 2 the \emph{sink population} ($\lambda_2 <1$). He also assumed that birds in excess of the number of territories in the source population emigrated from the source habitat to the sink habitat. Therefore, the source population held a constant density (all territories filled), but the size of the population in the sink depended on both its own growth rate $\lambda_2 < 1$ and also the number of immigrants.

\paragraph{A result}
One of his main theoretical findings was that \emph{population density can be a misleading indicator of habitat quality} (Fig. \ref{fig:ssP}). If we assume that excess individuals in the source migrate to the sink, then as habitat quality and reproduction increase in the source population, the source population comprises \emph{an ever decreasing proportion} of the total population! That is, as $\lambda_1$ gets larger, $n_1/(n_1+n_2)$ gets smaller. Thus, density can be a very misleading predictor of long-term population viability, if the source population is both productive and exhibits a high degree of emigration.

\paragraph{A model}
We can use a matrix model to investigate source-sink populations \cite{Berry:2008}. Let us mix up typical demographic notation (e.g., Chapter 2) with that of Pulliam \cite{Pulliam:1988ez}, so that we can recognize Pulliam's quantities in a demographic matrix model setting. Further, let us assume a pre-breeding census, in which we count adults. The population dynamics would thus be governed by \textbf{A}

\begin{equation}
\label{SSBerry}
\mathbf{A} = \left(
\begin{array}{cccc}
P_{A1} + P_{J1}\beta_1 &  M_{12}                    \\
M_{21}                &  P_{A2} +P_{J2}\beta_2 \\
\end{array} \right)
\end{equation}

where the upper left element (row 1, column 1) reflects the within-patch growth characteristics for patch 1. The lower right quadrant (row 2, and column 2) reflects the within-patch growth characteristics of patch 2. 

We then assume, for simplicity,  that migration, $M$, is exclusively from the source to the sink ($M_{21}>0$, $M_{12}=0$). We further assume that $\lambda_1 > 1$ but all excess individuals migrate to patch 2, so $M_{21} = \lambda_1 - 1>0$. Then \textbf{A} simplifies to 

\begin{equation}
\label{SSPulliam}
\mathbf{A} = \left(
\begin{array}{cccc}
1 &  0              \\
\lambda_1-1 &  \lambda_2 \\
\end{array} \right)
\end{equation}


\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{The spatial demographic Pulliam-like model}
We first assign $\lambda$ for the source and sink populations, and create a matrix.
<<>>=
L1 <- 2; L2 <- .4
A <- matrix(c(1,    0, 
              L1-1, L2), nrow=2, byrow=TRUE)
@ 
We can then use eigenanalysis, as we did in Chapter 2 for stage structured populations. The dominant eigenvalue will provide the long term asymptotic total population growth. We can calculate the stable ``stage'' distribution, which in this case is the distribution of individuals between the two habitats. 
<<>>=
eigen(A)
@
From the dominant eigenvalue, we see Pulliam's working assumption that the total population growth is set at $\lambda=1$. We also see from the dominant eigenvector that the sink population actually contains more individuals than the source population (0.51/(0.51+0.86) < 0.5). 

We could graph these results as well, for a range of $\lambda_1$. Here we let \texttt{p1} be the proportion of the population in the source. 
<<sourceP, fig=TRUE>>=
L1s <- seq(1, 3, by=0.01)
p1 <- sapply(L1s, function(l1) {A[2,1] <- l1-1;
                                eigen(A)$vectors[1,1]/sum(eigen(A)$vectors[,1])
                                })
plot(L1s, p1, type='l', ylab="Source Population", xlab=expression(lambda[1]) ) 
@
}
\end{boxedminipage}\medskip

\begin{figure}[ht]
  \centering
  \includegraphics[width=.5\linewidth]{sourceP.pdf}  
  \caption{The declining \emph{relative} abundance in the high quality habitat in a source-sink model. The proportion of the total population ($n_1/(n_1+n_2)$) in the \emph{source} population may decline with increasing habitat quality and growth rate $\lambda_1$ habitat.}
\label{fig:ssP} 
\end{figure}


\section{Two Types of \index{metapopulation}Metapopulations}
Our logistic model (Chapter 3) is all well and good, but it has no concept of \emph{space} built into it. In many, and perhaps most circumstances in ecology, space has the potential to influence the dynamics of populations and ecosystem fluxes \cite{Lehman1997,Leibold:2004fe,Loreau:2003hb}. The logistic equation represents a \emph{closed} population, with no clear accounting for emigration or immigration. In particular cases, however, consideration of space may be essential. What will we learn if we start considering space, such that sites are open to receive immigrants and lose emigrants?

First we consider ideas associated with different types of ``collections;'' we then consider a mathematical framework for these ideas.

\paragraph{A single spatially structured population} 
One conceptual framework that we will consider below is that of a single closed population, where individuals occupy sites in an implicitly spatial context (Fig. \ref{fig:metpop}). Consider a population in space, where a \emph{site} is the space occupied by one individual. One example might be grasses and weeds in a field. In such a population, for an individual within our population to successfully reproduce and add progeny to the population, the individual must first actually occupy a site. For progeny to establish, however, a propagule must arrive at a site that is unoccupied. Thus the more sites that are already occupied, the less chance there is that a propagule lands on an unoccupied site. Sites only open up at some constant per capita rate as individuals die at a per capita death rate.

\begin{figure}[ht]
  \centering
  \subfloat[A closed collection]{%
  \includegraphics[width=.38\linewidth]{metapop.pdf}\label{closed}}
  \subfloat[An open collection]{%
  \includegraphics[width=.5\linewidth]{metapop2b.pdf}\label{open}}
  \caption{Collections of sites. \subref{closed} Sites may be recolonized via internal propagule production and dispersal only, or \subref{open} sites  may receive immigrants from an outside source that is not influenced by the collection. Each site (A-F) may be a spot of ground potentially occupied by a single plant, or it may be an oceanic island potentially occupied by a butterfly population. Sites may also be colonized via both internal and external sources.}
  \label{fig:metpop}
\end{figure}

\paragraph{A metapopulation}
The other conceptual framework that we consider here is that of \emph{\index{metapopulation}metapopulations}. A metapopulation is a population of populations, or a collection of populations (Fig. \ref{fig:metpop}). Modeling metapopulations emerged from work in pest management when Levins \cite{Levins1969} wanted to represent the dynamics of the proportion of fields infested by a pest. He assumed that a field was either occupied by the pest, or not. The same models used to represent a population of individuals that occupy sites (above) can also be used to represent populations that occupy sites, with conceptually similar ecological interpretation. In this case, each \emph{site} is a location that either contains a population or not. In this framework, we keep track of the proportion of all populations that remain extant, that is, the proportion of sites that are occupied. As with a single population (above), the metapopulation is closed, in the sense that there exists a finite number of sites which may exchange migrants.

Whether we consider a single spatial population, or single metapopulation, we can envision a \emph{collection of sites} connected by dispersal. Each site may be a small spot of ground that is occupied by a plant, or it may be an oceanic island that is occupied by a population. All we know about a single site is that it is occupied or unoccupied. If the site is occupied by an individual, we know nothing of how big that individual is; if the site is occupied by a population, we know nothing about how many indiviuals are present. The models we derive below keep track of the proportion of sites that are occupied. These are known loosely as \emph{metapopulation models}. Although some details can differ, whether we are modeling a collection of spatially discrete individuals in single population or a collection of spatially discrete populations, these two cases share the idea that there are a collection of sites connected by migration, and each is subject to extinction. 

The most relevant underlying biology concerns colonization and extinction in our collection of sites (Fig. \ref{fig:metpop}). In this chapter, we will assume that all sites experience equal rates; when we make this assumption, we greatly simplify everything, and we can generalize across all sites. All of the models we consider are simple elaborations of what determines colonization and extinction. Another useful concept to consider is whether the collection of sites receives propagules from the outside, from some external source that is not influenced by the collection of sites (Fig. \ref{fig:metpop}).


\section{Related Models}
Here we derive a single mathematical framework to describe our two types of models.  In all cases, we will consider how total rates of colonization, $C$, and extinction, $E$, influence the the rate of change of $p$, the proportion of sites that are occupied,
\begin{equation}
  \label{eq:4}
    \frac{\D{p}}{\D{t}}=C-E.
\end{equation}
We will consider below, in a somewhat orderly fashion, several permutations of how we represent colonization and extinction of sites (e.g., \cite{Gotelli1993, Gotelli1991}). 

\subsection{The classic Levins model}
\index{metapopulation!Levins}Levins \cite{Levins1969}  proposed what has come to be known as the classic metapopulation model,
\begin{equation}
  \label{eq:intcol}
  \frac{\D{p}}{\D{t}} = c_{i}p\left(1-p\right) - ep. 
\end{equation}
This equation describes the dynamics of the proportion, $p$, of a set of fields invaded by a pest (Fig. \ref{fig:levins}). The pest colonizes different fields at a total rate governed by the rate of propagule production, $c_i$, and also on the proportion of patches that contain the pest, $p$. Thus, propagules are being scattered around the landscape at rate $c_ip$. The rate at which $p$ changes, however, is also related to the proportion of fields that are unoccupied, $(1-p)$, and therefore available to become occupied and increase $p$.  Therefore the total rate of colonization is $c_ip(1-p)$. The pest has a constant local extinction rate $e$, so the total extinction rate in the landscape is $ep$.  

The parameters $c_i$ and $e$ are very similar to $r$ of continuous logistic growth, insofar as they are dimensionless instantaneous rates. However, they are sometimes thought of as probabilities. The parameter $c_i$ is approximately the proportion of open sites colonized per unit time. For instance, if we created or found 100 open sites, we could come back in a year and see how many became occupied over that time interval of one year, and that proportion would be a function of $c_i$. The parameter $e$ is often thought of as the probability that a site becomes unoccupied per unit time. If we found 100 occupied sites in one year, we could revisit them a year later and see how many became \emph{un}occupied over that time interval of one year.

We use the subscript $i$ to remind us that the colonization is coming from within the sites that we are studying (i.e. internal colonization). With internal colonization, we are modeling a closed spatial population of sites, whether ``site'' refers to an entire field (as above), or a small patch of ground occupied by an individual plant \cite{Tilman1994}. 

\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{The Levins metapopulation model (Fig. \ref{fig:levins})}
A function for a differential equation requires arguments for time, a vector of the state variables (here we have one state variable, $p$), and a vector of parameters.
<<>>=
levins <- function(t, y, parms) {
       p <- y[1]
       with(as.list(parms), {
         ## Here is the equation
       dp <- ci * p * (1-p) - e * p
       return(list(dp))
       }) }
@ 
By using \texttt{with}, we can specify the parameters by their names, as long as \texttt{parms} includes names. The function returns a list that contains a value for the derivative, evaluated at each time point, for each state variable (here merely $\D p /\D t$). We then use \texttt{levins} in the numerical integration function \texttt{ode} in the \texttt{deSolve} package.
<<>>=
library(deSolve)
prms <- c(ci=0.15, e=0.05)
Initial.p <- 0.01
out.L <- data.frame(ode(y=Initial.p, times=1:100, func=levins, parms=prms ))
@ 
We then plot the result (Fig. \ref{fig:levins}).
<<levins, fig=TRUE>>=
plot(out.L[,2] ~ out.L[,1],   type='l', ylim=c(0,1), ylab="p", xlab="time")  
@
}
\end{boxedminipage}\medskip

Can we use this model to predict the eventual equilibrium? Sure --- we just set eq. \ref{eq:intcol} to zero and solve for $p$. This model achieves and equilibrium at,
\begin{align*}
  \label{eq:3}
  0 &=c_{i}p-c_{i}p^2 - ep\\
  p^* &= \frac{c_{i}-e}{c_{i}}=1-\frac{e}{c_{i}}.
\end{align*}
 When we do this, we see that $p^* > 0$ as long as $c_i > e$ (e.g., Fig. \ref{fig:levins}). When is $p^* = 1$, so that all the sites are filled? In principle, all sites cannot be occupied simultaneously unless $e=0$!

\begin{figure}[ht]
  \centering
  \subfloat[Levins]{\includegraphics[width=.45\linewidth]{levins.pdf}  \label{fig:levins} }
    \subfloat[Propagule Rain and Hanski]{\includegraphics[width=.45\linewidth]{IMH.pdf}  \label{fig:IMH} }
  \caption{Three metapopulation models, using similar parameters ($c_i=0.15$, $c_e=0.15$, $e=0.05$).}

\end{figure}

\subsection{Propagule rain}
From where else might propagules come? \index{metapopulation!propagule rain}\index{metapopulation!Gotelli}If a site is not closed off from the rest of the world, propagules could come from outside the collection of sites that we are actually monitoring. 

For now, let us assume that our collection of sites is continually showered by propagules from an external source. If only those propagules are important, then we could represent the dynamics as, 
\begin{equation}
  \label{eq:proprain}
  \frac{\D{p}}{\D{t}} = c_{e}\left(1-p\right) - ep 
\end{equation}
where  $c_{e}$ specifies rate of colonization coming from the external source. Gotelli \cite{Gotelli1991} refers to this model as a metapopulation model with ``propagule rain'' or the ``island--mainland'' model. He calls it this because it describes a constant influx of propagules which does not depend on the proportion, $p$, of sites occupied for propagule production. Extinction here is mediated only by the proportion of sites occupied, and has a constant per site rate. 

\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{The propagule rain metapopulation model (Fig. \ref{fig:IMH})}
A function for a differential equation requires arguments for time, a vector of the state variables (here we have one state variable, $p$), and a vector of parameters.
<<>>=
gotelli <- function( t, y, parms) {
       p <- y[1]
       with(as.list(parms), {
         ## Here is the equation
       dp <- ce *  (1-p) - e * p
       return(list(dp))
       })
       }

@ 
The function returns a list that contains a value for the derivative, evaluated at each time point, for each state variable (here merely $\D p /\D t$.
}
\end{boxedminipage}\medskip


We can solve for this model's equilibrium by setting eq. \ref{eq:proprain} equal to zero.
\begin{align}
0  &= c_{e}-c_{e}p - ep\\
p^*&=\frac{c_{e}}{c_{e}+e}.
\end{align}

Of course, we might also think that both internal and external sources are important, in which case we might want to include both sources in our model,
\begin{align}
  \label{eq:genmod}
  \frac{\D{p}}{\D{t}} &= (c_{i}p+c_{e})\left(1-p\right) - ep\\
\end{align}
As we have seen before, however, adding more parameters is not something we take lightly. Increasing the number of parameters by, in this case, 50\% could require a lot more effort to estimate.

\subsection{The rescue effect and the core-satellite model}

Thus far, we have ignored what happens between census periods. Imagine that we sample site ``A'' each year on 1 January. It is possible that between 2 January and 31 December the population at site A becomes extinct and then is subsequently recolonized, or ``rescued'' from extinction. When we sample on 1 January in the next year, we have no way of knowing what has happened in the intervening time period. We would not realize that the population had become extinct and recolonization had occurred.

We can, however, model total extinction rate $E$ with this \emph{rescue effect},
\begin{equation}
  \label{eq:reseff}
  E = - ep\left( 1-p \right).
\end{equation}
Note that as $p \rightarrow 1$, the total extinction rate approaches zero.  Total extinction rate declines because as the proportion of sites occupied increases, it becomes increasingly likely that dispersing propagules will land on all sites. When propagules happen to land on sites that are on the verge of extinction, they can ``rescue'' that site from extinction.

Brown and Kodric-Brown \cite{Brown1977a} found that enhanced opportunity for immigration seemed to reduce extinction rates in arthropod communities on thistles. They coined this effect of immigration on extinction as the ``\index{metapoopulation!rescue effect}rescue effect.'' MacArthur and Wilson \cite{MacArthur1963} also discussed this idea in the context of island biogeography. We can even vary the strength of this effect by adding yet another parameter $q$, such that the total extinction rate is $ ep\left( 1-qp \right)$  (see \cite{Gotelli1993}).

Assuming only internal propagule supply and the simple rescue effect results in what is referred to as the the \index{metapopulation!core-satellite}core-satellite model,
\begin{align}
  \label{eq:CoreSat1}
  \frac{\D{p}}{\D{t}} &= c_ip\left(1-p\right) - ep\left( 1-p \right)
\end{align}
This model was made famous by Illka \index{metapopulation!Hanski}Hanski \cite{Hanski1982}. It is referred to as the \emph{core-satellite} model, for reasons we explore later. 

\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{The core-satellite metapopulation model}
A function for a differential equation requires arguments for time, a vector of the state variables (here we have one state variable, $p$), and a vector of parameters.
<<>>=
hanski <- function( t, y, parms) {
       p <- y[1]
       with(as.list(parms), {
         ## Here is the equation
       dp <- ci * p *  (1-p) - e * p *(1-p)
       return(list(dp))
       })
       }
@ 
The function returns a list that contains a value for the derivative, evaluated at each time point, for each state variable (here merely $\D p /\D t$).
}
\end{boxedminipage}\medskip


\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{Graphing propagule rain and core-satellite  models (Fig. \ref{fig:IMH})}
First, we integrate the models using the same parameters as for the Levins model, and collect the results.
<<>>=
prms <- c(ci<- 0.15, ce <- 0.15, e=0.05)
out.IMH <- data.frame(ode(y=Initial.p, times=1:100, func=gotelli, parms=prms ))
out.IMH[["pH"]]<- ode(y=Initial.p, times=1:100, func=hanski, parms=prms )[,2]
@ 
We then plot the rsult (Fig. \ref{fig:levins}).
<<IMH, fig=TRUE>>=
matplot(out.IMH[,1], out.IMH[,2:3], type='l', col=1, ylab="p", xlab="time")  
legend("topleft", c("Hanski", "Propagule Rain"), lty=2:1, bty='n')
@
}
\end{boxedminipage}\medskip


\subsubsection{Core-satellite equilibria}
What is the equilibrium for the Hanski model (eq. \ref{eq:CoreSat1})? We can rearrange this to further simplify solving for $p^*$.
\begin{align}
  \label{eq:CoreSat2}
  \frac{\D{p}}{\D{t}} &= \left(c_i-e \right) p \left(1-p\right)
\end{align}
This shows us that for any value of $p$ between zero and one, the sign of the growth rate (positive or negative) is determined by $c_i$ and $e$. If $c_i>e$, the rate of increase will always be positive, and because occupancy cannot exceed 1.0, the metapopulation will go to full occupancy ($p^*=1$), and stay there. This equilibrium will be a \emph{stable attractor} or stable equilibrium. What happens if for some reason the metapopulation becomes globally extinct, such that $p=0$, even though $c_i>e$? If $p=0$, then like logistic growth, the metapopulation stops changing and cannot increase. However, the slightest perturbation away from $p=0$ will lead to a positive growth rate, and increase toward the stable attractor, $p^*=1$. In this case, we refer to $p^*=0$ as an \emph{unstable equilibrium} and a repellor.

If $c_i<e$, the rate of increase will always be negative, and because occupancy cannot be less than 0, the metapopulation will become extinct ($p^*=0$), and stay there. Thus $p^*=0$ would be a stable equilibrium or attractor. What is predicted to happen if, for some odd reason this population achieved full occupancy, $p=1$, even though $c_i<e$? In that case, $(1-p)=0$, and the rate of change goes to zero, and the population is predicted to stay there, even though extinction is greater than colonization. How weird is that? Is this fatal flaw in the model, or an interesting prediction resulting from a thorough examination of the model? How relevant is it? How could we evaluate how relevant it is? We will discuss this a little more below, when we discuss the effects of habitat destruction.

What happens when  $c_i=e$? In that case, $c_i - e = 0$, and the population stops changing. What is the value of $p$ when it stops changing? It seems as though it could be any value of $p$, because if $c_i-e=0$, the rate change goes to zero. What will happen if the population gets perturbed --- will it return to its previous value? Let's return to question in a bit.

To analyze stability in logistic growth, we examined the slope of the partial derivative at the equilibrium, and we can do that here. We find that the partial derivative of eq. \ref{eq:CoreSat1} with respect to $p$ is
\begin{equation}
  \label{eq:CoreSatPD1}
  \frac{\partial \dot{p}}{\partial p} = c - 2cp - e + 2ep
\end{equation}
where $\dot{p}$ is the time derivative (eq. \ref{eq:CoreSat1}). A little puzzling and rearranging will show
\begin{equation}
  \label{eq:CoreSatPD2}
  \frac{\partial \dot{p}}{\partial p} = \left(c_i - e\right)\left(1-2p\right)
\end{equation}
and make things simpler. Recall our rules with regard to stability (Chapter 3). If the partial derivative (the slope of the time derivative) is negative at an equilibrium, it means the the growth rate approaches zero following a perturbation, meaning that it is stable. If the partial derivative is positive, it means that the change accelerates away from zero  following the perturbation, meaning that the equilibrium is unstable. So, we find the following guidelines: 
\begin{itemize}
\item $c_i>e$
  \begin{itemize}
  \item $p=1$, $\partial \dot{p}/ \partial p < 0$, stable equilibrium.   
  \item $p=0$, $\partial \dot{p}/ \partial p>0$, unstable equilibrium.
  \end{itemize}
 \item $c_i < e$
  \begin{itemize}
  \item $p=1$, $\partial \dot{p}/ \partial p > 0$, unstable equilibrium.   
  \item $p=0$, $\partial \dot{p}/ \partial p < 0$, stable equilibrium.
  \end{itemize}
\end{itemize}

What if $c_i=e$?  In that case, both the time derivative ($\D p / \D t$) and the partial derivative ($\partial \dot{p}/ \partial p$) are zero \emph{for all values of \emph{p}}. Therefore, if the population gets displaced from any arbitrary point, it will remain unchanged, not recover, and will stay displaced.  We call this odd state of affairs a \emph{neutral equilibrium}. We revisit neutral equilibrium when we discuss interspecific competition and predation.


We can also explore the stability of one of these equilibria by plotting the mettapopulation growth rate as a function of $p$ (Fig. \ref{fig:hanski}). When we set $c_i>e$, and examine the slope of that line at $p^*=1$, we see the slope is negative, indicating a stable equilibrium.

\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{An equilibrium for the core-satellite metapopulation model (Fig. \ref{fig:hanski})}
We first create an expression for the growth itself, $\D p / \D t$. We then plot it, while we evaluate it, on the fly.
<<>>=
dpdtCS <- expression((ci-e)*p*(1-p))
ci <- 0.15; e <- 0.05
p <- seq(0,1, length=50)
<<hanski, fig=TRUE>>=
plot(p, eval(dpdtCS), type='l', ylab="dp/dt")
@ 
}
\end{boxedminipage} \medskip

\subsubsection{Levins $vs$. Hanski}
Why would we use Levins' model instead of Hanski's core-satellite model? To explore this possibility, let's see how the Hanski model might change gradually into the Levins model. First we define the Hanski model with an extra parameter, $a$,
\begin{equation}
  \label{eq:7}
  \frac{\D p}{\D t} = c_ip\left(1-p\right) - ep\left(1-ap\right).
\end{equation}
Under Hanski's model, $a=1$ and under Levins' model $a=0$. If we solve for the equilibrium, we see that 
\begin{equation}
  \label{eq:8hanski}
  p^*=\frac{c-e}{c-ae}
\end{equation}
so that we can derive either result for the two models. In the context of logistic growth, where $K=Hp^*$, this result, eq. \ref{eq:8hanski}, implies that for the Hanski model, $K$ fills all available habitat, whereas the Levins model implies that $K$ fills only a fraction of the total available habitat. That fraction results from the dynamic balance between $c_i$ and $e$.
\begin{figure}[ht]
  \centering
  \includegraphics[width=.5\linewidth]{hanski.pdf}
  \caption{Metapopulation growth rate as a function of $p$, in the core-satellite model ($c_i=0.15$, $e=0.05$). When we plot population growth rate for the core-satellite model, for arbitrary parameter values where $c_i>e$, we see that growth rate falls to zero at full occupancy (i.e., at $p^*=1$). We also see that the slope is negative, indicating that this equilibrium is stable.}
  \label{fig:hanski}
\end{figure}


\section{Parallels with Logistic Growth}
\index{metapopulation!parallels~with~logistic~growth}
It may have already occurred to you that the closed spatial population described here sounds a lot like simple logistic growth. A closed contiguous population, spatial or not, reproduces in proportion to its density, and is limited by its own density. Here we will make the connection a little more clear. It turns out that a simple rearrangement of eq. \ref{eq:intcol} will provide the explicit connection between logistic growth and the spatial population model with internal colonization \cite{Roughgarden1998}. 

Imagine for a moment that you are an avid birder following a population of Song Sparrows in Darrtown, OH, USA (Fig. \ref{fig:MeloDDG1}). If Song Sparrows are limited by the number of territories, and males are competing for territories, then you could think about male Song Sparrows as ``filling up'' some proportion, $p$, of the available habitat. You have already described this population with the logistic growth model ($\D N/\D t = rN(1-\alpha N)$). Lately, however, you have been thinking about how territories, spatially arranged in the landscape, may limit this population. You therefore decide that you would like to use Levins' spatially-implicit metapopulation model instead (eq. \ref{eq:intcol}). How will you do it? You do it by \emph{rescaling} logistic growth\index{scaling~(logistic~growth)}.

Let us start by defining our logistic model variables in other terms. First we define $N$ as
\begin{equation*}
  N=pH
\end{equation*}
where $N$ is the number of males defending territories,  $H$ is the total number of possible territories, and $p$ is the proportion of possible territories occupied at any one time. At equilibrium, $N^*=K=p^*H$, so $\alpha=1/(p^*H)$. Recall that for the Levins model, $p^*=(c_i-e)/c_i$, so therefore,
\begin{equation*}
  \alpha = \frac{c_i}{\left(c_i - e\right)H}.
\end{equation*}

We now have $N$, $\alpha$, and $K$ in terms of $p$, $H$, $c_i$ and $e$, so what about $r$? Recall that for logistic growth, the per capita growth rate goes to $r$ as $N \to 0$ (Chapter 3). For the Levins metapopulation model, the per patch growth rate is 
\begin{align}
  \frac{1}{p}\frac{\D p}{\D t}&=c_{i}\left(1-p\right) - e.
\end{align}
As $p \to 0$ this expression simplifies to $c_i-e$, which is equivalent to $r$. Summarizing, then, we have,
\begin{gather}
  \label{eq:5}
  r=c_i-e\\
  N=pH\\
  \alpha = \frac{1}{K} = \frac{1}{p^*H} = \frac{c_i}{H\left(c_i-e\right)}\\
\end{gather}

Substituting into logistic growth ($\dot{N} = rN(1-\alpha N)$), we now have 
\begin{align}
  \label{eq:6}
  \frac{\D (pH)}{\D t} &= \left(c_i-e\right)pH\left(1-\frac{c_i}{H\left(c_i-e\right)}Hp\right)\\
  &= \left(c_i-e\right)pH - \frac{ c_i-e }{c_i-e}c_ip^2H\\
  &= H\left(c_ip\left(1-p\right) - ep \right)
\end{align}
which is almost the Levins model. If we note that $H$ is a constant, we realize that we can divide both sides by $H$, ending up with the Levins model eq. \ref{eq:intcol}. 

\section{Habitat Destruction}
Other researchers have investigated effects of habitat loss on metapopulation dynamics \cite{Kareiva:1995kx, Nee:1992vn, Tilman1994}. Taking inspiration from the work of Lande \cite{Lande:1987uq, Lande:1988ad}, Karieva and Wennergren \cite{Kareiva:1995kx} modeled the effect of \index{habitat destruction}\index{metapopulation!habitat destruction}habitat destruction, $D$, on overall immigration probability.  They incorporated this into Levins' model as
\begin{equation}
  \label{eq:D}
\frac{\D p}{\D t}=c_ip(1-D-p) - ep
\end{equation}
where $D$ is the amount of habitat destroyed, expressed as a fraction of the original total available habitat. 

\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{Habitat destruction model }
To turn eq. \ref{eq:D} into a function we can use with \texttt{ode}, we have,
<<>>=
lande <- function(t, y, parms) {
       p <- y[1]
       with(as.list(parms), {
         ## Here is the equation
       dp <- ci * p * (1-D-p) - e * p
       return(list(dp))
       }) }
@
}
\end{boxedminipage} \medskip

Habitat destruction, $D$, may vary between 0 ($=$ Levins model) to complete habitat loss 1.0; obviously the most interesting results will come for intermediate values of $D$ (Fig. \ref{fig:lande}).

\begin{figure}[ht]
  \centering
  \includegraphics[width=.5\linewidth]{lande.pdf}
  \caption{Metapopulation dynamics, combining the Levins model and habitat destruction ($c_i=0.15$, $e=0.05$).}
  \label{fig:lande}
\end{figure}

\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{Illustrating the effects of habitat destruction (Fig. \ref{fig:lande})}
We can plot the dynamics for three levels of destruction, including none. We first set all the parameters, and time.
<<>>=
library(deSolve)
prmsD <- c(ci=0.15, e=0.05, D=0)
Ds <- c(0,.2,.5)
Initial.p <- 0.01
t <- 1:200
@ 
We then create an empty matrix of the right size to hold our results, and then integate the ODE.
<<>>=
ps <- sapply(Ds, function(d) {
  prmsD["D"] <- d
  ode(y=Initial.p, times=t, func=lande, parms=prmsD )[,2]
})
@ 
Last, we plot it and add some useful labels.
<<lande, fig=TRUE, width=5, height=5>>=
matplot(t, ps, type='l', ylab="p", xlab="time")
text(c(200,200,200), ps[200,], paste("D = ", Ds, sep=""),
     adj=c(1,0))
@
}
\end{boxedminipage} \medskip

What is the equilibrium under this model? Setting eq. \ref{eq:D} to zero, we can then solve for $p$.
\begin{gather}
  \label{eq:Deq}
  0=c_i-c_iD-c_ip - e\\
  p^* = \frac{c_i-c_iD - e}{c_i} = 1 -  \frac{e}{c_i} - D
\end{gather}
Thus we see that habitat destruction has a simple direct effect on the metapopulation.

\paragraph{A core-satellite habitat loss scenario}

Let us return now to that odd, but logical, possibility in the core-satellite model where $c_i < e$ and $p=1$. Recall that in this case, $p=1$ is an \emph{unstable} equilibrium ($p=0$ is the stable equilibrium for $c_i<e$). We discuss this in part for greater ecological understanding, but also to illustrate why theory is sometimes useful --- because it helps us explore the logical consequences of our assumptions, even when, at first blush, it seems to make little sense.

Imagine that at one time, a metapopulation is regulated by the mechanisms in the core-satellite model, including the rescue effect, and $c_i>e$. We therefore pretend that, the metapopulation occupies virtually every habitable site (let $p=0.999$). Now imagine that the environment changes, causing $c_i<e$. Perhaps human urbanization reduces colonization rates, or climate change enhances extinction rates. All of a sudden, our metapopulation is poised on an unstable equilibrium. What will happen and how might it differ with and without the rescue effect?

When $c_i>e$, we see that $p^*=1$ is the stable attractor (Fig. \ref{fig:lande2}). However, when $c_i < e$, we see the inevitable march toward extinction predicted by the Hanski model (core-satellite) (Fig. \ref{fig:lande2}). Last, when we compare it to the Levins model, we realize something somewhat more interesting. While the Levins model predicts very rapid decline, the Hanski model predicts a much more gradual decline \emph{toward extinction}. Both models predict extinction, but the rescue effect delays the appearance of that extinction. It appears that the rescue effect (which is the difference between the two models) may act a bit like the ``extinction debt'' \cite{Tilman1994} wherein deterministic extinction is merely delayed, but not postponed indefinitely. Perhaps populations influenced by the rescue effect might be prone to unexpected collapse, if the only stable equilibria are 1 and 0. Thus simple theory can provide interesting insight, resulting in very different predictions for superficial similar processes.

\medskip \noindent
\begin{boxedminipage}{\linewidth}
  {\footnotesize
\paragraph{The unexpected collapse of core populations}
Here we plot the dynamics of metapopulations starting at or near equilbrium. The first two use the Hanski model, while the third uses Levins. The second and third use $c_i < e$.
<<>>=
C1 <- ode(y=0.999, times=t, func=hanski, parms=c(ci=0.2, e=0.01))
C2 <- ode(y=0.999, times=t, func=hanski, parms=c(ci=0.2, e=0.25))
L2 <- ode(y=0.95, times=t, func=levins, parms=c(ci=0.2, e=0.25))
@ 
Next, we plot these and add a legend.
<<lande2, fig=TRUE>>=
matplot(t, cbind(C1[,2],C2[,2],L2[,2]), type='l', ylab="p", xlab="Time", col=1)
legend('right', c("c > e", "c < e", "c < e (Levins)"), lty = 1:3, bty='n' )
@
}
\end{boxedminipage} \medskip

\begin{figure}[ht]
  \centering
  \includegraphics[width=.5\linewidth]{lande2}
  \caption{Metapopulation dynamics, starting from near equilibrium for $c_i=0.20$ and $e=0.01$. If the environment changes, causing extinction rate to increase until it is greater than colonization rate, we may observe greatly delayed, but inevitable, extinction (e.g., $c_i=0.20, e=0.25$).}
\label{fig:lande2}
\end{figure}


\section{Core-Satellite Simulations}
\index{metapopulation!core-satellite~simulations}
 Here\footnote{This section relies extensively on code} we explore a simple question that Hanski posed long ago: what would communities look like if all of the populations in the community could be described by their independent own core-satellite model? To answer this question, he created communities as collections of independent (non-interacting) populations that behave according to his metapopulation model with internal colonization and the rescue effect \cite{Hanski1982}. He found that such simulated communities predict that many species will be in almost all sites (``core species''), and even more species will exist at very few sites (``satellite species''). This seems to be a relatively common phenomenon \cite{Collins:1991fk}, and an observation we described at the beginning of the chapter (Fig. \ref{fig:CandG}).

Hanksi's goal was to simulate simultaneously a substantive number of species, to create a community. Each species is assumed to be governed by internal propagule production only, and the rescue effect. Further, he assumed that the long term average density independent growth rate ($r=c_i-e$) was zero. That is, the populations were not \emph{systematically} increasing or decreasing. However, he allowed for stochastic year-to-year variation in probabilities $c_i$ and $e$.  

In these simulations here, we will select the mean for each parameter, $c_i$ and $e$, and the proportion, $\phi$ (``phi'') by which they are allowed to vary. The realized values of $c_{i,t}$ and $e_t$ at any one point in time are random draws from a uniform distribution within the ranges $i\pm \phi i$ and $e\pm \phi e$. (This requires that we do numerical integration at each integer time step since there is no obvious analytical solution to an equation in which the parameters vary through time. This will keep these parameters constant for an entire year, and yet also allow years to vary.)

We start by using the \texttt{args()} function to find out what arguments (i.e. options) are available in the simulation function, \texttt{MetaSim}.
<<>>=
args(MetaSim)
@ 
What options (or arguments) can you vary in \texttt{MetSim}? The `method' may equal \texttt{CoreSatellite}, \texttt{Levins}, \texttt{IslandMainland}, or \texttt{HabitatDestruction}. The default is \texttt{CoreSatellite}; if an argument has a value to begin with (e.g. \texttt{method='CoreSatellite'}), then it will use that value unless you replace it.

Let's start with an initial run of 10 simulations (produces dynamics for 10 populations) to reproduce Hanski's core-satellite pattern by using the rescue effect with equal $i$ and $e$.
<<CoreSatSim10, fig=TRUE>>=
out.CS.10 <- MetaSim(method="hanski", NSims=10)
matplot(out.CS.10$t, out.CS.10$Ns, type="l", xlab="Time", ylab="Occupancy", 
        sub=out.CS.10$method)

@ 
\begin{figure}[ht]
  \centering
  \includegraphics[width=.67\linewidth]{CoreSatSim10}
  \caption{Core-satellite species dynamics with stochasticity ($\bar{i}=\bar{e}=0.2$).}
\label{fig:cssim10}
\end{figure}
These dynamics (Fig. \ref{fig:cssim10}) appear to be completely random. A \emph{random walk} is a dynamic that is a random increase or decrease at each time step. Such a process is not entirely random because the abundance at time $t$ is related to the abundance at time $t-1$, so observations in random walks are correlated through time; they are temporally autocorrelated.

Does a single metapopulation growth rate appear related to $p$, the metapopulation size? What would a deterministic dynamic look like if $c_i>e$? It would increase rapidly at first, and then slow down as it approached 1.0. Can you detect that slow-down here? Similarly, as a metapopulation declines toward extinction, its progression toward $p=0$ slows down. As a result, we tend to accumulate a lot of common and rare species, for which $p$ is close to one or zero.

Now we will do more simulations (50 species), and run them for longer (500 time intervals $vs.$ 50). Doing many more simulations will take a little longer, so be patient\footnote{\texttt{system.time} merely times the process, in secs.}.
@ 
<<>>=
system.time(
out.CS.Lots <- MetaSim(method="hanski",NSims=50, Time=1000)
)
@ Now we can plot the results in two different ways. First, you could do the same
time series, although this may not tell you much. Alternatively, we can plot a histogram of the 50 species' final abundances, at
$t=500$.  
<<CoreSatSimHist, fig=TRUE, width=5, height=5>>=
hist(out.CS.Lots$Ns[501,], breaks=10, main=NULL,
     xlab=expression("Occupancy ("*italic("p")*")"), ylab="Number of Species",sub=paste(out.CS.Lots$method, " Model", sep=""))

@ 
\begin{figure}[ht]
  \centering
  \subfloat[With Rescue Effect]{%
  \includegraphics[width=.48\linewidth]{CoreSatSimHist.pdf}\label{cs}}
\subfloat[No rescue effect]{%
  \includegraphics[width=.48\linewidth]{LevinsSimHist.pdf}\label{lev}}
  \caption{The species-abundance distribution resulting from dynamics for 50 independent metapopulations with internal colonization. \subref{cs} includes the rescue effect (Hanski's model), and  note that most species are either common ($p>0.8$) or rare ($p<0.2$). Levins model \subref{lev} does not include the rescue effect, and there are very few core species ($p>0.8$).}
\label{fig:cssimhist}
\end{figure}
Our simulations (Fig. \ref{fig:cssimhist}) should be consistent with
the core-satellite hypothesis --- are they? In Hanski's model, we see that most metapopulations are either core species ($p>0.8$) or satellite species ($p<0.2$) (Fig. \ref{cs}). This is not to imply that there should be hard rules about what constitutes a core and satellite species, but rather merely shows we have a plethora of both common and uncommon species.

What does the Levins model predict? Let's run the simulations and find out.
<<>>=
system.time(
out.L.Lots <- MetaSim(NSims=50, Time=500, method="levins")
)
@ 
Now we plot a histogram of the 50 species' final abundances, at
$t=500$.
<<LevinsSimHist, fig=TRUE, width=5, height=5>>=
hist(out.L.Lots$Ns[501,], breaks=10,
     xlab=expression("Occupancy ("*italic("p")*")"), 
     ylab="Number of Species", 
     main=NULL, sub=paste(out.L.Lots$method, " Model", sep=""))
@ 
In contrast to the core-satellite model, the Levins model predicts that many fewer species are common (Fig. \ref{lev}). Thus these two population models make contrasting predictions regarding the structure of communities (i.e. relative species abundances), and provide testable alternatives \cite{Collins:1991fk}.
\index{core-satellite|see{metapopulation}}\index{propagule~rain|see{metapopulation}}\index{Levins|see{metapopulation}}

\section{Summary}
In this chapter, we have introduced \emph{space} as an important component of population dynamics. We provided a source-sink framework for linked populations, where population size depends on both intrinsic capacities of a habitat patch, and on immigration and emigration rates. We described (i) a population of individuals within a site, and (ii) a population of populations within a region in the ``metapopulation'' framework. We showed similarities and differences between related metapoopulation models, and between related metapopulation and logistic models. We investigated the response of metapopulations to habitat destruction. Last, we have shown how different population dynamics lead to different community structure.


\section*{Problems}
\addcontentsline{toc}{section}{Problems}
\begin{prob}
  \textbf{Equilibria}\\
  Derive expressions and calculate equilibria for the following metapopulation models, with $c_i=0.05,~e=0.01$. Show your work --- start with the differential equations, set to zero, and solve $p^*$; then substitute in values for $c_i,~e$.\\
(a) Levins model. \\
(b) Propagule rain model (\texttt{gotelli}).\\
(c) Propagule rain model that also includes both external and internal propagule production and dispersal.\\
(d) Hanski model.\\
(e) Lande (habitat destruction) model (with \texttt{D=0.1}).
\end{prob}

\begin{prob}
  \textbf{Habitat destruction}\\
  Compare different levels of habitat destruction.\\
(a) Use the habitat destruction model (\texttt{lande}) to compare 9 levels of destruction (\texttt{ds <- seq(0,.8, by=.1)}), using  $c_i=0.1,~e=0.01$. Plot of graph of the dynamics through time, and calculate the equilibria directly.\\
(b) Write an ODE function for a habitat destruction model with rescue effect. Let the ``rescue'' have an additional parameter, $a$, such that extinction rate is $ep(1-ap)$. \\
(c) Let $D=0.5,~c_i=0.1,~e=0.02$, and vary $a$ over five levels (including $a=0,~1$) to investigate the effects of ``relative rescue effect'' on the equilibria and dynamics of a metapopulation.
\end{prob}