\chapter{Modeling With Systems of Difference Equations}

\specialsection{Introduction}
In our discussion of the assumptions behind logistic models for population dynamics in Chapter 7, we mentioned
\index{logistic growth model}
\index{model!logistic growth}
that the assumption that all individuals in the population are essentially treated as interchangeable is quite unrealistic
 from the biological viewpoint.  Animal species (for which sexual reproduction is the most common strategy) have populations made up 
 of males and females.  Logistic models do not take that into account.  Moreover, 
 recall that we said in Assumptions~\ref{LogAss} that the population was
``unstructured by age.''   In particular, all live individuals are treated as being able to
reproduce.  Thus, there is no juvenile period where individuals are not yet sexually mature.
And similarly, there is no senescence period where individuals have aged out of the ability to reproduce.
This is not representative of the facts of life for many longer-lived animal species.  For instance
humans do not become sexually mature until they reach puberty.\footnote{There has been a 
decrease in the average age of onset of puberty over time.  For instance, a study reported that
the average age of onset of puberty in girls in the Germany was 16.6 years in 1860, but
this had fallen to 10.5 years in 2010.  Similar decreases have been observed in many 
other countries as well, including the U.S.  A parallel trend for boys has been observed, but with
ages about one or two years greater than for girls.  There is no clear consensus about the cause, although
improvements in nutrition and greater prevalence of certain chemicals in the environment
tied to production of sex hormones have been suggested.  See
{\tt https://www.theguardian.com/society/2012/oct/21/puberty-adolescence-childhood-onset}, accessed July 18, 2017.}
At the other end of the life cycle, human females typically live another 20 to 30 years after they reach menopause,
when they become infertile.\footnote{The average age of onset of menopause
is around 51 years for U.S. women.}  Human males typically remain
able to reproduce somewhat longer.  

If we consider what it would take to include features like a juvenile period or a senescence period
in a population model,  
then it is clear that the total population $P(n)$ at time $n$ should be stratified or structured into (at least) three groups,
the juveniles $J$, the mature individuals $M$ and the senescent individuals $S$ with 
$$P(n) = J(n) + M(n) + S(n).$$
In addition we might want to subdivide $J,M,S$ even further into the males and females within
each age cohort if there were significant differences between the two sexes.  This would 
require something like six dependent variables $J_m(n)$ for the male juveniles, $J_f(n)$ for the
female juveniles, etc.  
A population model would  then describe how each of these contribute to the corresponding
$J_m(n+1)$, $J_f(n+1)$, $M_m(n+1)$, $M_f(n+1)$, $S_m(n+1)$, and $S_f(n+1)$  at the next time step.  
Hence what we will be describing is 
a \emph{simultaneous system} of several first order difference equations for the age and gender strata of the
population as functions of the discrete time variable.  

In this chapter, we will study such structured population models in detail first, then
look at several other types of models given by systems of difference equations.
The models we will develop are also commonly called \emph{compartment models}
\index{model!compartment|(}
\index{compartment model|(}
since developing them will often involve analyzing how each of the different strata or compartments
at time $n$ contributes to the corresponding compartments at time $n + 1$.  

\specialsection{Age-structured Population Models}
\index{model!age-structured population|(}
\index{age-structured population model|(}
Let's begin with a very simple example that shows how to set up one of these compartment models.

\begin{example}
\label{FirstPopEx}
Instead of a long-lived animal species, let us consider a situation more typical of \emph{insects}
where the population is made up of only \emph{juveniles} and \emph{adults}.  We will also treat
only the \emph{female} members of the population.  We state our assumptions in words first,
then consider how to translate them into a more formal mathematical model.  We assume:
\begin{itemize}
\item Juveniles enter the juvenile population by way of hatching of eggs laid by the adults in the previous time step
\item Juveniles leave the juvenile population either by becoming adults, or by dying
\item Adults enter the adult population from the juvenile population when they mature
\item Adults leave the adult population by dying
\end{itemize}
At this level of detail we are not distinguishing between the possible causes of death.  This
could be due to predation, to ``natural causes,'' to diseases, or something else. 

It will be helpful to visualize all of this by drawing the following sort of \emph{compartment diagram},
where the boxes marked $J$ and $A$ represent the juvenile and the adult strata of the insect 
population, respectively.\footnote{You may notice a similarity between this sort of diagram and the
energy flow diagrams we studied in the Part I Summary Project from Chapter 3.  Those
were slightly different though, in that they did not contain the \emph{dynamic} aspect
of a change from one time step to the next.}

\setlength{\unitlength}{1mm}
\begin{figure}
\begin{picture}(120,30)
\put(10,5){\line(1,0){40}}
\put(50,5){\line(0,1){20}}
\put(10,5){\line(0,1){20}}
\put(10,25){\line(1,0){40}}
\put(70,5){\line(1,0){40}}
\put(110,5){\line(0,1){20}}
\put(70,5){\line(0,1){20}}
\put(70,25){\line(1,0){40}}
\put(20,5){\vector(0,-1){5}}
\put(80,5){\vector(0,-1){5}}
\put(50,10){\vector(1,0){20}}
\put(70,20){\vector(-1,0){20}}
\put(23,14){juveniles}
\put(84,14){adults}
\put(22,0){death}
\put(82,0){death}
\put(57,22){eggs}
\put(51.5,5){maturation}
\end{picture}
\caption{A schematic two-compartment population model.}
\label{fig:BasicDiagram}
\end{figure}

To turn Figure~\ref{fig:BasicDiagram} into a model, we need to make some assumptions.  
\begin{itemize}
\item Let us say that the deaths from the juvenile group make up some constant fraction $0< d_J < 1$ of whatever that population is. 
That is, the constant $d_J$ represents the {\it juvenile death rate per unit time}.  Therefore, the 
arrow marked death leading out of the $J$ box in Figure~\ref{fig:BasicDiagram} will represent a loss of $d_J\cdot J(n)$ individuals.
\item Similarly say the deaths from the adult group represent a loss of $d_A A(n)$ where $d_A$ is a constant adult death
rate per unit time.  Note: The juvenile death rate and the adult death
rate could conceivably be different\footnote{For instance, if the adults have better defenses against predators than the juveniles do.} so
we have given them different names to allow for this.
\item The arrow marked ``maturation'' represents the movement of (former) juveniles into the adult group when they mature.
Let us say that this proportion is some constant fraction $m$ of $J(n)$.  We do not necessarily want to assume that $d_J + m= 1$
since it might take more than one time step for the juveniles to mature; in other words the difference $1 - (d_J + m)$ will
represent the fraction of juveniles at time $t = n$ that remain juveniles at time $t = n+1$.  
\item Finally, the arrow marked ``eggs'' represents the number of juveniles produced by the adults laying eggs that hatch and 
produce new juveniles.\footnote{We are also assuming, in effect, that the time step is the time it takes for the eggs to hatch.}  We will assume
that each adult (female) lays some average number of eggs $b$ per unit time that hatch and produce juveniles.
\end{itemize}
Note that the $d_J, d_A, m, b$ parameters are all treated as \emph{constants} here.  That means in real-world terms
that we are assuming whatever challenges or opportunities these organisms are encountering from their habitat 
are \emph{not changing}.  Models of this form should probably \emph{only be used} when we can be reasonably 
sure that those assumptions
are not being violated or when projections over short time spans are sufficient.  

We record this more precise information in the updated diagram in Figure~\ref{fig:Diagram}
\setlength{\unitlength}{1mm}
\begin{figure}
\begin{picture}(120,30)
\put(10,5){\line(1,0){40}}
\put(50,5){\line(0,1){20}}
\put(10,5){\line(0,1){20}}
\put(10,25){\line(1,0){40}}
\put(70,5){\line(1,0){40}}
\put(110,5){\line(0,1){20}}
\put(70,5){\line(0,1){20}}
\put(70,25){\line(1,0){40}}
\put(20,5){\vector(0,-1){5}}
\put(80,5){\vector(0,-1){5}}
\put(50,10){\vector(1,0){20}}
\put(70,20){\vector(-1,0){20}}
\put(26,14){$J(n)$}
\put(86,14){$A(n)$}
\put(22,0){$d_J\cdot J(n)$}
\put(82,0){$d_A\cdot A(n)$}
\put(55,22){$b \cdot A(n)$}
\put(54,5){$m \cdot J(n)$}
\end{picture}
\caption{A simple two-compartment population model.}
\label{fig:Diagram}
\end{figure}
Now we are ready to turn our understanding of the flows into and out of the compartments into a model.  The idea
is that we total up \emph{all of the changes, both increases and decreases for each compartment}.  Referring to Figure~\ref{fig:Diagram},
we get the following
\begin{align}
\label{eq:TwoCChange}
\nonumber
J(n+1) - J(n) &= b\cdot A(n) - m\cdot J(n) - d_J\cdot J(n)\\
A(n+1) - A(n) &= m\cdot J(n) - d_A\cdot A(n).
\end{align}
This can be simplified and rearranged to the final form in \eqref{eq:TwoC}.
\begin{align}
\label{eq:TwoC}
\nonumber
J(n+1) &= (1 - m - d_J)\cdot J(n) + b\cdot A(n) \\
A(n+1) &= m\cdot J(n) + (1 - d_A)\cdot A(n) .
\end{align}
This is what we meant before by a system of simultaneous difference equations.  Each is linear and homogeneous and has 
constant coefficients according to the terminology introduced in Chapter 7.  Note that the 
values $J(n+1)$ and $A(n+1)$ both depend on both $J(n)$ and $A(n)$ in this case.  

As a final comment we mention that 
the constants $b,m,d_J,d_A$ will usually depend on the species that we are looking at and also on 
the habitat.  The constants $d_J, d_A$ would include contributions representing predation from 
other species, the average life span of the individuals, etc.  The constant $b$ is determined by how
successfully the species is reproducing; the constant $m$ is related to how long the
juvenile period lasts.  Experimental
data could be used to estimate these parameters and the process could be described as ``fitting'' a model of this form
to the experimental data.  This is somewhat analogous to what we did in Chapters 4, 5, and 6 to fit
linear, exponential, and power law models.  However, more sophisticated techniques are required
for this sort of fitting.    We will not discuss that aspect of developing such models in
detail.  $\triangle$
\end{example}

Generalizing what we did in Example~\ref{FirstPopEx}, we can incorporate any number of age strata.

\begin{example}
\label{MoreGeneralStructured}
A general population model for a population structured by age could be constructed in a similar fashion.
Say we wanted to model a human female population broken down into 7 age groups:
\begin{align*}
P_1(n) &= {\rm number\ between\ 0\ and\ 9\ years\ in\ age\ at\ time\ } n\\ 
P_2(n) &= {\rm number\ between\ 10\ and\ 19\ years\ in\ age\ at\ time\ } n\\
P_3(n) &= {\rm number\ between\ 20\ and\ 29\ years\ in\ age\ at\ time\ } n\\
P_4(n) &= {\rm number\ between\ 30\ and\ 39 \ years\ in\ age\ at\ time\ } n\\
P_5(n) &= {\rm number\ between\ 40\ and\ 49\ years\ in\ age\ at\ time\ } n\\
P_6(n) &= {\rm number\ between\ 50\ and\ 59\ years\ in\ age\ at\ time\ } n\\
P_7(n) &= {\rm number\ between\ 60\ years\ and\ older\ at\ time\ } n
\end{align*}
As in the previous example, let us assume that all birth and death rates are constant.
For simplicity, say the time step is 10 years, so if a woman counted in $P_5(n)$ survives
for that 10 years, she will be counted in $P_6(n+1)$.  Imagine a diagram as in Figure~\ref{fig:Diagram}, 
but now with 7 boxes, one for each age group, arranged in order from left to right.  For each box, we will have an arrow 
to the next one to the right corresponding to indivduals who survive and join the next age group.  Next, each box has an 
arrow representing deaths.  Say some fraction $d_i \cdot P_i(n)$ of the women in group $i$ \emph{do not survive}
to enter the $(i + 1)$st group, so $(1-d_i)\cdot P_i(n)$ do contribute to $P_{i+1}(n+1)$.  Finally, suppose 
that $b_i$ is the birth rate per unit time (10 years) for the women in group $i$.   This means that
we also have arrows back from the box for $P_i$ to $P_1$, given by $b_i P_i(n)$.  We expect
$b_1$ and $b_7$ to be very small, but probably not exactly zero.\footnote{With the decreasing average
age of onset of puberty for girls, births to mothers of age 9 and younger have been recorded(!)  Similarly,
a few women remain fertile past the usual average age of onset of menopause and techniques like
in vitro fertilization can be used to start pregnancies which then are brought to term after implantation of the
fetus in an older woman.  Births to mothers age 60 and above are possible, if extremely rare.} Finally, in the last
box for the 7th group, $(1 - d_7)\cdot P_7(n)$ represents the individuals who survive
from the 60 and older group.  We can think of that flow as an arrow from the box for $P_7(n)$ to itself.
Putting this all together we have
\begin{align}
\label{eq:GenStructModel}
\nonumber
P_1(n+1) &= b_1\cdot P_1(n) + b_2\cdot P_2(n) + b_3\cdot P_3(n) +  b_4\cdot P_4(n) + b_5\cdot P_5(n) \\ \nonumber
& \qquad + b_6\cdot P_6(n) + b_7\cdot P_7(n)\\ \nonumber
P_2(n+1) &= (1 - d_1)\cdot P_1(n)\\ \nonumber
P_3(n+1) &= (1 - d_2)\cdot P_2(n)\\
P_4(n+1) &= (1 - d_3)\cdot P_3(n)\\ \nonumber
P_5(n+1) &= (1 - d_4)\cdot P_4(n)\\ \nonumber
P_6(n+1) &= (1 - d_5)\cdot P_5(n)\\ \nonumber
P_7(n+1) &= (1 - d_6)\cdot P_6(n) + (1 - d_7)\cdot P_7(n).
\end{align}
Note that this is slightly different from the model in Example~\ref{FirstPopEx} because, except for the final group, no one
in group $i$ at time $t = n$ remains in group $i$ at $t = n+1$.  We have, for instance
$$P_2(n+1) - P_2(n) = -P_2(n) + (1 - d_1)\cdot P_1(n),$$
because every woman who was 10 to 19 years old at time $t = n$ has aged out of that group by $t = n+1$.  
The $-P_2(n)$ terms cancel, leaving the right-hand side from the second equation in \eqref{eq:GenStructModel}.  $\triangle$
\end{example}

Now that we have looked at techniques for formulating structured population models as systems of difference equations,
let's return to our first example and understand some of the behavior of solutions.

\begin{example}
\label{FirstPopExSols}
\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.6\linewidth]{DecayJ.png}
  \caption{The juveniles.}
  \label{fig:DecayJ}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.6\linewidth]{DecayA.png}
  \caption{The adults.}
  \label{fig:DecayA}
\end{subfigure}
\caption{The evolution of the juveniles and adults with $d_A = .85$ in Example~\ref{FirstPopExSols}.}
\end{figure}

(a)  To illustrate some possible behavior, let's consider the following (hypothetical) values of the 
constants in \eqref{eq:TwoC}:
$$b = 15, m = 0.05, d_J = .9, d_A = 0.85$$
As with a single difference equation, we need initial values to get the process started.  So
let's assume $J(0) = 0$ and $A(0) = 30$.\footnote{For instance, this could mean that 30 
mosquitoes enter a new habitat with no juveniles and start to breed there.}  
Starting from 
the initial conditions we compute:
\begin{align*}
J(1) &= (1 - 0.05 - 0.9)\cdot 0 + 15\cdot 30  = 450\\
A(1) &= 0.05\cdot 0 + (1 - 0.85)\cdot 30  = 4.5.
\end{align*}
Then with the values of $J(1)$ and $A(n)$ in \eqref{eq:TwoC} we find
\begin{align*}
J(2) &= (1 - 0.05 - 0.9)\cdot 4.5 + 15\cdot 4.5 +  = 90\\
A(2) &= 0.05\cdot 450 +  (1 - 0.85)\cdot 4.5 +  = 23.175,
\end{align*}
and so forth.  We can compute as many different steps as desired.
Let's not worry about the fact that we have fractional parts in some of the numbers here.  
We can compute the values out to $n = 40$  by using Excel in a way similar to what we did in 
\index{spreadsheet}
Chapter 7 to compute
approximate solutions of single difference equations.  
This can be done easily by setting up two parallel columns, say {\tt A} for the $J$ values, {\tt B} for the $A$ values.
We enter the initial conditions in row 1 in those columns.  Then in {\tt A2} enter the formula for computing
$J(1)$ in terms of $J(0)$ and $A(0)$:  {\tt = 0.05*A1 + 15*B1} and {\tt = 0.05*A1 + 0.15*B1}.  Then
copy and paste those formulas down the rest of those columns.  Since our addresses are not 
relative, the row numbers are increased by one in each subsequent row.  
You will see that 
there is an initial period of large oscillations from one time step to the next.  But that 
appears to die down, and the populations of juveniles and adults tend to lower and lower
values.  Figures~\ref{fig:DecayJ} and \ref{fig:DecayA} show the evolution over time in this case.\footnote{I used software other than 
Excel to generate these, of course.}  
Note that both portions of the population seem to be decaying to zero.

\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.6\linewidth]{GrowthJ.png}
  \caption{The juveniles.}
  \label{fig:GrowthJ}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.6\linewidth]{GrowthA.png}
  \caption{The adults.}
  \label{fig:GrowthA}
\end{subfigure}
\caption{The evolution of the juveniles and adults with $d_A = 0.75$ in Example~\ref{FirstPopExSols}.}
\end{figure}

\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.6\linewidth]{EquilJ.png}
  \caption{The juveniles.}
  \label{fig:EquilJ}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.6\linewidth]{EquilA.png}
  \caption{The adults.}
  \label{fig:EquilA}
