% SIMIODETemplate-TeacherORStudent.tex and simiode.cls for LaTeX2e
% Version 29 July 2015.
% 
% 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 suggestions of the 
% PRIMUS Editorial Board, especially John Thoo.
%
% This template is to be used with the  class file "simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
%
\documentclass{simiode}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Basically the Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review (but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
% This option was created  by Eric Sullivan, Carroll College, Helena MT.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\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.

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

\begin{document}

\newpage
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote. See above for instructions on toggling from Teacher Version to Student Version.
\title{\ifnum\TeacherEdition=1 TEACHER VERSION  \else STUDENT VERSION \fi\\INSECT COLONY SURVIVAL OPTIMIZATION}



\markboth{Insect Colony Optimization}{Insect Colony Optimization}


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

\makeStitlePDFLaTex



\begin{abstract}
 We present a system of nonlinear differential equations to model the  control of energy flow into producing workers or reproducers in an insect colony. With a set of given parameters and a number of different  energy  functions we permit students to explore options as to just how much of the input energy available from nature is diverted from worker production to reproducer production. Finally, we move students to consider a bang-bang solution for our control function and, in doing so, introduce them to  the notions of optimal control. We offer some running commentary which suggests alternate models from a bee hobbyist  referee.
\end{abstract}








\section*{SCENARIO DESCRIPTION}


We model insect colony propagation or survival from nature using differential equations. We ask you to analyze and report on what is going on and  develop an optimal strategy for success of a natural phenomena.

\vspace*{.3 true cm}
Your responses are to contain traditional elements of good reporting, complete statement of the problem - copied from the problem or reformulated in your own words, your analysis and results, and discussion of issues relevant to the analyses and results. 

\vspace*{.3 true cm}


In  considering the life cycle of colonies of insects it is  believed  that 
the colony appears to maximize the number of reproducer insects at the  end 
of  its  season.   This enables the colony to spread  its  genes  and  thus 
propagate  its own kind in an optimal fashion.  But how does the colony  do 
this? One could consider different models for hibernating (e.g., wasp) and non-hibernating (e.g., honey bees) insects.  The non-hibernating case is more complex due to the need to build up substantial resources and hives to survive the winter.

We consider an energy model for the colony in which we discuss two kinds of 
insects in the colony, workers and reproducers (or queens).

Let  $W(t)$ be the number of worker insects in the colony at time  $t$  (usually measuring days into  the season) and 
let  $Q(t$) be the number of reproducer insects or queens in the colony  at 
time $t$. We will not presume any age-class distinctions.

A proposed model \cite{MaceviczOster1976} for rates of growth of each species is (can we defend/explain  this model?):

\begin{eqnarray}\label{eq:model}
\frac{d W(t)}{dt} &=& b \cdot u(t) \cdot  R(t)  \cdot W(t) - \mu  \cdot W(t)\\ \nonumber
\frac{ d Q(t)}{dt} &=& b \cdot c  \cdot (1 - u(t))  \cdot R(t)  \cdot W(t) - \nu \cdot  Q(t)
\end{eqnarray}

Here  we  start with $W(0) = 1$ and $Q(0) = 0$ as we  are  counting  the 
founding queen as a worker, since she forages to feed the first brood.

Another possible population to consider  is drones. One  needs them to fertilize the queens and they typically come from other hives.  A good hive during the mating season may have 1-10 percent drones that gather daily in clouds that queens seek out to be fertilized.  Queens mate with 5-20 drones and each drone dies after mating. Unfortunately drones do nothing else but eat honey.  They don't gather resources, defend, make wax, take care of brood, etc.  Typically, all the drones are chased out of the hive at the end of the season. In  (\ref{eq:model}) we do not consider drones.

\begin{itemize}
\item $R(t)$  is  the  return function (per unit member of the colony) which  stands  for  the  resource 
abundance and the foraging capability of the colony at time $t$.  Resource abundance is very important for honey bees and is not a constant at all. So this is a first approximation. This  may 
be a measure of the floral diversity in the case of a colony of bees.  Often 
it  can  be  assumed to be constant, but the models are  quite  robust  and 
apply when $R(t)$ is more reflective of broader changes in the ecosystem.
\item $-\mu W(t)$ is an exponential decay term which says, in the absence of  any 
energy input ($R(t)$) the worker population will die off.
\item $u(t)$  is that utilization percentage (between 0 and 1) which  the  colony 
applies to worker production. Turning  to the equation for $Q(t)$ we see that $(1 - u(t))$ then represents  that 
portion  (between  0  and 1) of the resource which the  colony  applies  to 
reproducer or queen production.
\item $b$ stands for a biomass conversion factor - 
i.e.  from the resource energy, how much actually becomes worker  insects per unit time.  
\item $c$  stands for a differential conversion factor, i.e. if $c > 1$ it  takes 
more energy to produce a reproducer insect than to produce a worker insect and if   
  $c  = 1$ the energy needed to produce worker and reproducer  insect are  the  same.  
Often this latter case is assumed. While it  takes more energy to make a queen, e.g., more wax for the cell (larger) and more royal jelly (probably twice as much), for a  queen is larger,  as much as twice as large as a worker insect.  But the energy is not a lot in either case.
\item $-\nu  Q(t)$ represents exponential decay of the reproducer population  in 
the absence of any further energy inputs. NB:  For honey bees the reproducer population is basically 1 or 2.  They do produce several queens at the same time (2 to 8 about) when they get ready to swarm.  But the Queens ``fight'' and only one survives. It is not known  why they do this but it is believed that  it is to manage risk for they must have one (but not more) survive or the hive will die. This is because the mother of these queens typically took off in a swarm a week or so before this event.  Hence you have 2 surviving queens - one home and one gone to start a new hive.  This can repeat several times in a season so there might be up to 4 surviving queens (3 swarms) by season's end.
\end{itemize}

Another factor is honey inventory. When there is a lot then the bees begin to change behavior. They get ready for a swarm and raise more bees.  So maybe one could also  include an ``Inventory'' of resources that is added to with available resources from outside the colony. Perhaps there is an inventory threshold above which behavior changes. One other factor is  that a swarm needs enough bees and honey to take with them to start a successful new colony and leave behind enough bees (and brood) to be viable.  A rough estimate is 50\% of the bees go in a swarm and good colony size before a swarm is about 50,000 bees. There is a LOT of variation here.  
\vspace*{.3 true cm}

Now  the question is this, 

\begin{center}
\parbox{9.6 true cm}{``What  $u(t)$ should the colony use so that the  number 
of reproducers, $Q(T)$, at  season's end (when $t = T$) is maximized?''}
\end{center}

Determining  this function $u(t)$ is in the province of a discipline  known 
as {\it optimal control theory\/}.  This  is one of a class of problems which seeks to 
find the value of $T$ so that  some functional value, in this case $Q(T)$,  is optimized while subject  to 
some dynamic system, in this case differential equations (\ref{eq:model}) with initial conditions.

Consider  these values of the parameters: $b = 0.0013$, $c = 1$, $\mu  =  0.022$, 
$\nu = 0.005$, $T = 205$ - days, and $R(t) = 50$. Recall $u(t)$ is a function (in this case, called an {\it admissible control function\/}, actually)  that is between 0 and 1, where $u(t) = 0$, means at time $t$ the colony puts 0\% of its energies into worker production; $u(t) = 1$, means at time $t$ the colony puts 100\% of its energies into worker production; and, for example, $u(140) = 0.31$, means at time $t=140$ the colony puts 31\% of its energies into worker production.


\begin{enumerate}
\item  {\bf Model - try outs!}

For each  model below, i.e.~$u(t) = u_1(t)$,  $u(t) = u_2(t)$, and $u(t) = u_3(t)$:
\begin{itemize}
\item [a)] plot $u(t)$ in the time interval $[0, 205]$;
\item [b)] describe  what the colony is doing;
\item [c)] compute the number of Queen insects at the end of the season, i.e. $Q(205)$; and
\item [d)] plot both $Q(t)$ and $W(t)$ in the time interval $[0, 205]$ - with variables identified.

