% SIMIODESample.tex and simiode.cls for LaTeX2e
% Version 11 March 2013.
% 
% This class (simiode.cls) is a modification of that used by 
% the journal PRIMUS - Problems, Resources, and Issues in 
% Mathematics Undergraduate Studies. 
% 
% See http://www.tandfonline.com/loi/upri .
% 
% The class was developed by  Marc Frantz with the  suggestions of the 
% PRIMUS Editorial Board, especially John Thoo, and is much appreciated.
%
% This template "SIMIODEPaper.tex" is to be used with the 
% class file "simiode.cls".  You may rename the template.  
% The latest version of this class and template files 
% can be found at
% 
%                <http://www.SIMIODE.org> .
%
% For an explanation of using the class file "simiode.cls", 
% see the pdf ``Description and Illustration of Use: 
% A LaTeX Class File for SIMIODE'' found at  
% 
%                <http://www.SIMIODE.org> . 
%
% A "Quick Start Guide" appears below.
%
% For information on using TeX or LaTeX, go to %                 <http://www.tug.org>.

% This file must be processed with version LaTeX2e or higher.
%
\documentclass{simiode}

\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}

\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{hyperref}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding double spacing %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manuscripts should be submitted doubled spaced.  This is 
% accomplished via the \baselinestretch command below.
% However, you may obtain an idea of how your article will 
% look in final form, so that you can better plan your 
% presentation.  To see your manuscript formatted very 
% similarly to the on-line published format of  SIMIODE, 
% comment the following line out by inserting a percent sign % (%) at its beginning.

\renewcommand{\baselinestretch}{1.3}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding macros %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.


\newcommand{\AmSLaTeX}{$\cal A$\kern-.1667em\lower.5ex\hbox{$\cal
M$}\kern-.125em $\cal S$-\LaTeX}
\newcommand{\BibTeX}{{\sc Bib}\TeX}
\newcommand{\scls}{{\tt simiode.cls}}
\newcommand{\stemp}{{\tt SIMIODEPaper.tex}}
\newcommand{\SI}{SIMIODE}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Top matter is above and 
% the beginning of the document below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\newpage

\title{\ifnum\TeacherEdition=1 TEACHER VERSION \else STUDENT VERSION \fi\\Disease Spread}




\markboth{Modeling the Spread of Disease}{Modeling the Spread of Disease}





\author{Brian Winkel, Director SIMIODE\\
Chardon OH 44024 USA\\
BrianWinkel@simiode.org}

\makeStitlePDFLaTex


\begin{abstract}
We offer a physical situation, using a grid and M\&M candies, to simulate the spread of disease. Students conduct the simulation and collect the data which is used to estimate parameters (in several ways)  in a differential equation model for the spread of the disease. Students  produce the standard logistic equation when trying to model the interaction between those who have the disease and those who do not. Reasoned assumptions need to be articulated by students.
\end{abstract}


\section*{SCENARIO DESCRIPTION}


We offer a modeling opportunity in which the phenomenon of the spread of disease  can be described by one differential equation. 

First, let us consider growth in which we have a quantity which grows by a constant amount, say, $k$, i.e.,
\begin{equation}\label{eq:straight}
y'(t) = k\,\, , \,\, y(0) = y_0 \,\, 
\end{equation}

\noindent In this case we can integrate both sides to obtain $y(t) = k*t + c$ where the initial condition $y(0) = y_0$ determines $c$ so that we have a complete solution, $y(t) = k*t + y(0)$.

A more interesting growth equation comes from assuming that the growth rate of some quantity is proportional to  the amount present, e.g., money earns interest and populations breed. Here


\begin{equation}\label{eq:expGrowth}
y'(t) = k \, y(t)\,\, , \,\, y(0) = y_0 \,\, 
\end{equation}

\section*{Exponential Growth Population Modeling\label{sec:exp}}


Exponential growth (and its related decay issues) stem from the following differential equation:

\begin{equation}
y'(t) = k y(t)\,\, , \,\, y(0) = y_0 \,\, .\label{eq:ExpGrowth}
\end{equation}

\noindent $k$ is referred to the {\it rate constant\/} and if $k>0$ we call it a growth constant, while if $k<0$ we call it a decay constant.  