\end{subfigure}
\caption{The evolution of the juveniles and adults with $d_A = 0.7895$ in Example~\ref{FirstPopExSols}.}
\end{figure}

(b) On the other hand, suppose we change $d_A$ to $0.75$ (meaning that more adults survive
to the next time step and hence more eggs are produced).   Figures~\ref{fig:GrowthJ} and 
\ref{fig:GrowthA} show the evolution over time in this case.  We see
what looks like exponential growth (after an initial oscillatory period).  

(c)  Finally, with $d_A = 0.7895$, the numbers of juveniles and adults appear
\index{equilibrium!solution}
to be settling down to equilibrium values  (again, after an initial oscillatory period).  
Figures~\ref{fig:EquilJ} and \ref{fig:EquilA} show the evolution over time in this case.  
We conclude this example by examining what we can say about equilibrium solutions of the system with
$d_A = 0.7895$. 
As we know from Chapter 7, equilibrium solutions should come by taking the difference equation and putting it in
the rearranged form giving the change in $Q$:
$$Q(n+1) - Q(n) = \text{\rm right-hand\ side}.$$
Then we setting the right-hand side equal to zero and solve.  If we do that with \eqref{eq:TwoCChange}, 
we get the simultaneous system:
\begin{align}
\label{eq:kernel}\nonumber
0 &= - 0.95\cdot J(n) + 15\cdot A(n)\\
0 &= 0.05\cdot J(n) - 0.7895\cdot A(n).
\end{align}
defining the equilbrium solutions.\footnote{These two equations are \emph{nearly} proportional.  
In fact given the values for $m, d_J, b$ that we were
using, there is a unique $d_A = \frac{15}{19} \doteq 0.7895$ for which this is true exactly.  In this 
situation, there are infinitely many different solutions of a system of equations like
\eqref{eq:kernel}.}  If we examine the computed solutions we see
$$J(100) \doteq 258.3 \quad {\rm and} \quad A(100) \doteq 16.4$$
are close to the apparent equilibrium values from Figures~\ref{fig:EquilJ} and \ref{fig:EquilA}.