\begin{eqnarray*}
u_1(t) &=&0.5 + 0.5 \sin \left(\frac{\pi}{205} t\right)\\
u_2(t) &=&\frac{1}{205} t\\
u_3(t) &=&1- \frac{1}{205} t\\
\end{eqnarray*}
\end{itemize}
\item  Now pick one {\bf strategy of your own}  - different in form and substance from those offered in (1) - and model it with a function $u(t)$. Proceed to analyze it in the same manner as (a) - (d) in (1) above. Explain what you were trying to model in your selection of $u(t)$.





\item {\bf Bang-Bang Solution} 
 Quite  often the optimal solution for $u(t)$ is called a Bang-Bang  Theory, 
because $u(t)$ is all or nothing, (on or off) in this case 1 or 0, respectively.

There is a switching that occurs from not raising any queens to raising a few queens.  However, there are way too many bees to have them devote all their energy to raising queens. You only need a few bees to raise queens, say 5 to 100 out of 50,000. It is probably  driven by the rate at which a worker can produce royal jelly and how much a queen needs. The remainder of the bees deal with gathering more nectar, making wax cells, fending off intruders, gathering pollen, cleaning out the hive, etc.  Thus another model might consider that  there are two levels of $u(t)$: low ($u(t)=0$) and high ($u(t) \approx 0.0001\cdot W(t)$).  It may be  that the energy to make a queen is not a big issue. It is  more about number of bees, honey inventory, brood inventory, natural resource levels, season (you need lots of drones out there from other hives), temperature, etc.
\vspace*{.3 true cm}
Now can we determine an optimal switching time $t_s$ when the colony should  go 
from producing all worker insects ($u(t) = 1$ for $t < t_s$, so $u(t) = 1 - {\rm UnitStep}[t-t_s]$) to producing all 
reproducers (queens) ($u(t) = 0$ for $t > t_s$)?