This exponential growth model can be used to predict the growth of money where $k$ is the decimal value of the percent interest on continuously compounded interest situations or the growth of organisms where $k$ is the per unit time growth rate.  In terms of parameter estimation the problem of estimating $k$ is universally the same and the various approaches (1) linearizing the solution by taking natural logarithms of the non-time variable and fitting the linearized data using linear regression, (2) estimating growth rates and fitting a linear model relating $y'(t)$ and $y(t)$, or (3) direct nonlinear fitting of the analytic solution to (\ref{eq:ExpGrowth}) to the data, will serve well to estimate $k$. We seek more appropriate models for growth in biological situations and this demands we do a better job of modeling the reality of finite resources or populations. This leads to a limited growth population model. 






\section*{Limited Growth\label{sec:logistic} Population  Model}

\begin{table}[h]
\centering

\begin{tabular}{|c|c||c|c|} \hline
\multicolumn{1}{|c|}{ Time (days)}&\multicolumn{1}{|c||}{ Population (\#)}& \multicolumn{1}{|c|}{ Time (days)}&\multicolumn{1}{|c|}{ Population (\#)}\\ \hline \hline
 0 & 200 & 26 & 783\\ \hline
 2 & 240&28 & 792\\ \hline
 4 & 281&  30 & 837\\ \hline
 6 & 317 &  32 & 800\\ \hline
 8 & 362 & 34 & 866\\ \hline
 10 & 392 & 36 & 933\\ \hline
 12 & 437& 38 & 942\\ \hline
 14 & 496& 40 & 941\\ \hline
 16 & 542&  42 & 986\\ \hline
 18 & 601& 44 & 964\\ \hline
 20 & 615& 46 & 940\\ \hline
 22 & 662 & 48 & 1009\\ \hline
 24 & 758 &  50 & 984\\ \hline
\end{tabular}


\parbox{7.4 true cm}{\caption{\label{table:WinkelPopData} Population data with units of number (\#) of members of the population at time $t$ in days.}}
\end{table}

In Table~\ref{table:WinkelPopData} we have data from a population  which does not actually grow exponentially.  There is a limit to growth. Often this is a natural limit, such as regional resource limitations, geometry of the environment, or borders. Just consider growing microorganisms in a Petrie dish which has only a finite space.





\begin{figure}[h]
\centering
\includegraphics[width=11cm]{1-17-ModelingTwoLogDataGenPlot}

%\vspace*{-1.7 true cm}

\parbox{10.6 true cm}{\caption{\label{fig:1-17-ModelingTwoLogDataGenPlot} The plot of the population data in Table~\ref{table:WinkelPopData}.}}

\end{figure}

\subsubsection*{Conjecturing a growth model for Table~\ref{table:WinkelPopData} data}

\begin{itemize}
\item [1)] Take time to consider how we could modify the differential equation (\ref{eq:ExpGrowth}) for exponential growth so as to actually model growth that would level off or whose growth rate would eventually get smaller and smaller. Discuss your models with colleagues and decide if any of the models under discussion are suitable. Look for characteristics of the differential equation you offer. For example, does your differential equation appear to come to a point where growth is 0 (leveling off)? Does it have any flaws?  For example, does it create some population out of nothing?
\item [2)]Produce some conjectures and  ``talk'' about other models being proposed by your colleagues. Probe each model, test it under obvious conditions, e.g., ``If we have no population do we get no growth in our model?''
\end{itemize}

While not ignoring your efforts we are going to turn - as you turn to the next page - to a well-established model for a while and we will offer up a chance for you to try out your model shortly. So hold on to your ideas. See the last question in this paper in which you can take your model as far as you wish.




The differential equation (\ref{eq:logistic}) is referred to as the {\it logistic differential equation\/}. It is used in such areas as modeling limited growth population biology as well as the study of the spread of rumors, disease, and technologies:
\vspace*{.15 true cm}

\begin{equation}\label{eq:logistic}
\displaystyle{y'(t) = r y(t)\frac{K-y(t)}{K},\qquad y(0) = y_0}\, . 
\end{equation}

\vspace*{.15 true cm}
\noindent Here $y(t)$ is the population in appropriate units. 

\subsubsection*{Scrutinize the logistic differential equation model}

\begin{itemize}
\item [3)] Take some time to give the logistic model the scrutiny you and your colleagues gave to each other's models. Try to push it around a bit, play some ``what if" games with it. Test  it, interpret it, challenge it.  Well, what do you think about this model?
\end{itemize}



In different contexts the logistic differential equation  (\ref{eq:logistic}) could be the percent of the population that has heard a specific  rumor or acquired knowledge or uses  a specific  technology. Or it could be the number in the population that have had a given disease. $r$ represents the intrinsic growth rate (in population models $r$ is birth rate less death rate) with $K$  as the limiting value of $y(t)$.  $K$ is  called the {\it carrying capacity\/}. The latter comes from ecology and refers to the number of the given species that the environment can support or carry.
For  excellent histories of the logistic curve in the development of population biology see \cite{Kingsland1982,Hutchinson1978}.

Now, let us return to the data as offered in Table~\ref{table:WinkelPopData}. This particular data is what we call a {\it toy\/} set of data and is not from any real experiment. In later work with the logistic model we shall use actual data from \cite{Gause1934}. Indeed, in this activity you will generate data from a physical simulation using M\&M candies as members of your population.


Upon plotting the data (see Figure~\ref{fig:ModelingTwoLogDataGenPlot}) and developing a model, presumably the logistic differential equation (\ref{eq:logistic}), based on the  sigmoidal or ``S" shape of the solution of the logistic equation  (\ref{eq:logistic}), we ask  how we might confirm or deny that (\ref{eq:logistic}) is a good model for our data.  Throughout this early discussion we may never actually solve  (\ref{eq:logistic}) by hand - although one could do so by using the separation of variables and partial fraction strategies described elsewhere. Rather, when we do need an analytic solution  we rely on {\it Mathematica\/}'s symbolic solving command {\tt DSolve} or we may go to its numerical solving command {\tt NDSolve} as appropriate. 

In assessing how good a model we may have developed in (\ref{eq:logistic}) we might consider determining criteria for estimating the parameters, i.e., finding the parameters which when we substitute them in  (\ref{eq:logistic}) and solve it will permit us to compare how our model predicts the data in Table~\ref{table:WinkelPopData}.  

\begin{itemize}
\item [4)] Given the data in  Table~\ref{table:WinkelPopData} and as presented in Figure~\ref{fig:ModelingTwoLogDataGenPlot}  how would you confirm or deny that (\ref{eq:logistic}) is a good model for our data, for the phenomena of limited growth.
\item [5)] Before going on, consider how you would estimate the parameters $r$ and $K$ and offer up your ideas and approaches. Discuss the pros and cons  of your approaches and your colleagues' approaches. Write them down. State the steps you might use. Do it here, do it now.
\end{itemize}








\subsection*{Developing strategies for estimating parameters\label{subsection:Strategies}}

Here are our several approaches that usually crop up.  How do your approaches compare?
\begin{enumerate}
 
\item Can we estimate the parameters by transforming the logistic differential equation  model and then determining the  least sum of square differences between transformed data and  model to determine which parameters give the best fit - nonlinear fit of $\frac{dP}{dt}$ vs. $r P\frac{K-P}{K} = a P - b P^2$, where $a = r$ and $b = \frac{r}{K}$? 
\item Can we estimate the parameters by transforming the ODE model and then determining the  least square differences between transformed data and  model to determine which parameters give the best fit - linear fit of $\frac{1}{P}\frac{dP}{dt}$ vs. $r \frac{K-P}{K} = a  - b P$, where $a = r$ and $b = \frac{r}{K}$? 
\item  Can we estimate the parameters by analytically  solving the logistic differential equation model in (\ref{eq:logistic}) and then determining the  least sum of square differences between data and analytic model to determine which parameters which yield   the least sum of square differences, i.e. nonlinear fit of the analytic solution to the data at the time observation points of the data?
\item  Can we estimate the parameters by numerically solving the logistic differential equation model for a wide variety of parameters and determining the least sum of square differences between data and model values for EACH set of parameters and then determine the parameters $r$ and $K$ which give rise to the smallest least sum of square differences? 
\end{enumerate}
 
Approach (4) is a bit strange and rather computationally intense, but it is still doable, namely try lots, and we mean may LOTS, of parameter sets $r$ and $K$ in a reasonable vicinity of what you believe might work to see if you can get systematically lucky and get pretty good estimate. 

However, let us concentrate on the first two approaches (1) and (2).   Even inside approaches  (1) and (2) we have several approaches to transforming the data. For example in each case we can take a first difference approach to estimating the derivative from the data or a symmetric first difference approach to estimating the derivative. Thus, we shall have approaches (1a) and (1b) as well as (2a) and (2b) with the (a) denoting the plain first difference approach and (b) denoting the symmetric first difference approach. 


In all approaches we shall use the same measure of how good the method is, namely the sum of square errors, i.e., the squares of the difference between the value of the  solution to the logistic differential equation~(\ref{eq:logistic}) evaluated at the time value for each data point and the corresponding population value at that time. This would look like: $$\sum_{i=1}^n (\hat{y}(t_i,r,K) - P_i)^2 \, , $$ where $\hat{y}(t_i,r,K)$ is the value of our  solution to our logistics differential equation~(\ref{eq:logistic}) with parameters $(r,K)$ at time $t_i$ and $P_i$ is the value of the population at time $t_i$ from our data set for $i=1, \ldots, n$, while  $n$ is the size of our data set.


 
\section*{Parameter estimation opportunities}

\begin{itemize}
\item[6)] Use approach (1a) and (1b) to estimate parameters $r$ and $K$ using Trendline in EXCEL and fitting a parabola with intercept zero to the transformed data. Use your parameters in each case and return to {\it Mathematica\/}. Enter the data as a list of pairs of numbers and {\tt ListPlot} that data, calling the plot, dataPlot. Now solve the differential equation and grab the solution with this line of {\it Mathematica\/} code:

{\tt ys[t\_] = y[t] /.DSolve[{y'[t] == r y[t] (K - y[t])/K, y[0] == y0}, 

\qquad \qquad y[t], t][[1]]}

\noindent with the parameters from your EXCEL adventure. Your solution should look like this, only with your numerical estimates for $r$ and $K$ from your EXCEL adventure, instead of a general from equation with the letters $r$ and $K$ as in (\ref{eq:logSolved}):  

\begin{equation}\label{eq:logSolved}
ys(t) = \frac{K y_0 e^{r t}}{K+y_0 e^{r t}-y_0}\, .
\end{equation}


Now after substituting $y_0 = 200$ (see Table~\ref{table:WinkelPopData}) into
(\ref{eq:logSolved}) we proceed to evaluate the solution at each time $t = 0, 2, 4, \ldots , 50$ from the data; subtract the observed corresponding data values, 200, 240, 281, \ldots, 984; square these differences; and add the square differences or errors to obtain: 
\begin{equation}\label{eq:ss}
ss(r,K) = \sum_{i=1}^n (\hat{y}(t_i,r,K) - P_i)^2 \, . 
\end{equation}

We seek the values of $r$ and $K$ that will minimize the sum in (\ref{eq:ss}) and so we call up {\it Mathematica\/}'s {\tt FindMinimum} command. From this  output we can read our best $r$ and $K$ values. Take these and substitute them into (\ref{eq:logSolved}) along with $y_0 = 200$ and we have a fully developed model which we can plot (call it ModelPlot) and  compare with our data plot by {\tt Show}ing them together in {\it Mathematica\/}. 
 
Discuss your results from (1a) and (1b) approaches.


\item[7)] Use approaches (2a) and (2b) to estimate parameters $r$ and $K$ using Trendline in EXCEL and fitting a linear function  to the transformed data. Use your parameters in each case and return to {\it Mathematica\/}. Enter the data as a list of pairs of numbers and {\tt ListPlot} that data, calling the plot, dataPlot. Then solve the differential equation and grab the solution with this line of {\it Mathematica\/} code:

{\tt ys[t\_] = y[t] /.DSolve[{y'[t] == r y[t] (K - y[t])/K, y[0] == y0}, 

\qquad \qquad y[t], t][[1]]}

\noindent with the parameters from your EXCEL adventure. Your solution should look like this:

\begin{equation}\label{eq:logSolved}
ys(t) = \frac{K y_0 e^{r t}}{K+y_0 e^{r t}-y_0}\, .
\end{equation}


After substituting $y_0 = 200$ (see Table~\ref{table:WinkelPopData}) into
(\ref{eq:logSolved}) we proceed to evaluate the solution at each time $t = 0, 2, 4, \ldots , 50$ from the data; subtract the observed corresponding data values, 200, 240, 281, \ldots, 984; square these differences; and add the square differences or errors to obtain: 

\begin{equation}\label{eq:ss}
ss(r,K) = \sum_{i=1}^n (\hat{y}(t_i,r,K) - P_i)^2 \, . 
\end{equation}

We seek the values of $r$ and $K$ that will minimize the sum in (\ref{eq:ss}) and so we call up {\it Mathematica\/}'s {\tt FindMinimum} command. From this  output we can read our best $r$ and $K$ values. Take these and substitute them into (\ref{eq:logSolved}) along with $y_0 = 200$ and we have a fully developed model which we can plot (call it ModelPlot) and  compare to our data plot by {\tt Show}ing them together in {\it Mathematica\/}. 

Discuss your results from (2a) and (2b) approaches.

\item[8)] Finally, directly estimate $r$ and $K$ by first evaluating the solution, (\ref{eq:logSolved}), to the logistic differential equation~(\ref{eq:logistic}) at each time value $t = 0, 2, 4, \ldots , 50$ along with $y_0 = 200$ and form the sum of square errors, as before. Except that now we can simply {\tt FindMinimum} this sum of square errors to find directly which $r$ and $K$  make the solution ((\ref{eq:logSolved})) with these $r$ and $K$ values our best fitting model.  