It is easy to check that the right-hand sides of \eqref{eq:kernel} evaluate to numbers very close to zero for 
these $A$ and $J$ values.\footnote{Technical note:  It is possible to derive exact formulas for, and a complete quantitative and qualitative description of,
these solutions $J(n)$ and $A(n)$ using \emph{linear algebra}.  
\index{linear algebra}
What is needed is information about the \emph{eigenvalues} $\alpha_1,\alpha_2$ of the
matrices $\begin{pmatrix} .05 & 15 \\ .05 & 1-d_A\end{pmatrix}$ from the right-hand sides of \eqref{eq:TwoC}.  It turns out that
$J(n)$ and $A(n)$ have the form
$$J(n) = c_1 \alpha_1^n + c_2 \alpha_2^n \quad {\rm and} \quad A(n) = d_1 \alpha_1^n + d_2 \alpha_2^n,$$
where $c_1,c_2,d_1,d_2$ are constants depending on the initial conditions and on $d_A$.  As long as $\begin{pmatrix} c_1 \\ d_1\end{pmatrix}$
is an eigenvector for $\alpha_1$ and $\begin{pmatrix} c_1 \\ d_1\end{pmatrix}$ is an eigenvector for $\alpha_2$, then these
formulas for $J(n)$ and $A(n)$ define solutions of the system of difference equations.
When $d_A = 0.75$, one eigenvalue $\alpha_1 > 1$ and the other $\alpha_2 < 0$ with $|\alpha_2| < 1$, so we get approximate exponential growth for
large $n$.  The oscillatory behavior observed above comes from the alternating signs of $\alpha_2^n$.
When $d_A = 0.85$, $0 <\alpha_1<1$ and $\alpha_2 < 0$ with $|\alpha_2|<1$, so both $J(n)$ and $A(n)$
tend to $0$ for large $n$ and we get oscillations again.  Finally when $d_A = \frac{15}{19}$, $\alpha_1 = 1$ and $\alpha_2 <0$ with $|\alpha_2| < 1$, which
shows that there is an equilibrium which is approached as the oscillations from $\alpha_2^n$ die away.}   $\triangle$
\end{example}
\index{model!age-structured population|)}
\index{age-structured population model|)}

\specialsection{Epidemiological Models}
\index{epidemiological models|(}
\index{model!epidemiological|(}
In this section we will consider some models for the spread of diseases through a population of some species of organisms.
What we will see is that the general compartment setup can also be used to develop models of this sort.  Most of the 
models in this section will be \emph{non-linear}, though, because in order to spread, infectious diseases require some contact 
between an already-infected individual (or something that has been in contact with an infected individual\footnote{Spoiler Alert:  like the
doll in the film {\it War for the Planet of the Apes} from summer 2017.}) and an individual who is not yet infected.  We will see 
how modeling those sorts of interactions between populations can lead to non-linear terms in difference equation models.

\subsection*{The SIR Model}
\index{SIR model|(}
\index{model!SIR|(}
We begin with a simple example called the \emph{SIR model}.\footnote{The original form of this, a differential equation model,
was proposed by W.O. Kermack and A.G. McKendrick, ``A Contribution to the Mathematical Theory of Epidemics''. Proceedings of the Royal Society A. 115 (1927): 700--721.}
This epidemic model applies primarily to viral infections where the time scale for the infection and resolution (either by recovery or death)
is much shorter than the life times of the affected individuals and infection \emph{confers lasting immunity} in those who survive it.  Typical time
steps might be days or weeks, not years.  Examples
include measles, mumps, or rubella (``German measles'').  These are \emph{generally} non-fatal 
(especially in children).  The basic model also applies to infections
where death occurs in some fraction of those infected, but does not distinguish between recovery and death as an outcome.
It does not apply well to diseases like the common cold, which do not confer immunity.

\begin{assumptions}
\label{SIRAss}
The assumptions of the SIR model include the following.  
\begin{itemize}
\item The habitat is closed (recall from Chapter 7 that this means there is no in- or out-migration).  
Moreover, there is no important geographic or age structure of the population.  Everyone can come into 
contact with everyone else.\footnote{This basic model
was first used, in fact, to understand outbreaks of disease in environments such as boarding schools, so these
assumptions are more reasonable in that type of setting.  But treating a portion of the human population of the world
in isolation is unrealistic in most situations. With the greater mobility provided by easy, cheap air travel to virtually
all parts of the world, there are new challenges to controlling disease outbreaks!}
\item  The name comes from the fact that the population is divided into three groups, the \emph{susceptible population}, $S$, the \emph{infected population}, $I$,
and the \emph{removed population}, $R$.  Individuals are removed if they contract the disease and either recover or die from it.  
It is a feature of the model that the total population $P = S(n) + I(n) + R(n)$ does not depend on $n$, which means that any 
individuals who do die are still being counted in the value of $P$.   
\item Susceptible individuals are those who have not contracted the disease yet.  The basic SIR model does not
attempt to distinguish between the individuals in the $S$ group on the basis of better- or worse-functioning immune systems
or other resistance to the disease.  
\item Infected individuals are assumed to be contagious the whole time they are infected, but they remain in contact with
susceptibles.\footnote{That is, the model does not include any segregation or quarantine of infected individuals.}
 \item New infections occur because of contacts between infected and susceptible individuals.  The product $S(n) \cdot I(n)$ 
is a measure of the number of possible contacts.  Some constant fraction of them produce new infections.  That is, there is 
a constant $\beta$, $0 < \beta < 1$, such that the susceptibles are reduced by $\beta\cdot S(n) \cdot I(n)$ and 
the infecteds are increased by the same number.
\item \emph{Recovered} individuals acquire immunity and cannot be infected a second time.  (The last part is also true for
any individuals who die from the infection, of course!)
\item  The fraction of the infected group that is removed in each time unit is some constant $\gamma$ with $0 < \gamma < 1$.  
\end{itemize}
\end{assumptions}

These assumptions are summarized in the compartment diagram from Figure~\ref{fig:SIRDiagram}.
\setlength{\unitlength}{.6mm}
\begin{figure}
\begin{picture}(180,30)
\put(12,5){\line(1,0){36}}
\put(48,5){\line(0,1){20}}
\put(12,5){\line(0,1){20}}
\put(12,25){\line(1,0){36}}
\put(72,5){\line(1,0){36}}
\put(108,5){\line(0,1){20}}
\put(72,5){\line(0,1){20}}
\put(72,25){\line(1,0){36}}
\put(132,5){\line(1,0){36}}
\put(168,5){\line(0,1){20}}
\put(132,5){\line(0,1){20}}
\put(132,25){\line(1,0){36}}
\put(48,15){\vector(1,0){24}}
\put(108,15){\vector(1,0){24}}
\put(28,14){$S$}
\put(88,14){$I$}
\put(148,14){$R$}
\put(50,17){$\beta \cdot S\cdot I$}
\put(114,17){$\gamma \cdot I$}
\end{picture}
\caption{The SIR model in schematic form.}
\label{fig:SIRDiagram}
\end{figure}
Translating this into difference equations we get
\begin{align}
\label{eq:SIRModel}  \nonumber
S(n+1) &= S(n) - \beta\cdot S(n)\cdot I(n)\\
I(n+1) &= I(n) + \beta\cdot S(n)\cdot I(n) - \gamma\cdot I(n)\\ \nonumber
R(n+1) &= R(n) + \gamma\cdot I(n)
\end{align}

\begin{example}
\label{SIRExample}
Figure~\ref{fig:SIRPlot} shows the course of one epidemic modeled with SIR.  To generate this, I 
used the hypothetical parameter values $\beta = 0.004$ and $\gamma = 0.11$.  
\begin{figure}
\begin{center}
\includegraphics[width=3.0in,width=3.0in]{SIR.png}
\caption{A computed solution of SIR Model--$S(n)$ in red, $I(n)$ in black, $R(n)$ in blue. }
\label{fig:SIRPlot}
\end{center}
\end{figure} 
The initial conditions were $S(0) = 99$, $I(0) = 1$, $R(0) = 0$.   We have a large susceptible group
to start but the infection is introduced with $I(0)  = 1$.  With 
this combination of values for $\beta,\gamma$ the susceptible population $S(n)$ steadily decreases to what appears
to be an equilibrium value at a small positive number.  The infected population $I(n)$ increases to
a maximum at about $n = 21$, then decreases to values close to zero.  Finally, the removed
population $R(n)$ steadily increases until it levels off at something a bit than $100$.  The
epidemic has died out but a small number of susceptible individuals are left.  $\triangle$
\end{example}

The size of the removed
population $R(n)$ is also effectively the cumulative number of infections in Figure~\ref{fig:SIRPlot}.  The cumulative
number of infections in the 2014-2015 Ebola outbreak studied in Chapter 5 looks something like this.
\index{Ebola outbreak 2014-15}

\subsection*{Observations About SIR}
Because of the non-linear $S(n)\cdot I(n)$ terms in two of the equations in \eqref{eq:SIRModel}, we should
not expect to find exact formulas for the solutions of SIR models.  However, we can make several qualitative
observations based on the form of the equations.  
\begin{itemize}
\item The constants $\beta$ and $\gamma$ have definite
interpretations so their values could be estimated by observing the properties of an infectious disease 
even without some of the more advanced model fitting techniques to which we alluded above.   
\item The constant $\gamma$ measures the \emph{rate at which individuals move from the infected
group to the removed group}.  Hence $\frac{1}{\gamma}$ is a measure of \emph{how long an individual
stays infected} -- the total time from infection to entry into the removed group.
\item The constant $\beta$ measures \emph{how infectious the disease is} -- how easily it spreads from contact.
A high value of $\beta$ would indicate a very infectious disease that would spread quickly; a 
small value of $\beta$ would indicate a disease that is harder to catch in the sense that a smaller
fraction of contacts produce infections.
\item The second equation in \eqref{eq:SIRModel} can be 
rearranged to the form:
$$I(n+1) - I(n) = (\beta\cdot S(n) - \gamma)\cdot I(n)$$
This says that $I(n+1) > I(n)$ exactly when $\beta\cdot S(n) - \gamma > 0$, or equivalently
\begin{equation}
\label{eq:BasicReproductionNumber}
\framebox{$\displaystyle{I(n+1) > I(n) \quad \Longleftrightarrow \quad S(n) > \frac{\gamma}{\beta}}$.}
\end{equation}
This has several important consequences.  If $S(0) > \frac{\gamma}{\beta}$ 
then the number of infected individuals will increase at first and this is the formal way
\index{epidemic (definition)}
an \emph{epidemic} is defined.  If $S(0) \le \frac{\gamma}{\beta}$, on the other hand, then the number of infected
individuals is decreasing from the start and no epidemic will occur.   For this reason, epidemiologists
refer to the constant 
$$R_0 = \frac{\beta \cdot S(0)}{\gamma}$$
as the \emph{basic reproduction number} for an infectious disease outbreak.  Values $R_0 > 1$ 
\index{basic reproduction number (SIR)}
indicate there will be an epidemic because $R_0 > 1$ implies $S(0) > \frac{\gamma}{\beta}$ as in 
\eqref{eq:BasicReproductionNumber}.  Values $R_0 \le 1$ mean the outbreak will die out without producing
an epidemic.   Similarly, since $S(n)$ is steadily decreasing, we expect that as long as new
infections occur, it will eventually be true that  $S(n) = \frac{\gamma}{\beta}$.  At this point, 
the change from $I(n)$ to $I(n+1)$ becomes \emph{negative} and the outbreak begins
to subside.  
\end{itemize}

\begin{example}
\label{SIRExampleComment}
In the epidemic modeled in Example~\ref{SIRExample}, in Figure~\ref{fig:SIRPlot} we can
see that $I(n)$ begins to decrease at about $n = 21$.  We used $\beta = 0.004$ and $\gamma = 0.11$
there, so $\frac{\gamma}{\beta} = 27.5$.  In Exercise~\ref{ex:SIRExampleVerify}, you will check
that $S(20) \doteq 29.8$ and $S(21) \doteq 25.2$.  Therefore, $n = 21$ is indeed the first
time that $S(n)$ drops below the threshold value $\frac{\beta}{\gamma} = 27.5$ and
that is where the ``corner is turned'' and this epidemic begins to subside. 
\index{threshold} $\triangle$
\end{example}

\subsection*{Extensions and Modifications}
The basic SIR model can be modified to include other features and take other factors into account.
This yields a large and growing class of different epidemiological models that are used to 
study the progress of disease outbreaks and manage public health responses.

\begin{itemize}
\item One very basic modification would be to treat the recovered individuals separately from 
the individuals who die as a result of a potentially fatal infectious disease.  You will 
develop a modified SIR model of that form in Exercise~\ref{ex:SIRRecoveredDead}.
\item The SIR model could also be modified to account for population dynamics from births
and deaths from other causes.  This just amounts to adding new arrows in and out of the compartments
in Figure~\ref{fig:SIRDiagram}.
\item To remedy some of the unrealistic features of SIR from Assumptions~\ref{SIRAss}, we 
know some diseases strike different age groups within a population in different ways because
of differences in immunity, differences in overall health, etc.  The basic SIR model can be 
combined with the techniques for age-structured populations that we developed in \S 2 of this chapter
to study epidemiological questions within structured populations.  You will see a simple example in Exercise~\ref{ex:SIRAge}.
\item The SIR model can also be used to study the effects of vaccination on the dynamics
of disease outbreaks.  Recall that these models apply mainly to diseases where, in many cases, 
humans have learned to trick our immune systems into mounting the same defense responses
acquired through actual exposure to the disease without undergoing the full exposure.  
You will see some examples in Exercise~\ref{ex:SIRVaccination}.
\index{vaccination}
\item More realistic models of the interaction between susceptibles and infecteds than the 
simple product term $S(n) \cdot I(n)$ have also be studied.  
\end{itemize}

Epidemiological modeling is a huge subject and we have only scratched the surface.
However, this limited exposure should convince you that the techniques we have
developed have applications of great interest in public health questions.  
\index{SIR model|)}
\index{model!SIR|)}
\index{epidemiological models|)}
\index{model!epidemiological|)}