\begin{itemize}
\item [a)] Show what  the curve  of 
$Q(T=205)$ vs.~$t_s$ looks like, i.e. try out a number of switching times, $t_s$ in the time interval $[0, 205]$ (over a broad range and over a narrow range about your optimal value, collect the resulting $(t_s, Q(205))$  data points into a set,  called data, and plot these  points to see what an optimal switching time might be. 

\item [b)] Using the colony's optimal switching time you have discovered in (3), show how this compares to all the results from the  Model - try outs! section above and your own strategy in (2).
\end{itemize}

\item   Explore the effect on our optimization strategy from (3)  if we change our $c$ value - recall what the significance of $c$ is - see above. For instance, suppose we alter our problem from $c = 1$ to, say, $c = 2$, or $c = \frac{1}{2}$ (be sure to explain what each of these means in terms of the colony and the model) then what effect do such changes have on our optimal switching time? 

\item   Explore the effect on our optimization strategy from (3)  if we change our $R(t)$ function to more closely model the availability of resources as modeled by the seasonal (in this case 205 days from early spring to late fall)  return function. Propose a function for $R(t)$  which is different from $R(t)=50$  which you believe might be one seasonal resources model. Then compare your model results with those of (3).  Explain any differences or similarities.

\item Some species (e.g., bees) are known to swarm or get up and leave the hive, as the reproducers start new colonies in the middle of the season. Model this behavior with the parameters, $b = 0.0013$, $c = 1$, $\mu  =  0.022$, 
$\nu = 0.005$, $T = 205$ - days, and $R(t) = 50$. Then compare this with your results (i.e. number of reproducers available at end of the season) in your optimal switching time discovered in the bang-bang Solution section above.  Make any additional assumptions you need.
\end{enumerate}