Take these values of $r$ and $K$ and substitute them into (\ref{eq:logSolved}) along with $y_0 = 200$ and we have a fully developed model which we can plot (call it ModelPlot) and  compare to our data plot by {\tt Show}ing them together in {\it Mathematica\/}. 

Discuss your result  and compare it with those from (1a) and (1b) as well as from (2a) and (2b) approaches.
\end{itemize}

\section*{Modeling spread of disease with M\&M simulations}

While we will examine some real data from the literature and attempt to model it with  a logistic model, we take some time here to develop our own data of the spread of disease using this logistic model approach. Namely, if in a population of size $K$ we have $y(t)$ members who have a communicable disease and can spread it upon contact with susceptible members of which there are $(K - y(t))$ then a possible model for the rate of the spread of the disease could be the logistic equation:

\begin{equation}\label{eq:spread}
y'(t) = \alpha y(t) (K-y(t))
\end{equation}

\noindent where the product $y(t)*(K-y(t))$ represents the number of interactions that could take place per unit time and $\alpha$ is a constant which reflects just how many interactions take place AND haw many result in the spread of the disease. Of course we could convert to percentage of the population whereupon $K = 1$ stands for 100\% of the population, $y(t)$ stands for the percentage of the population who have the disease, and   $ K-y(t) = 1 - y(t)$ stands for the population which does not have the disease and is susceptible.

We recognize (\ref{eq:spread}) as an over simplified model of the spread of disease, because it does not reflect several realities, among them the need to accommodate immune individuals, the various stages of gestation of the disease before it might be transmittable, effects on various age groups, etc. Nevertheless, we shall proceed to work with the model of (\ref{eq:spread}), picking up on more realistic variations as our modeling and differential equations skills increase!

\begin{itemize}
\item [9)] Articulate the assumptions and oversimplifcations we need to make to postulate (\ref{eq:spread}) as a model  for the spread of disease.
\end{itemize}