\specialsection{Predator-Prey Models}
\index{predator-prey model|(}
\index{model!predator-prey|(}
The difference equation models that we have been studying in this chapter 
are very versatile and powerful tools for studying all sorts of dynamical problems.  In this final section, we
want to illustrate some of their power by modeling interactions between two different species.   A well-known
real world data set is the 90-year record of pelts of snowshoe hare and lynx bought by the Hudson's Bay 
Company in Canada.  Human hunters caught animals of both species and sold their fur to the company, which
recorded these numbers.   But in the wild, the lynx (large cats) are predators who rely on the hares (large rabbits)
for much of the food they consume and the hares are herbivores.  The data is shown 
in Figure~\ref{fig:LynxHareData}.\footnote{From {\tt http://jan.ucc.nau.edu/lrm22/lessons/predator$\_$prey/predator$\_$prey.html}, 
accessed July 20, 2017.}
\begin{figure}
\begin{center}
\includegraphics[width=3.0in,width=3.0in]{hare_lynx_graph.png}
\caption{Hare and Lynx Pelts Bought by the Hudson's Bay Company.}
\label{fig:LynxHareData}
\end{center}
\end{figure}
Now the idea is that the numbers here should reflect something about the \emph{total populations} of the 
hares and the lynxes since the human hunting efforts were roughly comparable over the whole time.  
Note that there seem to \emph{approximately cyclical variations} (together with a lot of ``noise''
that can probably be explained by variations in weather, random variations in the hare and lynx
catches, and so forth) going on in this data set.  

Our goal is to model the interaction between predators and prey.  In rough terms, we can see that
the presence of large numbers of prey animals should stimulate the growth of the population of predators
(probably with some time delay in the response).  Then in the next phase, large numbers of predators should drive down
the prey population to the point where the large numbers of predators cannot be supported any longer.  At this
point the predator population should hit a maximum and then decrease.  
This will allow the prey population to recover and start another cycle.   Any model we develop should
ideally show something like this sort of cyclical variation, at least in some circumstances.  We will consider 
a pair consisting of a \emph{predator} species and the \emph{prey} species that they feed on.  Models of the same form
can also be used for the interactions between parasites and hosts, between herbivores and the plant species they
feed on, etc.\footnote{The model we study here goes back to differential equation models introduced by Alfred Lotka and
Vito Volterra in the 1920's.}
\index{Lotka, Alfred}
\index{Volterra, Vito}

\begin{assumptions}
\label{LVAss}
To formulate the basic form of our model, we make the following assumptions.  We will denote by $P(n)$ the 
population of predators at time $t = n$ and by $Q(n)$ the population of prey individuals at time $n$.  
\begin{itemize}
\item The predator species is entirely dependent on the prey species for food.  The prey species, on the 
other hand, is entirely supported by other resources available in the habitat.
\item The predator species population would follow an \emph{exponential decay model} in the 
\index{exponential decay}
absence of any of the prey.
\item The habitat has ample resources for the use of the prey and the 
population of the prey species would follow a \emph{exponential growth model}
\index{exponential growth}
as in Chapter 7 in the absence of any of the predators.\footnote{As we know from Chapter 5, this is not realistic in the long run.  However, 
these models can also be modified easily enough to make the prey species follow a logistic model if we wish.  We will consider
that sort of modification shortly.}
\item Contacts between predators and prey are modeled, as in the SIR model, by product terms $P(n)\cdot Q(n)$.
\end{itemize}
\end{assumptions}

These assumptions correspond to the following pair of first order, nonlinear difference equations where $0 < a < 1$ and $b,c,d > 0$ are constants:
\begin{empheq}[box=\fbox]{align}
\label{eq:LV}
\nonumber
P(n+1) &= (1 - a) \cdot P(n) + b\cdot P(n)\cdot Q(n)\\
Q(n+1) &= (1 + c)\cdot Q(n) - d\cdot P(n)\cdot Q(n)
\end{empheq}
We will call this system of difference equations the \emph{basic predator-prey model}.

Note that if $Q(n) = 0$ for all $n$, then there are none of the prey.  In that case $P(n+1) = (1 - a)\cdot P(n)$ is an exponential
decay model and the predator population will go to zero as $n$ increases.  The closer $a$ is to $1$, the 
faster that exponential decay will be.  On the other hand, if $P(n) = 0$ for all $n$, 
then $Q(n+1) = (1 + c)\cdot Q(n)$.  Since we assume $c > 0$, this gives exponential growth of the prey population and the 
larger $c$ the faster the growth will be.  The constant $b$ in effect
measures the value of the prey to the predators.  Other things being equal, the larger $b$ is, the more the prey
will contribute more to the growth of the predator population.   The constant $d$, on the other hand, measures 
how effective the predators are as predators.  Other things being equal, the larger $d$ is, the more the predators
will contribute to the reducing the growth of the prey population.

\begin{example}
\label{FirstLV}
We consider the predator-prey model from \eqref{eq:LV} using the parameters $a = 0.4$, $b = 0.004$, $c = 0.02$, $d = 0.0006$
and initial conditions $P(0) = 4$, $Q(0) = 100$.\footnote{Technical note: For those familiar with the solutions of the Lotka-Volterra differential equation
model, the sensitivity of the difference equation version to changes in these parameters and the variety of different sorts of 
solutions one can obtain can be quite surprising!}  The solutions are plotted together in Figure~\ref{fig:LVPlot}.
\begin{figure}
\begin{center}
\includegraphics[width=3.0in,width=3.0in]{LVPlot.png}
\caption{A computed solution of the predator-prey model--predators $P(n)$ in red, prey $Q(n)$ in blue. }
\label{fig:LVPlot}
\end{center}
\end{figure}
We can note several features of the solutions.  First the overall trend seems to be a \emph{cyclic variation} superimposed
on an increasing trend in the heights of the peaks of both populations.  This is perhaps somewhat surprising the first
time you encounter solutions of the systems of this form.  However, it does capture the some aspect of the cyclic
variations that are observed in predator-prey interactions in the real world as in Figure~\ref{fig:LynxHareData}(!) 
 
 A qualitative description of 
why the solutions of \eqref{eq:LV} exhibit this behavior is rather easy to provide.  When the number of predators is low, 
the growth of the prey is nearly exponential.  However, as the prey population grows, so does the $b\cdot P(n)\cdot Q(n)$
term in the first equation in \eqref{eq:LV}.  This produces a rise in the population of the predators in turn as we see
in the ``spikes'' of the predator near $n = 40, 125$, and $225$ in Figure~\ref{fig:LVPlot}.  Around those
spikes the $-d\cdot P(n)\cdot Q(n)$ terms are ``outweighing'' the exponential growth term in the second equation
in \eqref{eq:LV} and the prey population is decreasing.  But then the predators also decline to near zero and a new
phase of growth in the prey population begins.  $\triangle$
\end{example}

Even when the behavior of one of these solutions is relatively well-behaved over the short run (as in this example), 
it is possible for $Q(n)$ to attain negative values at the low end of one of the oscillations, and then the solution can ``go haywire'' 
in spectacular ways.  This happens in fact for $n = 875$ in our example and the results stop being realistic, or even
realistically computable, shortly thereafter(!))  To avoid behavior of this type, and also (and probably more importantly!)
to include more realistic assumptions in these models, it is also possible to build a \emph{logistic growth model} for
the prey into the set-up by considering \emph{modified predator-prey models} of the form:
\begin{empheq}[box=\fbox]{align}
\label{eq:LVLog}
\nonumber
P(n+1) &= (1 - a) \cdot P(n) + b\cdot P(n)\cdot Q(n)\\
Q(n+1) &= (1 + c)\cdot Q(n) - \frac{c\cdot Q(n)^2}{M} - d\cdot P(n)\cdot Q(n)
\end{empheq}
When $P(n) = 0$ for all $n$, the second equation becomes equivalent to \eqref{eq:Logistic}, with $M$ representing
the carrying capacity for the prey.
\index{logistic growth model!in predator-prey model}
\index{carrying capacity}

\begin{example}
\label{FirstLVLog}
We use the same parameter values $a = 0.4$, $b = 0.004$, $c = 0.02$, $d = 0.0006$
and initial conditions $P(0) = 4$, $Q(0) = 100$ as in Example~\ref{FirstLV}, but introduce the new term $-\frac{0.02\cdot Q(n)^2 }{300}$
with $M = 300$
in the second equation as in \eqref{eq:LVLog}.  The computed solutions are plotted in Figure~\ref{fig:LVLogisticPlot}.
\begin{figure}
\begin{center}
\includegraphics[width=3.0in,width=3.0in]{LVLogisticPlot.png}
\caption{A computed solution of the modified predator-prey model \eqref{eq:LVLog}--predators $P(n)$ in red, prey $Q(n)$ in blue. }
\label{fig:LVLogisticPlot}
\end{center}
\end{figure}

As before, when the number of predators is low, 
the prey population is growing (but approximately according to a logistic model rather than exponentially).  
However, as the prey population grows, so does the $b\cdot P(n)\cdot Q(n)$
term in the first equation in \eqref{eq:LVLog}.  This produces a rise in the population of the predators in turn as we see
in the ``spikes'' of the predator near $n = 50, 140$, and $230$ in Figure~\ref{fig:LVLogisticPlot}.  Around those
spikes the $-d\cdot P(n)\cdot Q(n)$ terms are ``outweighing'' the growth term in the second equation
in \eqref{eq:LVLog} and the prey population is decreasing.  But then the predators also decline to near zero and a new
phase of growth in the prey population begins.
Here the oscillations in the predator and prey populations are ``tamped down'' and they even appear to be decreasing in 
amplitude as the cycles continue.  

A different way to visualize what is happening here is to plot the data points
\index{phase portrait}
$(P(n),Q(n))$ for $n = 0, \ldots, 300$.  The result is called the \emph{phase portrait} of the system of difference equations
and is shown in Figure~\ref{fig:LVLogisticPhase}.
\begin{figure}
\begin{center}
\includegraphics[width=3.0in,width=3.0in]{LVLogisticPhase.png}
\caption{Phase portrait of the modified predator-prey model, solution from Figure~\ref{fig:LVLogisticPlot}.}
\label{fig:LVLogisticPhase}
\end{center}
\end{figure}
In the phase portrait, the explicit dependence on the time $t = n$ is not shown, 
but we can see that as $n$ increases
the point $(P(n),Q(n))$ is spiraling clockwise and moving inward in the first quadrant of the $P,Q$-plane.  $\triangle$
\end{example}

Considering what could be happening in Figure~\ref{fig:LVLogisticPhase}, it is natural to ask:  Will the phase portrait
continue to spiral, and if so, what is it spiraling toward?  It is at least somewhat intuitively clear that if the 
spiraling continues for all $n$, then the phase portrait form of the solution is either tending toward some 
oval-shaped closed curve inside the region we plotted in Figure~\ref{fig:LVLogisticPhase}, or else
toward a point $(P_0,Q_0)$ defined by an equilibrium solution of the system \eqref{eq:LVLog}. 
\index{equilibrium!solution}
In fact, we can see that there is an equilibrium point in the interior of the first quadrant for all systems of the
form either \eqref{eq:LV} or \eqref{eq:LVLog}.  To find equilibria of \eqref{eq:LV}, we would solve
\begin{align*}
0 &= -a\cdot P(n) + b\cdot P(n)\cdot Q(n) = P(n)\cdot (-a + b\cdot Q(n))\\
0 &= c\cdot Q(n) - d\cdot P(n)\cdot Q(n) = Q(n) \cdot (c - d\cdot P(n)).
\end{align*}
The first equation is zero if either $P(n) = 0$ or $Q(n) = \frac{a}{b}$.  The second is zero if either $Q(n) = 0$ or $P(n) = \frac{c}{d}$.
This means there are two equilibrium points:  $(0,0)$ where both predator and prey populations are zero, and the nonzero 
equilbrium
$$
\framebox{Equilibrium of \eqref{eq:LV} at $\displaystyle{(P_0,Q_0) = \left(\frac{c}{d},\frac{a}{b}\right)}$}
$$
in the first quadrant.   You will carry out the corresponding computation for the modified system \eqref{eq:LVLog} in Exercise~\ref{ex:LVLogEquil}.
The result is that there are now three equilibria: two at the equilibria of the logistic model for the prey in the absence of the 
predators: $(0,0), (0,M)$, and 
\begin{equation}
\label{eq:LVLogEquil}
\framebox{Equilibrium of \eqref{eq:LVLog} at $\displaystyle{(P_0,Q_0) = \left( \left(1 - \frac{a}{b\cdot M}\right)\cdot\frac{c}{d}, \frac{a}{b}\right)}$.}
\end{equation}

\begin{example}
\label{FirstLVLogEquil}
In Example~\ref{FirstLVLog}, using the values of the parameters given there:
 $a = 0.4$, $b = 0.004$, $c = 0.02$, $d = 0.0006$ and $M = 300$, according to \eqref{eq:LVLogEquil}, the
equilibrium point in the interior of the first quadrant is located at 
$$\left( \left(1 - \frac{0.4}{0.004\cdot 300}\right)\cdot \frac{0.02}{.0006}, \frac{0.4}{0.004}\right) \doteq (22.2,100).$$
This is located inside the beginning of the spiral in Figure~\ref{fig:LVLogisticPhase}.  If we plot more points, 
we see that it is apparently the case the solution is tending toward that equilibrium point as $n$ increases
without bound.\footnote{Technical note:  This is rather awkward to explore using software like Excel since we need to go
out to about $n = 1000$ to really start to see what is happening; the spirals
are moving in quite slowly.  Moreover that doesn't in itself show that the 
solution really does approach the equilibrium point, at least not to the level of evidence required
by (pure) mathematicians(!)  There is a theory of \emph{linearization} around equilibria of nonlinear systems
and asymptotic stability criteria that imply this really is what is happening in this case. The eigenvalues of 
the linearized system at the equilibrium point are complex numbers with negative real parts, which implies that the solutions
tend toward the equilibrium.}  $\triangle$
\end{example}

Other sorts of interactions between two or more species can be modeled using similar techniques.  We
will not pursue those models in this elementary course, though.  
\index{predator-prey model|)}
\index{model!predator-prey|)}


\specialsection{Chapter Project}

\subsection*{Background}  Our chapter project for this chapter involves modeling the (short-term) global \emph{carbon cycle}--the 
\index{carbon cycle|(}
movements of carbon among the various living and nonliving components of the terrestrial environment.  This 
includes processes acting over both long and short time scales.  For 
instance, there are large stores of carbon in the rocky subsurface layers of the Earth's crust (in the forms of coal, oil
and natural gas) that were created over long geological time spans by chemical processes acting on the tissues of plants
that lived millions of years ago.  We will not consider the process by which such deposits are formed because
they take place over such long time scales.  We will consider processes that take place over time spans of years 
or decades, though.  This is the \emph{short-term} part of the cycle.    

The main reason this is important is that this short-term cycling of carbon is necessary for
the continuation of life on Earth.  This is true because carbon makes up about $50\%$ of the \emph{dry weight} of most plant and animal tissues 
(that is, what is left after the water is taken away).  Anyone who has taken high school biology has surely marveled at
the way plants and animals have co-evolved to use carbon in complementary ways.  Animals breathe in oxygen from the atmosphere and exhale
carbon dioxide as a product of their metabolisms from consuming plant tissues as food; they also return carbon to the atmosphere in other ways.   
On the other hand, plants take in carbon dioxide and incorporate it into their tissues by means of photosynthesis, while releasing oxygen back 
into the atmosphere as a byproduct.  

As we mentioned in the Chapter Project on the Mauna Loa atmospheric carbon dioxide measurements data set from Chapter 4,
\index{Mauna Loa $CO_2$ data}
in addition to being the primary constituent of living tissue, carbon is also present in the atmosphere (mostly in the form of carbon dioxide).
Carbon dioxide is also an important \emph{greenhouse gas} and the current life forms present could not exist on Earth without its warming effect.  However,
for the past 200 years or so humans have been increasing the atmospheric carbon dioxide content by means of fossil fuel burning
and other activities.  How will that additional carbon dioxide affect the existing mechanisms of the carbon cycle, Earth's climate and in turn its life forms? 
Exactly how the carbon cycle interacts with climate is a very difficult question and we will not attempt to model that because
there are too many components and their interactions are only starting to be understood by scientists.  

We will concentrate on understanding the carbon cycle itself and starting to understand how some human activities 
can create changes. A BIG disclaimer is certainly in order here:  This is definitely a ``toy" model that is much simpler than the real world 
and at the same time much simpler than the climate models that scientists are currently using to try to understand the evolution of the 
Earth's climate under the influence of anthropogenic sources of atmospheric carbon dioxide.  For that reason, you should not take any of the 
computed values as especially realistic predictions.  

\subsection*{The Pre-Industrial-Revolution Carbon Cycle}
Figure~\ref{fig:PreIndustrial} shows an estimate of what the yearly carbon cycle looked like before about 1750 CE
(this is an arbitrary date meant to predate the beginnings of large-scale industrial activity by humans making use of fossil fuel burning).   The
boxes show the major \emph{carbon reservoirs} in the Earth system. The numbers in the boxes are total carbon contents in units of 
Gigatons $ =10^9$ metric tons
$ = 10^{12}$ kilograms $=10^{15}$ grams.\footnote{This is also called a ``Petagram."}  The arrows show the \emph{flows}, or ``fluxes'' of 
carbon between the different component systems, in units of $10^{12}$ kilograms \emph{per year}.  
\setlength{\unitlength}{1mm}
\begin{center}
\begin{figure}
\begin{picture}(130,100)
\put(12,80){\line(1,0){100}}
\put(12,80){\line(0,1){20}}
\put(12,100){\line(1,0){100}}
\put(112,80){\line(0,1){20}}

\put(12,40){\line(1,0){40}}
\put(12,40){\line(0,1){20}}
\put(12,60){\line(1,0){40}}
\put(52,40){\line(0,1){20}}

\put(72,40){\line(1,0){40}}
\put(72,40){\line(0,1){20}}
\put(72,60){\line(1,0){40}}
\put(112,40){\line(0,1){20}}

\put(12,0){\line(1,0){100}}
\put(12,0){\line(0,1){20}}
\put(12,20){\line(1,0){100}}
\put(72,0){\line(0,1){20}}
\put(112,0){\line(0,1){20}}

\put(47,89){Atmosphere -- 600}
\put(15,49){Land Life/Soils -- 2190}
\put(79,49){Oceans -- 40,000}
\put(25,9){Lithosphere -- 10,000,000}
\put(75,9){Oil,Coal,Gas -- 10,000}

\put(25,60){\vector(0,1){20}}
\put(40,80){\vector(0,-1){20}}
\put(3,71){Respiration/}
\put(3,66){decay -- 60}
\put(42,74){Photo-}
\put(42,69){synthesis -- }
\put(42,64){60}
\put(85,60){\vector(0,1){20}}
\put(100,80){\vector(0,-1){20}}
\put(103,71){Chemical exchanges --}
\put(103,66){90 (each way)}
\put(65,20){\vector(0,1){60}}
\put(68,29){Weathering -- .05}
\put(32,40){\vector(0,-1){20}}
\put(9,29){Burial -- .05}
\end{picture}
\caption{(Simplified) Carbon Cycle as of roughly 1750.}
\label{fig:PreIndustrial}
\end{figure}
\end{center}
This is a simplified schematic picture\footnote{Diagrams inspired by similar figures from the web site for a course at U. Arizona -- 
{\tt http://www.atmo.arizona.edu/students/courselinks/fall16/atmo336/lectures/sec3/carbon.html}, accessed July 21, 2017.} 
because it does not take into account the fact that 
there are effectively two ``layers'' in the oceans, an upper layer where most aquatic
life is located and where some carbon is dissolved into the water from the atmosphere, plus a lower layer
containing much less life, but also large carbon-bearing sediments on the ocean bottoms.  There
are fluxes between those two layers as well that are thought to be possibly important for the dynamics of the 
carbon cycle over longer time scales.  We are effectively ``lumping those together.''

\subsection*{Questions}  
\begin{enumerate}
\item[(A)] Assume the fluxes in Figure~\ref{fig:PreIndustrial} corresponding to Respiration/Decay, Photosynthesis, and the 
two-way Chemical Exchanges between the Atmosphere and the Oceans reservoirs, are all \emph{proportional to the amounts in
the source reservoir}, with proportionality constants determined by the values given.  However the fluxes corresponding to Weathering and Burial
are constant.   Formulate a system of four difference
equations equivalent to the diagram.  It is thought that the carbon cycle was very nearly \emph{at an equilibrium} before 1750.  Is
\index{equilibrium!solution}
that true for the given data for your model? \vskip 5pt

\item[(B)]  Now consider Figure~\ref{fig:PostIndustrial} which shows estimates for the amounts of carbon in the various
reservoirs and the fluxes \emph{as of about 1990 CE}, so almost up to the present.  Note this includes additional fluxes related to 
human activities that were not present before 1750 CE.\footnote{The numbers for pre-1750 and for 1990 are from different sources and were probably estimated
using different methods.  Hence they might not be completely consistent with each other.}
\setlength{\unitlength}{1mm}
\begin{center}
\begin{figure}
\begin{picture}(130,100)
\put(12,80){\line(1,0){100}}
\put(12,80){\line(0,1){20}}
\put(12,100){\line(1,0){100}}
\put(112,80){\line(0,1){20}}

\put(12,40){\line(1,0){40}}
\put(12,40){\line(0,1){20}}
\put(12,60){\line(1,0){40}}
\put(52,40){\line(0,1){20}}

\put(72,40){\line(1,0){40}}
\put(72,40){\line(0,1){20}}
\put(72,60){\line(1,0){40}}
\put(112,40){\line(0,1){20}}

\put(12,0){\line(1,0){100}}
\put(12,0){\line(0,1){20}}
\put(12,20){\line(1,0){100}}
\put(72,0){\line(0,1){20}}
\put(112,0){\line(0,1){20}}

\put(47,89){Atmosphere -- 750}
\put(15,49){Land Life/Soils -- 1975}
\put(79,49){Oceans -- 42,500}
\put(25,9){Lithosphere -- 10,000,000}
\put(75,9){Oil,Coal,Gas -- 9,350}

\put(30,60){\vector(0,1){20}}
\put(40,80){\vector(0,-1){20}}
\put(10,71){Respiration/}
\put(10,66){decay -- 60}
\put(42,74){Photo-}
\put(42,69){synthesis -- }
\put(42,64){60}
\put(73,60){\vector(0,1){20}}
\put(110,80){\vector(0,-1){20}}
\put(74,71){Chemical exchanges --}
\put(74,66){90 (each way)}
\put(65,20){\vector(0,1){60}}
\put(68,29){Weathering -- .05}
\put(32,40){\vector(0,-1){20}}
\put(9,29){Burial -- .05}

\put(112,10){\line(1,0){4}}
\put(116,50){\oval(8,80)[r]}
\put(116,90){\vector(-1,0){4}}
\put(122,33){\rotatebox{90}{fossil fuel burning -- 6.5}}

\put(8,90){\vector(1,0){4}}
\put(8,70){\oval(8,40)[l]}
\put(12,50){\line(-1,0){4}}
\put(0,56){\rotatebox{90}{deforestation -- 1.5}}

\end{picture}
\caption{(Simplified) Carbon Cycle, as of roughly 1990.}
\label{fig:PostIndustrial}
\end{figure}
\end{center}

Construct a new difference equation model for this regime, assuming that the fluxes (except for Burial and Weathering) 
that were present
before are all proportional to the amount in the source reservoir as of 1990 with proportionality constants
computed from the 1990 values of the fluxes.   \emph{However,  the new, anthropogenic fluxes are determined like this:
The fossil fuel burning and deforestation fluxes have been increasing by about $3\%$ per year on average during
this period}.  
For instance, this means that the difference equation for the Atmosphere reservoir should contain terms like
this corresponding to the flux in from the Oceans reservoir and the contribution from fossil fuel burning:
$$A(n+1) = A(n) + \frac{100}{41000}\cdot O(n) + \cdots + 6.5 \cdot (1.03)^n,$$
where $n$ is the number of years after 1990.  
There will be other terms as well, indicated by the $\cdots$.  
Using the 1990 values of carbon content of the reservoirs as initial conditions, 
\index{spreadsheet}
compute a solution for your model out to 2020 using Excel or suitable software.  \vskip 5pt

\item[(C)]  An atmospheric carbon amount of 600 Petagrams corresponds to an atmospheric carbon dioxide concentration
of about 283 ppm (parts per million), and the atmospheric amounts of carbon are always proportional to the carbon dioxide concentration with the 
same ratio.  From your computed solution in part (B), determine what your model says about the atmospheric carbon dioxide 
concentrations in ppm out to 2020.  How does this compare with the Mauna Loa measurements we discussed in Chapter 4?  
\index{Mauna Loa $CO_2$ data}
Is the model prediction greater or smaller than the real-world data? \vskip 5pt

\item[(D)]  (A ``research/thought question.'')  Is the assumption that the fluxes are proportional to the amounts in the source reservoir somewhat
reasonable (at least in the short-run)?  One possible argument for arguing that the fluxes should at least be increasing as the 
amount of carbon in the Atmosphere reservoir increases would be this:  More atmospheric carbon dioxide means
more nutrients and better growing conditions for many plants, which means more food for animals to eat, which means more respiration and decay.
Similarly, unless the Oceans are already saturated with carbon dioxide, the more there is in the Atmosphere, the more they can absorb
and the more dissolved carbon dioxide there is, the more will be released.  Does that seem to be consistent with the changes between
Figure~\ref{fig:PreIndustrial} and Figure~\ref{fig:PostIndustrial}?
Are there other factors that might affect this, though? \vskip 5pt

\item[(E)]  What is past is past and we cannot change it.  But we can use models\footnote{Preferably more realistic ones, of course.} 
to try to evaluate the effects of changes we might make in the future.  
\begin{enumerate}
\item[(1)] For instance, we might ask:  Suppose we were able to limit fossil 
fuel burning to a constant level of 10 Gigatons per year and eliminate carbon emissions due to deforestation and other changes in land use
entirely.  (This is still significantly higher than the Kyoto Accords target of $5.7$ Gigatons per year.)
What would happen?   To see,  the easiest way (if you have not done this already)
will be to create a new columns in your spreadsheet for the contributions due to those factors and use the values from
them to compute the other quantities.  You can manually change all values to 10 for the fossil fuel burning starting from any 
point.
\item[(2)]  What if we continue tracking the model into the future keeping fossil fuel burning at 10 Gt per year?  
Does the CO2 concentration look like it will ever return to current levels?  If so how long does it take? 
\item[(3)]  What is the largest constant (nonzero!) level for fossil fuel burning starting in 2017 that would still yield decreasing atmospheric
CO2 levels by the year 2100?  (This will require some experimentation!)
\item[(4)]  Model validation ``thought question'' -- Does your answer to (3) seem reasonable?  Is there something going on in this model that 
might not be that realistic over these time scales?  (Hint:  Look at the Land Life/Soils figures.  What would it take in real terms to have the 
amount of CO2 taken up by plants in photosynthesis increase by a factor of 10?)
\item[(5)]  What are some of the potentially important features of the real world carbon cycle that are being left out of this model and 
that would need to be taken into account in order to produce a more realistic model?  For example, might some of the fluxes in the 
cycle depend on other things like the global temperature, which we have not tried to address?   (You may wish to look up information 
about this.) 
\end{enumerate}
\end{enumerate}
\subsection*{Assignment}
Write up the results of your investigations including the spreadsheet results and your answers to the questions above in
a separate document.
\index{carbon cycle|)}

\begin{xcb}{Chapter Exercises}
\vskip 5pt
\begin{enumerate}
\item In Example~\ref{FirstPopExSols}, we noticed that the model 
\begin{align*}
J(n+1) &= 0.05\cdot J(n) + 15\cdot A(n)\\
A(n+1) &= 0.05\cdot J(n) - 0.7895\cdot A(n)
\end{align*}
appeared to have a solution settling down to an equilibrium value when we used the initial conditions $J(0) = 0$, $A(0) = 30$.
(See Figures~\ref{fig:EquilJ} and \ref{fig:EquilA}.)  \vskip 5pt
\begin{enumerate}
\item[(a)]  Does the solution tend to an equilibrium if we change the initial conditions?  If so, do the equilibrium values of $J,A$
depend on what the initial conditions are?  If not, what happens?\vskip 5pt
\item[(b)]  Now, keeping the initial conditions as before, investigate the effect of changing $b = $ the coefficient of $A(n)$
in the first equation in the system.  What happens if $b > 15$?  What happens if $b < 15$?  Does that seem reasonable
on biological grounds?  Explain.
\index{equilibrium!solution}
\end{enumerate}

\vskip 5pt
\item In Example~\ref{MoreGeneralStructured}, we developed a model for the dynamics of a human population stratified into 7 groups by age.\vskip 5pt
\begin{enumerate}
\item[(a)]  What would  be reasonable values for the constants $b_i$?  Recall these are birth rates per unit time for mothers
in each of the age strata, but the time unit is 10 years.  What would be reasonable values for the constants $d_i$?  Recall these are death rates
per unit time.\vskip 5pt
\item[(b)]  Using Excel or other appropriate software, investigate the behavior of the solutions of the model \eqref{eq:GenStructModel}  if 
$$b_1 = 0.001, b_2 = 0.08, b_3 = 0.2, b_4 = 0.15, b_5 = 0.1, b_6 = 0.06, b_7 = 0.0001$$
and
$$d_1 = 0.003, d_2 = 0.004, d_3 = 0.006, d_4 = 0.008, d_5 = 0.05, d_6 = 0.1, d_7 = 0.25.$$
Track the solutions to $n = 20$ (at least) and describe the behavior qualitatively.  In particular, is the total population 
undergoing growth, decline, or does it seem to be reaching an equilibrium?  \vskip 5pt
\item[(c)]  Repeat part (b), but changing the value of $b_3$ to $0.8$.  (Is that value unrealistically high?  Think about
the fact that the time step is 10 years(!))
\end{enumerate}

\vskip 5pt
\item Convert the information from the compartment diagram in Figure~\ref{fig:CompartmentEx} into a system of difference
equations.   For each compartment, assume that any portion of that compartment that does not flow out remains there
in the next time step.  For example, this means that the difference equation for $A$ should have the form
$$A(n+1) = 0.4\cdot A(n) + \cdots .$$
\setlength{\unitlength}{.6mm}
\begin{figure}
\begin{picture}(180,45)
\put(10,5){\line(1,0){40}}
\put(50,5){\line(0,1){20}}
\put(10,5){\line(0,1){20}}
\put(10,25){\line(1,0){40}}
\put(70,5){\line(1,0){40}}
\put(110,5){\line(0,1){20}}
\put(70,5){\line(0,1){20}}
\put(70,25){\line(1,0){40}}
\put(130,5){\line(1,0){40}}
\put(170,5){\line(0,1){20}}
\put(130,5){\line(0,1){20}}
\put(130,25){\line(1,0){40}}
\put(20,5){\vector(0,-1){5}}
\put(80,5){\vector(0,-1){5}}
\put(140,5){\vector(0,-1){5}}
\put(50,16){\vector(1,0){20}}
\put(70,12){\vector(-1,0){20}}
\put(110,15){\vector(1,0){20}}
\put(27,14){$A$}
\put(87,14){$B$}
\put(147,14){$C$}
\put(22,-1){$.2 \cdot A$}
\put(82,-1){$.15 \cdot B$}
\put(142,-){$.17 \cdot C$}
\put(54,19){$.4 \cdot A$}
\put(52,6){$.25 \cdot B$}
\put(112,18){$.2 \cdot B$}
\put(88,37){$.5\cdot C$}
\put(90,29){\oval(120,8)[t]}
\put(30,29){\vector(0,-1){4}}
\put(150,29){\line(0,-1){4}}
\end{picture}
\caption{A three-compartment model.}
\label{fig:CompartmentEx}
\end{figure}

\item This problem is adapted from information from the U.S. National Park Service.\footnote{From {\tt https://www.nps.gov/articles/denali-caribou-herd.htm},
\index{Denali caribou herd}
accessed July 19, 2017.}  Bulls (males) in the Denali caribou herd in Alaska can be divided into different age cohorts with different 
survival rates.  Assume there are 1500 cows (females) in the herd in addition to the bulls\footnote{Of course the real number of 
cows is also changing from year to year.  But for the purposes of this problem we will make a simplifying assumption.} and 
each of them produces one calf each year.  But only $50\%$ of them are males and only $20\%$ of them survive as long as one year because
the very young male calves are vulnerable to wolves, eagles, and other predators.  Young bulls between 1 and 4 years of age have
grown to a size where they can protect themselves much more succesfully and have an annual survival rate about about $90\%$.  They are not generally sexually active.
Bulls between 5 and 7 years old are fully mature.  For a period of 3 to 4 weeks each fall, during
the annual ``rut,'' they do not eat and spend much of their energy on an intensely competitive combat with other bulls.  If they are
successful, they may spend additional energy on mating with cows.  As a result, they can 
lose one-third of their normal body mass each year and they spend the rest of the time after the rut season each year getting back
to that normal mass to repeat the cycle.  Because of this annual stress, a bull's annual survival rate is reduced to 
about $78\%$ during this phase of their lives.  After this, in years 8 and 9, survival declines steeply to about $15\%$ per year.  
For the purposes of this problem, we will assume no bulls survive past age 10. \vskip 5pt
\begin{enumerate}
\item[(a)]  Make a diagram with boxes representing the populations of bulls each age from 1 through 10 years.  Put in arrows
representing the additions and subtractions from each of those population groups each year.  \vskip 5pt
\item[(b)]  Convert your diagram from part (a) into a system of difference equations.  Unlike most of the other models
in this section, one of the difference equations in this system will \emph{not be homogeneous}.  Which one is it and why?  \vskip 5pt
\item[(c)]  Using Excel or other software, investigate the solutions of you system from part (b).  If everything remained constant
so that this model continued to hold, would the population of bull caribou be growing, declining, or tending to an equilbrium?
\end{enumerate}


\vskip 5pt
\item (For readers of \S 3 of Chapter 7.)  Show that a homogeneous 
linear 2nd order difference equation with constant coefficients, $Q(n+2) + r\cdot Q(n+1) + s\cdot Q(n) = 0$, is equivalent to
the system of (first order) difference equations
\begin{align*}
P_1(n+1) &= P_2(n)\\
P_2(n+1) &= -r\cdot P_2(n) - s\cdot P_1(n),
\end{align*}
if we let $P_1(n) = Q(n)$ and $P_2(n) = Q(n+1)$.\footnote{This idea can be developed to show that all
second order and higher order 
difference equations can also be studied by means of the techniques introduced in this chapter. }

\vskip 5pt
\item Explain how the system of difference equations for the SIR model given in \eqref{eq:SIRModel} is consistent
with the assumption that the population is constant.  (Hint:  What is $S(n+1) + I(n+1) + R(n+1)$?)

\vskip 5pt
\item
\label{ex:SIRExampleVerify}
In this exercise you will verify the claims we made about the solution of the SIR model with $\beta = 0.004$ and $\gamma = 0.11$
studied in Examples~\ref{SIRExample} and \ref{SIRExampleComment}. \vskip 5pt
\begin{enumerate}
\item[(a)]  Using Excel or other software, generate the solution of the SIR model with these parameters and the 
initial conditions $S(0) = 99$, $I(0) = 1$, $R(0) = 0$.  Generate plots like the ones from Figure~\ref{fig:SIRPlot} and compare.  \vskip 5pt
\item[(b)]  Examine the values $S(20)$ and $S(21)$ you computed in part (a).  Explain why the number of infected
individuals begins to decrease at $n = 21$.  \vskip 5pt
\item[(c)]  Repeat part (a) but changing the initial conditions to $S(0) = 20$, $I(0) = 1$, $R(0) = 0$.  How
is this solution different?
\end{enumerate}

\vskip 5pt
\item
\label{ex:SIRRecoveredDead}  
\index{model!SIR}
\index{SIR model}
In this exercise, you will develop a modified SIR model that tracks individuals who recover separately from those who die.  \vskip 5pt
\begin{enumerate}
\item[(a)]  The basic idea here is quite simple -- we just want to split the $R$ box in Figure~\ref{fig:SIRDiagram} into two boxes, with
one $R_1$ box representing the \emph{recovered} individuals, and one $R_2$ box representing the individuals who \emph{die} as a result of the 
infection.  Develop an appropriate compartment diagram along these lines.  In most cases the rate $\gamma_1$ at which individuals leave the $I$ group
and enter the $R_1$ group (recovered individuals) will be different from the rate $\gamma_2$ at which individuals leave the $I$ group and enter the $R_2$ 
group of individuals who die.  Your diagram should allow for that. \vskip 5pt
\item[(b)]  Translate your diagram into a system of \emph{four} difference equations for the functions $S(n), I(n), R_1(n), R_2(n)$. \vskip 5pt
\item[(c)]   The first two equations in your model should be the same as those in \eqref{eq:SIRModel}, except the $- \gamma\cdot I(n)$
term in the second is replaced by something like $ - (\gamma_1 + \gamma_2)\cdot I(n)$.  Take the values $\beta = 0.04$ and 
$\gamma_1 = 0.09$, $\gamma_2 = 0.02$, $S(0) = 99$, $I(0) = 1$, $R_1(0) = 0$ and $R_2(0) = 0$.  Using Excel or other software, 
generate the solution of the modified SIR model with these parameters.  Generate plots like the ones from Figure~\ref{fig:SIRPlot} and
describe the course of the infectious outbreak.  \vskip 5pt
\item[(d)]  Refer  to \eqref{eq:BasicReproductionNumber}.  What is the threshold value for epidemics with this model?
\end{enumerate}

\vskip 5pt
\item
\label{ex:SIRAge}  
In this exercise, you will develop a simple age-structured SIR model.  Suppose for simplicity that we have a 
population $P$ split into two strata by age (roughly, think ``young'' and ``old'' individuals).  Each age stratum
will contain susceptible, infected, and removed subgroups as in the basic SIR model.  Call these
$S_1,S_2,I_1,I_2,R_1,R_2$  The new feature
is that infected people in either age stratum can infect others within either age stratum.\footnote{If not, if for instance young
people could only infect other young people, we would essentially have two separate SIR models.}\vskip 5pt
\begin{enumerate}
\item[(a)]  Develop an appropriate compartment diagram for this situation.  There will be 6 boxes in all.  To allow
for differences in resistance to the disease, etc. allow for two different constants $\gamma_1$ and $\gamma_2$
governing the arrows from the infected boxes to the removed boxes.  Also use four different constants
$\beta_{11}$, $\beta_{12}$, $\beta_{21}$, and $\beta_{22}$ in the arrows from the susceptible boxes to the infected boxes
because contacts between people in different age strata might produce infections at different rates.  For instance, 
you can think of $\beta_{ij}$ as the constant for infections between susceptibles in stratum $i$ and infecteds in stratum $j$. \vskip 5pt
\item[(b)]  Translate your compartment diagram into a system of 6 difference equations.
For example, the equation for $S_1$, the susceptibles in stratum 1 will look like this:
$$S_1(n+1) = S_1(n) - \beta_{11} \cdot S_1(n)\cdot I_1(n) - \beta_{12}\cdot S_1(n)\cdot I_2(n).$$
\item[(c)]  Using Excel or other software, compute some solutions of these models.  Use $\beta_{ij} = 0.004$ and $\gamma_i = 0.11$
for all $1 \le i,j \le 2$ and $S_1(0) = 99$, $S_2(0) = 100$, $I_1(0) = 1$, $I_2(0) = 0$, and $R_1(0) = R_2(0) = 0$.  \vskip 5pt
\item[(d)]  Change the parameters in part (c) letting $\beta_{21} = .01$ and $\beta_{22} = .007$ instead.  What changes?
\end{enumerate}

\vskip 5pt
\item 
\label{ex:SIRVaccination}
\index{vaccination}
With the basic SIR model, or modifications thereof, there are several ways we can study the effects of vaccination for an 
infectious viral disease.  Imagine the situation of an isolated human population of size $P$.\vskip 5pt
\begin{enumerate}
\item[(a)]  If we have the luxury of immunizing before any outbreaks occur, but we know the constants $\beta$ and $\gamma$
characteristic of the disease, how should public health authorities plan a vaccination campaign that will create enough
``herd immunity'' so that no epidemic outbreaks of the disease are possible?  How many $V$ out of the total population $P$
need to be vaccinated?  Not everyone needs to be vaccinated; how
many vaccinations are necessary? \vskip 5pt
\item[(b)]  Once an outbreak is underway, it is still possible to attempt to slow it down or drive it to conclusion via a vaccination 
campaign.  One approach would be to vaccinate some fixed number $v_0$ of the susceptible group each time step.   How
could the basic SIR model be modified to take this into account?  Construct a new modified compartment diagram representing
this situation and translate it into a system of difference equations.\footnote{Eventually the number of susceptible individuals will 
be reduced to negative values with this method.  When that happens, the solutions of the model stop being realistic.}
 (Hint:  This is just adding one more arrow to the compartment diagram.)  \vskip 5pt
 \item[(c)]  Another approach besides part (b) would be to vaccinate some fraction of the susceptible population each
 time step.  How could the basic SIR model be modified to take this into account?  Construct a new modified compartment diagram representing
this situation and translate it into a system of difference equations.
\end{enumerate}

\vskip 5pt
\item 
\index{equilibrium!solution}
\label{ex:LVLogEquil}  Show that the modified predator-prey system \eqref{eq:LVLog} has equilibria at $(0,0)$, $(0,M)$ from
the logistic model for the prey when no predators are present,  \emph{and} at 
$$(P_0,Q_0) = \left( \left(1 - \frac{a}{b\cdot M}\right)\cdot\frac{c}{d}, \frac{a}{b}\right)$$
in the interior of the first quadrant.

\vskip 5pt
\item
Refer to Example~\ref{FirstLVLog}. \vskip 5pt
\begin{enumerate}
\item[(a)]  Verify the calculations presented in that example using Excel or other suitable software.
\item[(b)]  What happens to the computed solutions if the value of the coefficient $a$ is changed to 
$a = .8$ (meaning that the predators die out faster if no prey are present)?
\item[(c)]  What happens to the computed solutions if $a = .1$?
\end{enumerate}

\end{enumerate}
\end{xcb}
\index{model!compartment|)}
\index{compartment model|)}

\endinput