\teacher{
\section*{NOTES FOR TEACHERS}

We were  tempted to call this modeling scenario, ``Optimal Control Applied to Insect Colonies -- Differential Equation Modeling at its Best?'' with emphasis on ``Best'' since this is an optimization modeling opportunity. However, we stayed more basic in Insect Colokny Survival Optimization.

In the paper \cite{MaceviczOster1976} the authors offer a model for utilization of energy input by a social insect colony with the objective of maximizing the number of reproducer insects at the end of the season subject to dynamic constraints of a system of differential equations. This is an example of the Maximality Principle of Optimal Control applied to differential equations in ecology to fit with the evolutionary theory of the selfish gene.

We have assigned this modeling scenario  for a number of years in several different approaches.   We have gone in cold turkey, i.e. having the class build the model entirely from scratch to giving them the model and asking them to play with it. We could let them build the terms for the energy diversion options into worker and reproducer insects toward the goal of maximizing the end of season number of reproducers or give them the terms and ask them to play with them. We do the latter in this presentation.

For honey bees the objective may be to  maximize the number of colonies that survive the winter, so there must be  considerations for honey and pollen inventory, number, size, and timing of swarms (they need time to build a hive), number of bees (divided among worker, drones, queen), and natural resources (honey, water, pollen) availability, and mating success.

We believe it is important to permit students in a differential equations class to see beyond the course and to see where this material can be applied in various venues. This work in optimal control and ecology is perfect for meeting those goals.

We usually spend about a period in class, especially if students are to be lead to discovery of (\ref{eq:model}), but principally to get them started on their inquiry so they can try out some strategies for different control functions, $u(t)$. 


We offer a file 6-5-Discovery-InsectColonyOpt-Article.pdf which is a reprint of our article about broader issues in discovering applications for our classroom \cite{Winkel2013}, including this study, as well as comments on student responses.


\noindent From \cite{Winkel2013} we quote,

\begin{quote}
We present Pontryagin's Maximal Principal \cite{PontryaginEtAl1962}  as a necessary condition for an optimal control of our
system (\ref{eq:model}), much like the condition that an optimum occurs for a function of one or more variables
where the derivatives are zero for suitably defined functions. However, we offer no proof. Indeed, we
just indicate that this principal can be applied to prove the existence of a bang-bang solution, i.e. 
the control function $u(t)$ is either 0 or 1 with it being 0 first and then at some optimal time (called
the switching time) switching to 1. Furthermore, we note for our students that optimal control is a
very powerful concept and tool in engineering and science, and would be worthwhile pursuing in their
coursework.

We do not derive the bang-bang solution, but rather we point to the literature. An appendix to
\cite{MaceviczOster1976} offers an explicit derivation of the bang-bang solution's switching time for this insect colony model.
An excellent, more general exposition of optimal control theory, with rich examples, can be found in
\cite{Berkovitz1974}. Students in our classes will do several projects and some will want to revisit the notion of optimal control. We do not send them to original sources such as Pontryagin's Maximal Principle \cite{PontryaginEtAl1962} which guarantees the existence of an optimal control function, indeed one which is bang-bang. We remind
them of what they have just inferred experimentally in their own models. We finally discuss the role
of optimization in evolutionary biology modeling.

We believe it is incumbent upon us as teachers to share our discoveries with our students, be it our
own research or that which we learn in our own continuing education. Moreover, it is very important
to point beyond where we are in every mathematics class to give students a sense of the direction of
mathematics and of its elegant capabilities to solve real problems in our world. We believe this short
narrative does all this and we encourage you to take your own personal journey of discovery, to bring
your students into that purview, and then to share your narrative through publication.
\end{quote}
}

\begin{thebibliography}{99}

\bibitem{Berkovitz1974} Berkovitz, L. D. 1974. {\it Optimal Control\/}. New York:  Springer-Verlag.

\bibitem{MaceviczOster1976} Macevicz,  S.  and G. Oster. 1976.  Modeling Social Insect Populations II: Optimal Reproductive Strategies in Annual Eusocial Insect Colonies. {\it Behavioral Ecology and Sociobiology\/}. 1: 265-282.

\bibitem{PontryaginEtAl1962} Pontryagin, L., V. Boltyanskii, R. Gamkrelidaze, and E. Mishchenko. 1962. {\it The Mathematical Theory of Optimal Processes\/}. New York:  John Wiley and Sons. 


\bibitem{Winkel2013} Winkel, B. J. 2013. Discover for Yourself - An Optimal Control Model In Insect
Colonies. {\it PRIMUS\/}. 23(5): 459-466.



\end{thebibliography}


\end{document}