So let us simulate a disease spread among M\&M's!!!  We shall need a regular size bag of M\&M chocolates - not the peanut kind (there should be about 55 pieces in this size bag); a small  cup or glass to hold the M\&M's; and a grid (copy in Figure~\ref{fig:ModelingTwoMMGridForLogisticGrowth} blown up to a size (300\%) so that each grid is 1.32 cm high and wide - the radius of an M\&M is 1.32 cm. Indeed, in the folder for this scenario there is a pdf file, M\&MLayout.pdf which, if printed, is the exact size of your population surface which you will need. Otherwise, follow instructions in caption of Figure~\ref{fig:ModelingTwoMMGridForLogisticGrowth}.

\begin{figure}[h]
\centering
\includegraphics[width=7cm]{1-17-ModelingTwoMMGridForLogisticGrowth.eps}

%\vspace*{-1.7 true cm}

\parbox{12.2 true cm}{\caption{\label{fig:1-17-ModelingTwoMMGridForLogisticGrowth} A grid used for our simulation is to be reproduced so that each grid is 1.32 cm high and wide when expanded to 300\% - the radius of an M\&M is 1.32 cm. To make the necessary grid sheet for the simulation copy this actual grid on a copy machine at 300\% - making sure it exactly fills an 8.5" x 11" sheet.}}

\end{figure}

\subsubsection*{The Simulation Experiment}

We are going to simulate the spread of a disease. Prepare and then  cut out a grid of squares  each 1.32 cm by 1.32 cm where 1.32 cm is the diameter of an M\&M on an 8.5'' x 11'' sheet of paper.  See instructions in Figure~\ref{fig:ModelingTwoMMGridForLogisticGrowth}. Each grid box represents an individual M\&M.  We  put bold X's in 8 of the boxes to indicate that the disease initially has 8 individual M\&M's infected.  Select these randomly. If you are not confident in doing so then use these {\it Mathematica\/} instructions to select randomly
\vspace*{.15 true cm}


\centerline{\tt Table[\{Random[Integer, \{1, 16\}], Random[Integer, \{1, 21\}]\}, \{i, 1, 8\}]}
\vspace*{.15 true cm}

\noindent where there will   be 16 M\&M's width and 21 M\&M's height.

Now tape one edge of the 8.5'' x 11'' sheet of paper to a desk or table  flat surface and then line the sheet with four vertical walls on each edge, e.g., four reams of  8.5'' x 11'' paper to make a bounded fenced region.

Take the 55 M\&M's in a standard M\&M bag and place them in a cup. These form our healthy population, the susceptibles. Thus we have 63 M\&M's in our population - 8 infected and 55 susceptibles - at the start of our simulation.

\begin{figure}[h]
\centering
\includegraphics[width=7cm]{1-17-ModelingTwoExp1MMDataCollection.eps}

%\vspace*{-1.7 true cm}

\parbox{12.9 true cm}{\caption{\label{fig:1-17-ModelingTwoExp1MMDataCollection} A grid used for our first simulation of infecting M\&M's. Notice the 8 original infectives with an X in each box and then the consecutive infectives numbered 1-55.}}

\end{figure}


Gently toss the M\&M's from the cup ``randomly" on the grid in the bounded fenced region. Each time we toss the M\&M's we will consider a generation or cycle or unit of time. If an M\&M ``touches" (overlaps in any part) a box with an X in it we mark the box in which most of the M\&M resided with a number (starting with 1 for our first infected M\&M beyond the original 8, and then a 2 for our second, and so forth). As an M\&M becomes infected we remove it from the healthy population, noting that our number of infecteds is in addition to the original or current infected population. The newly numbered boxes stand for infected M\&M's who stand ready (and willing!) to infect other susceptible M\&M's should one of these susceptible M\&M's touch one of the infected M\&M squares.  See Figure~\ref{fig:1-17-ModelingTwoExp1MMDataCollection} for a summary simulation.
 

Figure~\ref{fig:MMCount} (in Black and White) shows the result of the simulation at some point in time along the process.

\begin{figure}[h]
\centering
\includegraphics[width=7cm]{1-17-MMCount.eps}

%\vspace*{-1.7 true cm}

\parbox{12.9 true cm}{\caption{\label{fig:1-17-MMCount} Photograph of actual grid and M\&M's at some point in time of the simulation.  Can you identify the newly infected M\&M's from this picture?}}

\end{figure}











At each generation we toss the cup, determine those M\&M's who get infected, increment the number of newly infecteds, and write that number in a square nearest to  the M\&M that just got infected. We also keep a tally sheet of the total number of infected at each new generation or cycle or time interval. We show such numerical results in Table~\ref{table:Exp1}.

\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|} \hline \hline
 Generation &  \# Infected \\ \hline \hline
 0 & 8 \\ \hline
 1 & 11 \\ \hline
 2 & 16 \\ \hline
 3 & 22 \\ \hline
 4 & 28 \\ \hline
 5 & 37 \\ \hline
 6 & 46 \\ \hline
 7 & 48 \\ \hline
 8 & 54 \\ \hline
 9 & 58 \\ \hline
 10 & 61 \\ \hline
 11 & 62 \\ \hline
 12 & 62 \\ \hline
 13 & 62 \\ \hline
 14 & 62 \\ \hline
 15 & 63\\ \hline
\end{tabular}

\parbox{12.2 true cm}{\caption{\label{table:Exp1} Data from our first experiment in simulating infection of M\&M's. Yours should look somewhat the same, don't you think?}}
\end{center}

\end{table}


The ``touching" or overlap simulate contact and we presume contact is sufficient to transmit the disease, i.e. absolute contagion. Of course, we could toss a die, for example, to determine for each touch what the chance of spread of the disease could be, using some preassigned rule such as if face of die shows 5 or 6 then there is no infection in this case; doing this for each touched situation. This would slow up the process, but make it more realistic in some sense as most diseases are not automatically transmitted upon touch.

\begin{itemize}
\item [10)] Either from the data in Table~\ref{table:Exp1} or from one of your own experiments analyze the data by assuming the logistic differential equation~(\ref{eq:logistic1}) could be used to  describe the spread of the disease between the susceptible ($y(t)$) and infected ($K - y(t)$) M\&M's. 

\begin{equation}
\displaystyle{y'(t) = r y(t)\frac{K-y(t)}{K},\qquad y(0) = y_0}\, . \label{eq:logistic1}
\end{equation}

Then use some of our previously discussed methods to determine the parameters $r$ and $K$.
\item [11)] Discuss the fact that if you do not assume $K = 63$, i.e.~there are eventually at most $8 + 55 = 63$ M\&M's in our population (the original 8 infectives and the 55 susceptibles from your bag of yummy M\&M's) you might get a different model than if you do NOT assume a value for $K$.
\item [12)] Discuss (and do!) variations of the model, such as 
\begin{itemize} 
\item making the squares on your grid larger or smaller;
\item increasing the number of original infectives;
\item do not put the new infectives where they were infected but rather randomly put them in some other position (they could have wandered a while before dying);    
\item making a grand sheet of taping four 8.5" x 11" sheets to form a large 17" x 22" sheet with the same number of initial infecteds and the same number of susceptibles; or
\item  making your own changes.
\end{itemize}
\item [13)] What is the value of $r$ that you get in each case?  Does it change in any way? What does $r$ mean? What is the significance of $r$ without the limiting value $K$ (i.e.  unlimited bag of M\&M's - yummy, yummy!!!)
\item [14)] A time intensive (possibly project) activity is to use some random feature to determine that IF an uninfected or susceptible M\%M comes in contact with an infected one then with some probability the disease is passed on -- not all contact results in the spread of the disease in this case. For example, use a six sided die and  consider each contact, one by one. Toss the die, and if, say, it is a 3 or a 5 (you make up your own rule), then an infection actually takes place; otherwise no infection takes place and the potentially infected M\&M goes back into the susceptible pool. Change your probability rule.  Does increasing the chance of infection at each step do anything to the data, to the plot, to $r$? 
\item [15)] Take your model or one of those you develped before we ``spilled the beans'' (M\&M's actually!) and shared the logistic differential equation and address questions (3) - (7) above with that model, addressing the issues as best you can.
\end{itemize}

\teacher{
\section*{NOTES FOR TEACHERS}


We have never done all of these activities in one ``setting.'' Rather, we have used bits and pieces, broad swatches of this material in various situations and students are very good at addressing how we might estimate parameters in the logistic differential equation, probably because they are familiar with EXCEL and differencing and applying functions (e.g., logarithm) to whole columns of numbers is quite familiar to them.

When it comes to building models this is an early activity to try to build in limited growth and the students really struggle with this. Often, I give them lots of time, say 30 minutes with some individual and small group interaction with me, as well as total class time in which we examine some candidates and teach how to test at some obvious points, e.g., when the population is 0 does the model in question say there is growth - not good! Most times, indeed, almost always the class will come around to some form of the logistic differential equation and feel pretty confident about the group effort. This is good, for when we actually turn the page and offer them the same model and tell them that historically many mathematicians and ecologists have come to the same conclusion they are pleased with themselves.

This is an early activity in which we go after parameters through estimation techniques and so we take time to compare our results and methods. We depend heavily on technology, Mathematica and EXCEL and we are often teaching students the syntax en route to modeling success, so patience is the order of the day. However, since we do lots of parameter estimation efforts in our course this pays off often.

Instead of offering comments here we address each of the Activities offered in a separate Mathematica notebook, both in .nb and .pdf format and offer detailed comments on  mathematics, technology use, and student efforts.
}

\begin{thebibliography}{99}




\bibitem{Gause1934} Gause, G.~F. 1971. {\it The Struggle for Existence}. New York:  Dover Publications, Inc. First published in 1934 by The Williams \& Wilkins Company. 



\bibitem{Hutchinson1978} Hutchinson, G. Evelyn. {\it Introduction to Population Ecology\/}. New Haven CT:  Yale University Press.

\bibitem{Kingsland1982} Kingsland, Sharon. 1982. The Refractory Model:  The Logistic Curve and the History of Population Ecology. {\it The Quarterly Review of Biology\/}. 57: 29-52.





\end{thebibliography}
\end{document}
