\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
Eighth Mississippi State - UAB Conference on Differential Equations and
Computational Simulations.
{\em Electronic Journal of Differential Equations},
Conf. 19 (2010),  pp. 85--98.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2010 Texas State University - San Marcos.}
\vspace{8mm}}

\begin{document} \setcounter{page}{85}
\title[\hfilneg EJDE-2010/Conf/19/\hfil
Importance of brood maintenance terms]
{Importance of brood maintenance terms in simple
models of the honeybee - Varroa destructor - acute bee paralysis
virus complex}

\author[H. J. Eberl, M. R. Frederick, P. G. Kevan \hfil EJDE/Conf/19 \hfilneg]
{Hermann J. Eberl, Mallory R. Frederick, Peter G. Kevan}  % in alphabetical order

\address{
Department of Mathematics and Statistics\\
University of Guelph, Guelph, ON, N1G 2W1, Canada}
\email[Hermann J. Eberl]{heberl@uoguelph.ca}
\email[Mallory R. Frederick]{mfrederi@uoguelph.ca}

\address{Peter G. Kevan \newline
School of Environmental Sciences\\
University of Guelph, Guelph, ON, N1G 2W1, Canada}
\email{pkevan@uoguelph.ca}

\thanks{Published September 25, 2010.}
\subjclass[2000]{92D25, 92D30}
\keywords{Honeybees; varroa destructor; acute bee paralysis virus;
\hfill\break\indent
 colony collapse disorder; wintering losses; mathematical model}

\begin{abstract}
 We present a simple mathematical model of the infestation of
 a honeybee colony by the Acute Paralysis Virus, which is
 carried by parasitic varroa mites (\emph{Varroa destructor}).
 This is a system of nonlinear ordinary differential equations
 for the dependent variables: number of mites that carry the virus,
 number of healthy bees and number of sick bees.
 We study this model with a mix of analytical and computational
 techniques. Our results indicate that, depending on model parameters
 and initial data, bee colonies in which the virus is present can,
 over years, function seemingly like healthy colonies before they
 decline and disappear rapidly (e.g. Colony Collapse Disorder,
 wintering losses). This is a consequence of the fact that a
 certain number of worker bees is required in a colony to maintain
 and care for the brood, in order to ensure continued production
 of new bees.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{lemma}{Lemma}[section]
\allowdisplaybreaks

\section{Introduction}

A famous folklore quote says that four years after honeybees vanish
from the Earth, mankind will vanish too. Commonly attributed to Albert
Einstein, there seems to be little evidence that it is
authentic \cite{S:2007}. Nevertheless, the
message is clear: No bees, no pollination, no crops, no people.

Besides cows, pigs, and poultry, western honeybees (\emph{Apis
mellifera}) are among the economically most important domestic
livestocks in Europe and Northern America, not to mention their
ecological importance \cite{T:2008}.   Due to the demand of honey
and other hive products, but in particular because of the demand
for pollination services, apiculture has developed to become a
profitable business. Honeybees (\emph{Apis spp.: Hymenoptera:
Apidea}) live in highly integrated and complex social colonies to
the extent that a single colony can be thought of as a super
organism \cite{T:2008}.  Like all animals, honeybees are hosts to
parasites and pathogens.  For the temperate zone biotypes of
western honeybees, the devastating attacks of parasitic mites,
especially \emph{Varroa destructor}, have created major problems
for beekeepers around the world.  The damage caused by the mites
has been exacerbated by mite-borne infections of viruses
\cite{KHOG:2006}, notably the Acute Paralytic Bee Virus (APV,
ABPV). In recent years, beekeepers in the USA, and now other parts
of the world, have reported huge losses of colonies to the extent
that there is national, continental and international alarm over
their demise, e.g. \cite{NRC:2007}.  The combination of parasitic
mites and virus infections have been indicated as a major cause of
Colony Collapse Disorder (CCD) as the syndrome is diagnosed in
USA. Elsewhere, symptoms are not exactly the same as in CCD, but
huge and worrisome losses, in particular wintering losses, are
also recorded with the same parasite/virus complex suggested as a
major cause.  It has been suggested that various other stresses,
such as hard management, pesticides, weather, and poor diet may
also be involved to make matters worse \cite{KGSvE:2007}.

The Acute Bee Paralysis virus kills larvae, pupae, and adults only
in association with varroa mites \cite{KHOG:2006,M:2001}. In
contrast to other viruses, such as the \emph{Deformed Wing Virus},
larvae and pupae that are affected by the virus die before they
develop into bees \cite{M:2001, SM:2004}; i.e., all bees that are
infected with the APV virus must have acquired it as adults.


In this small study we present and analyze a simple model of
three ordinary differential equations for the disease dynamics of
the bee-mite-APV complex, that follows the classical SIR-like
approach of mathematical epidemiology. It will be formulated for
the dependent variables: (i) number of mites that carry the virus,
(ii) number of healthy bees and (iii) number of infected bees.
Our work builds firmly on the previous study \cite{SM:2004}, but
introduces a simple but important extension: In their original
model \cite{SM:2004}, the authors assumed that the birth rate of
bees might be affected by the presence of the virus in the colony,
but it is not directly linked to the size of the worker bee
population. In other words, even if all healthy bees are gone, the
pupae in the brood will develop into adults. Moreover, since it is
assumed that the queen bee is not affected by the virus, this
would allow a constant birth of new bees during spring, summer,
and autumn. However, carrying for and maintaining the brood is an
important task in bee hives, in which many workers are involved
\cite{K:2010,T:2008,W:1987}. If the bee population becomes too
small, larvae and pupae cannot be reared and die in the brood
cells before they can develop into adult bees. To reflect this in
the model, we introduce an additional brood maintenance term in
the birth rate which depends on the current size of the worker
population. From a mathematical point of view, including such a
maintenance term affects the structure of the omega limit set of
the system, because it allows for additional equilibria. Moreover,
it can change the stability of the existing equilibria. Like
\cite{SM:2004}, we will study our model separately for the four
seasons. Additionally we will use computer simulations to
investigate the fate of the bee colony through the years according
to our model.

\section{Mathematical Model}

The mathematical model is formulated in terms of the dependent
variables
\begin{trivlist}
\item{$m$:} number of mites that carry the virus,
\item{$x$:} number of honeybees that are virus free,
\item{$y$:} number of honeybees that are infected with the virus.
\end{trivlist}
It reads
\begin{gather}
\frac{dm}{dt} = \beta_1 (M-m) \frac{y}{x+y}
  - \beta_2 m \frac{x}{x+y} \label{mode}\\
\frac{dx}{dt} = \mu g(x)h(m) - \beta_3 m \frac{x}{x+y}
 - d_1 x \label{xode}\\
\frac{dy}{dt} = \beta_3 m \frac{x}{x+y} - d_2y \label{yode}
\end{gather}
In equation \eqref{mode}, the parameter $M$ denotes the number of
mites in the bee colony. With the same arguments as in
\cite{SM:2004}, we will treat $M$ as a parameter rather than a
dependent variable; in particular, $M$ is independent of $m$. By
treating $M$ as a parameter we implicitly assume that the mite
population reaches its carrying capacity very rapidly; i.e., we
make a quasi-steady state argument.

The parameter $\mu$ in \eqref{xode} is the maximum birth rate,
specified as the number of worker bees born per day.

The parameter $\beta_1$ in \eqref{mode} is the rate at which mites
that do not carry the virus acquire it. The rate at which infected
mites lose their virus to an uninfected host is $\beta_2$. The
rate at which uninfected bees become infected is $\beta_3$, in
bees per virus carrying mite and time.

Finally, $d_1$ and $d_2$ are the death rates for uninfected and
infected honeybees. We can assume that infected bees live shorter
than healthy bees, thus $d_2 > d_1$.

The function $g(x)$ expresses that a sufficiently large number of
healthy worker bees is required to care for the brood. The
inclusion of this term is the only difference between the model
studied here and the model in \cite{SM:2004}, on which this study
is based.
 We think of $g(x)$ as a switch function. If $x$ falls below a
critical value, which may seasonally depend on time, essential
work in the maintenance of the brood cannot be carried out anymore
and no new bees are born. If $x$ is above this value, the birth of
bees is not hampered. Thus $g(0,\cdot)=0$, $\frac{dg(0)}{dx}\geq
0$, $\lim_{x\to \infty} g(x)=1$. A convenient formulation of such
switch like behavior is given by the sigmoidal Hill function
\begin{equation}\label{gswitch}
g(x)=\frac{x^n}{K^n+x^n}
\end{equation}
where the new parameter $K$ is the size of the bee colony at which
the birth rate is half of the maximum possible rate and the
exponent $n>1$, see also Figure \ref{g_fig}. If $K=0$ is chosen,
then the original model of \cite{SM:2004} is obtained. In this
case the brood is always reared at maximum capacity, independent
of the actual bee population size, since $g(x)\equiv 1$.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig1}
\end{center}
\caption{Brood maintenance coefficient $g(x)=\frac{x^n}{K^n+x^n}$
as a function of bee population size for various choices of $K$
and $n$}\label{g_fig}
\end{figure}

The function $h(m)$ in \eqref{xode} indicates that the birth rate
is affected by the presence of mites that carry the virus. This is
in particular important for viruses like APV that kill infected
pupae before they develop into bees. The function $h(m)$ is
assumed to decrease as $m$ increases, $h(0)=1$,
$\frac{dh}{dm}(m)<0$ and $\lim_{m\to \infty} h(m)=0$;
\cite{SM:2004} suggests that this is an exponential function $h(m)
\approx e^{-m k}$, where $k$ is a non-negative function. We will
use this in the computer simulations later on.

The parameters $M$,$\mu$, $\beta_i$, $d_i$, $g(x)$, $h(m)$ are
assumed to be non-negative. They can change with time.  In
particular major differences may be observed between seasons. For
example, the life span of a worker bee in summer is much shorter
than in winter \cite{AO:2002,O:1988}; the birth rate for bees is
higher in summer than in spring and autumn, and it drops down to 0
in winter \cite{T:2008}. Seasonal averages for the model
parameters $\beta_{1,2,3}, \mu, d_{1,2}$ can be derived from the
data in \cite{SM:2004}. These are summarized in Table
\ref{para_table}. We will treat the remaining parameters $M$ and
$K$ as unknown in the next section of this study. We will
investigate the long term behavior and fate of the colony in
dependence of these two parameters. The size of the mite
population $M$ is a an obvious parameter to choose for such a
study. This will give us insight in how strong an infestation can
be fought off by a bee colony. The brood maintenance coefficient
$K$ is chosen as a free parameter, because we do not have a good
estimate. We will assume it to be in the order of several
thousands; in \cite{SM:2004} it is mentioned that at least 3000
worker bees are required in the beginning of spring to maintain
the brood. In particular, we are interested how the maximum mite
population size $M$ that can be tolerated by the bee colony
depends on $K$.


\begin{table}
\caption{Seasonal averages of model parameters, derived from
the data presented in \cite{SM:2004}.}\label{para_table}
\begin{center}
\begin{tabular}{ccccc}
Parameter & Spring & Summer & Autumn & Winter \\
\hline
$\beta_1$ &0.1593&0.1460&0.1489& 0.04226\\
$\beta_2$ &0.04959&0.03721&0.04750&0.008460\\
$\beta_3$ &0.1984&0.1460&0.1900&0.03384\\
$d_1$ &0.02272&0.04&0.02272&0.005263\\
$d_2$ &0.2&0.2&0.2 & 0.005300\\
$\mu$ &500&1500& 500& 0\\
\hline
\end{tabular}
\end{center}
\end{table}

For the discussion of the model in the subsequent sections,  the
following result will be helpful.

\begin{lemma}\label{pos_inv_lemma}
There exists a $\tilde R >0$ such that the set
\[
\mathcal{K}_R := \{ (m,x,y) \in \mathbb{R}^3 : \;  0 < m< M_\infty,\; x>0, \; y>0, \;x+y<R \}
\]
is positively invariant under \eqref{mode}-\eqref{yode} for all
$R>\tilde R$, where by $M_\infty$ we denote $M_\infty:= \max_t
M(t)$.
\end{lemma}

\begin{proof}
We apply the usual tangent criterion \cite{W:2000}. At $m=0$ we
have $\frac{dm}{dt}=\beta_1M\frac{y}{x+y}>0$ and at $m=M$ we have
$\frac{dm}{dt}=-\beta_2 M \frac{x}{x+y}<0$ for $x>0$, $y>0$. Thus
solutions with $0< m(0) < M$ will satisfy $0<m(t)<M$ for all
$t>0$. Similarly, for $y=0$ we have
$\frac{dy}{dt}=\beta_3m\frac{x}{x+y}>0$ and for $x=0$ we have
$\frac{dx}{dt}=0$. Thus our system \eqref{mode}-\eqref{yode} is
positivity preserving.

The remaining boundary of the set $\mathcal{K}_R$ is the plane $R-x-y=0$, 
which has the outer normal vector $n=(0,1,1)^T$.
To apply the tangent criterion we have to show that along $R-x-y=0$
the inequality
\[
\Big(\frac{dx}{dt} + \frac{dy}{dt}\Big)_{R-x-y=0}
=\mu g(x)h(m)  - d_1 x  - d_2 y \leq 0
\]
holds for all large enough $R$; i.e., that
\begin{equation}\label{pos_inv_crit}
\mu g(x)h(m)  +(d_2-d_1) x  - d_2 R \leq 0,\quad\text{for all }
 x\leq R.
\end{equation}
As noted above, we can assume $d_2>d_1$. Thus we obtain with $x<R$
and $0 \leq g(x)\leq 1$, $0\leq h(m)\leq 1$ from
\eqref{pos_inv_crit} the inequality
\[
\mu g(x)h(m)  +(d_2-d_1) x  - d_2 R < \mu -d_1R.
\]
Hence the assertion \eqref{pos_inv_crit} holds certainly for
$R>\tilde R= \mu/d_1$.
\end{proof}

Note, however, that this inequality is not sharp; i.e.,
 one can even find a smaller $\tilde R$ with the same property.

Thus, all solutions of \eqref{mode}-\eqref{yode} that start in
$\mathcal{K}_R$, $R>\tilde R$, remain there. In particular they
are positive and bounded from above by constants. Moreover, even
if initially $m(0)>M$, the solution will be attracted by
$\mathcal{K}_R$.

\section{The Autonomous Case}\label{auto_sec}

We start our investigation of the disease dynamics model
\eqref{mode}-\eqref{yode} with the simple case of constant
coefficients as in \cite{SM:2004}, pointing out that these results
only have quantitative meaning if the characteristic time scale of
the dynamics is shorter than the period over which the parameters
can be assumed to be constant. However, in many cases, e.g. Table
\ref{para_table}, parameter values are only available in form of
seasonal averages, and this is the situation that we have in mind.
We have to distinguish between three types of equilibria of
\eqref{mode}-\eqref{yode}.

\subsection{Equilibrium points}

{\bf Disease free equilibrium.}
This is the equilibrium that is attained by an entirely
healthy population,
\begin{equation}
m=0, \quad x=x^\ast, \quad y =0
\end{equation}
where $x^\ast$ is a root of the function
\begin{equation}\label{xast_root}
F(x)=\mu g(x) -d_1 x.
\end{equation}
For maintenance functions of the form $g(x)= \frac{x^n}{K^n + x^n}$
this reduces to finding the positive roots of the polynomial
of degree $n$, $G(x)=x^n -\frac{\mu}{d_1} x^{n-1} + K^n$.
According to Descarte's rule of signs, there are at most two such
roots.
In the case $n=2$ these are easily calculated as
\begin{equation}\label{xast_eq}
x^\ast_{1,2} = \frac{1}{2}\Big(\frac{\mu}{d_1} \pm
\sqrt{ \frac{\mu^2}{d_1^2} - 4 K^2}\Big).
\end{equation}
They exist  if $\mu/d_1 > 2K$. The Jacobian of the  system
\eqref{mode}-\eqref{yode} in general form reads
\begin{equation}\label{Jacobian}
J(m,x,y)=
 \begin{bmatrix}
-{\frac {\beta_1 y+\beta_2 x}{x+y}}&
-{\frac {y \left(\beta_1(M-m)+\beta_2 m\right) }{ \left(
x+y \right) ^{2}}}&{\frac {x
\left( \beta_1(M-m)+\beta_2 m \right) }{ \left( x+y \right) ^{2}}}
\\
\mu g(x) h'(m)-{\frac {\beta_3 x}{x+y}}&
\mu g'(x) h(m) -\frac{\beta_3 my}{(x+y)^2} -d_1
&{\frac {\beta_3 mx}{\left( x+y \right) ^{2}}}\\
 {\frac {\beta_3x}{x+y}}&
 {\frac {\beta_3 my}{ \left( x+y \right) ^{2}}}&
 -{\frac{\beta_3 mx}{(x+y)^2} -d_2}
 \end{bmatrix}.
\end{equation}
In the case of the disease free equilibrium this reduces to
\begin{equation}
J(0,x^\ast,0)=\begin{bmatrix}
-\beta_2&  0&  \frac{\beta_1M}{x^\ast} \\
d_1 {x^\ast} h'(0)-\beta_3&
\mu g'(x^\ast) -d_1 &0\\
\beta_3 & 0 &  -d_2
 \end{bmatrix}.
\end{equation}
Its eigenvalues are
\begin{equation}
\lambda_1 = \mu g'(x^\ast)-d_1,\quad
\lambda_{2,3} = -\frac{1}{2}\Big(\beta_2+d_2 \mp
\sqrt {(\beta_2-d_2)^2 +4\frac{\beta_1 \beta_3 M}{x^\ast}}\Big).
\end{equation}
We note that the eigenvalues depend on the specific form of the
brood maintenance term $g(x)$ via $x^\ast$ and \eqref{xast_root}.
For the special case $g(x)=\frac{x^2}{K^2+x^2}$, the eigenvalue
$\lambda_1$ takes the simpler form
\begin{equation}
\lambda_1= \mp d_1 \sqrt{1-4 K^2\big(\frac{d_1}{\mu}\big)^2}.
\end{equation}
Particularly it is real. The disease free equilibrium with
$x=x^\ast_2$ is always unstable, because $\lambda_1>0$. On the
other hand, for $x^\ast_1$,  the associated eigenvalue $\lambda_1$
is negative. The eigenvalues $\lambda_{2,3}$ are always real,
$\lambda_3$ is always negative. Therefore, it is the eigenvalue
$\lambda_2$ that decides the stability of the disease free
equilibrium with $x=x^\ast_1$. Stability is achieved if it is
negative, i.e. if
\begin{equation}\label{Mcrit}
M < \frac{\beta_2 d_2 x^\ast_1}{\beta_1 \beta_3}=: M_{\rm crit},
\end{equation}
where $x^\ast_1$ depends on birth and death rate of the bee
population as well as on $K$. We plot $M_{\rm crit}$ as a function of
$K$ in Figure \ref{KMplot_fig} for the parameters in Table
\ref{para_table}. We note that the critical value for the model
\cite{SM:2004} is obtained for the special case $K=0$. Thus, if
one includes the consideration  that in order to rear the brood, a
critical bee population size is needed, a drastic reduction in the
maximum number of mites that can be fought off by the colony is
observed.

\begin{figure}[ht]
\includegraphics[width=0.6\textwidth]{fig2}
\caption{Critical mite population $M_{\rm crit}$ as a function of 
the half maintenance population size $K$.}\label{KMplot_fig}
\end{figure}


\begin{figure}[ht]
\includegraphics[width=0.7\textwidth]{fig3a}
\includegraphics[width=0.7\textwidth]{fig3b}
\caption{Long-term behavior of model \eqref{mode}-\eqref{yode} in
dependence of number of mites $M$ and maintenance size of the
population $K$. Plotted are healthy bees $x$ (green) and infected
bees $y$ (red) at steady state.}\label{longterm_fig}
\end{figure}


{\bf Collapse equilibrium.}
This equilibrium describes the vanishing of the bee population.
These are the points
\begin{equation}
m=m^\ast, \quad x=0, \quad y =0
\end{equation}
where $m^\ast>0$ can be any positive number. Moreover,
since for $m(0)<M$ we will always have $m(t)\leq M$.
Therefore, we are only interested in $0<m^\ast<M$.  Of course,
from a practical point of view, the absence of bees will imply
the absence of mites, and thus $m=0$ as a consequence.

The right hand side of \eqref{mode}-\eqref{yode} does not satisfy
a Lipschitz condition for $x=y=0$, and some of the coefficients of
the Jacobian \eqref{Jacobian} blow up. Using the same Lyapunov
argument as in Lemma \ref{pos_inv_lemma}, we note that inequality
\eqref{pos_inv_crit} is also satisfied for all sufficiently small
$r>0$ and all $x<r$. To show this, we need to show again for
sufficiently small $r$ that
\begin{equation}\label{pos_inv_crit_small}
\mu g(x)h(m)  +(d_2-d_1) x  - d_2 r \leq 0,\quad\text{for all }
 x\leq r.
\end{equation}
By continuity and monotonicity of $g(x)$, and because
of $h(m)\leq1$, we obtain that this is satisfied if
\[
\mu g(r)-d_1 r \leq 0.
\]
The last inequality is certainly satisfied for all $r$ with
$0<r\leq x^\ast_2$ from \eqref{xast_eq}\footnote{Note: with the
same argument, a better estimate for $\tilde R$ in Lemma
\ref{pos_inv_lemma} than $\tilde R=d_1/\mu$ is obtained as $\tilde
R:=x^\ast_1$}. Thus, $\mathcal{K}_r$ is positively invariant for
$0<r\leq x^\ast_2$. Therefore, after small perturbations, the
solutions remain close to the collapse equilibrium, which is
therefore stable. We point out that the unconditional stability of
the collapse equilibrium stems from the fact that a large enough
bee population is required to care for the brood in our model. In
models that do not make this assumption this is not the case. The
model of \cite{SM:2004} is the special case $K=0$ in \eqref{xode},
in which  $x_2^\ast$ vanishes.


If \eqref{Mcrit} is satisfied, we have several stable equilibria
with $y^\ast=0$ (as well as the always unstable steady state
$(0,x^\ast_2,0)$ above). Whether one of these equilibria is
attained, and if so, which one, therefore, depends on the initial
strength of the colony. If the initial bee population lies within
$\mathcal{K}_{\tilde r}$, which depends on the parameters, then a
collapse equilibrium is attained. Otherwise the disease can be
fought off.

{\bf Endemic equilibrium.} In this state, the colony contains
healthy and sick bees. The equilibrium generally takes the form
\begin{equation}
m=m^\ast, \quad x=x^{\ast\ast}, \quad y=y^{\ast\ast}.
\end{equation}
The algebraic treatment of this steady state leads to involved
expressions that do not lend themselves to an insightful analysis.
While the eigenvalues of a three-dimensional system can be
explicitly calculated, this involves look-up tables to distinguish
between a multitude of possible cases. Instead of presenting these
calculations here, which can be performed by computer algebra
software, we demonstrate the long-term behavior of the system in
computational simulations, where we use for $\beta_{1,2,3},
d_{1,2}, \mu$ the parameters in Table \ref{para_table}, for spring
and summer conditions.

\subsection{Computer simulations}

In Figure \ref{longterm_fig}, we plot the populations size of
healthy and infected bees, $x$ and $y$, at steady state (all
solutions do converge to a steady state), in dependence of the
number of mites $M$ and the population size that is required to
maintain the population, $K$, for the maintenance function
$g(x)=\frac{x^2}{K^2+x^2}$. This is for the initial conditions
$m(0)=0.2 M$, $x(0)=3000$, $y(0)=0$. We show the results for
spring and summer conditions. As can be seen from Table
\ref{para_table}, the fall parameter values are very similar to
those in spring and one obtains very similar results.

In both cases, spring and summer, we observe that for small values
of $K$ (i.e. higher birth rates of healthy bees)  and $M$ (i.e.
less risk to catch the disease) the disease free equilibrium is
attained, as expected from criterion \eqref{Mcrit}. Equilibria of
the type $x=x^{\ast\ast}>0$, $y=y^{\ast\ast}>0$ are found for
small values $K$ and large values $M$. However, in this case, the
bee population $x+y$ is small compared to the size of an entirely
healthy population. Thus, we conclude that such endemic colonies,
while in a stable equilibrium, are not functioning like healthy
populations. For large values of $K$ (i.e. a strong colony is
required to maintain the brood) and $M$ (i.e. high risk to be
infected), the population dies out and attains the collapse
equilibrium. Moreover, the unconditional stability of the collapse
equilibrium also implies that colonies which are initially very
small will never develop into healthy colonies, because they will
never reach the required population strength for brood maintenance
(not shown in our simulations). This is in particular important
for the spring season, because in winter no bees are born but the
colony naturally declines due to bee death. If at the end of the
winter the colony is not strong enough to rear the brood, it will
die out.

We note that under winter conditions, the only and asymptotically
stable equilibrium is the trivial equilibrium $x^\ast=y^\ast=0$.
However, winter worker bees live much longer than spring, summer
or fall bees, namely several months compared to several weeks, cf
also \cite{AO:2002}. Hence the winter bee death rate is small,
therefore, decline in the population is slow, compared to the
duration of the season. Therefore, long term analysis of the model
for winter conditions, as we conducted here for the case of
constant parameters, does not lead to any useful insight. A bee
colony that is large enough at the beginning of winter that, by
the end of winter it is big enough to recover, will survive. A
healthy population will decline exponentially according to
\[
\frac{dx}{dt}= -d_1x.
\]
If infected bees are present at the beginning of the winter, this,
however, gives only an upper estimate. A lower estimate is
obtained from
\[
\frac{dx}{dt}+\frac{dy}{dt}=-d_1x-d_2y > -d_2 (x+y).
\]

\section{Bee colonies through the seasons: the non-autonomous case}

A rigorous qualitative analysis of the model for time varying
parameters is essentially more complicated than the autonomous
case, even if we assume that the parameters are time periodic with
a period of one year. The quantitative data that are available to
us are seasonal averages. In principle we can construct an
infinite number of non-negative periodic functions with the same
seasonal averages. However, on an objective basis, it is not
possible to decide which of these interpolated functions is more
realistic or better than other ones. Therefore, we simply assume
the model parameters to be constant across four seasons of 91 days
each, at their average values, but varying between seasons, cf
Table \ref{para_table}. The model parameters are repeated in every
year. In order to explore the dynamics of
\eqref{mode}-\eqref{yode} through the seasons we conduct computer
simulations. The inclusion of the colony maintenance term $g(x,t)$
can cause the model to become stiff for small bee colony sizes. In
order to deal with this situation we implemented a simple first
order Rosenbrock-Wanner method for the integration of our
governing equations.

We assume all parameters but $K$ and $M$ to be given and we
investigate the effect of these two parameters on the disease
dynamics. Moreover, the analysis above has shown several stable
steady states can be found for the autonomous system. Therefore,
it can be expected that the initial data might play a role for the
fate of the bee colony in the simulations.

A numerical exploration cannot give a complete picture of the
dynamics of the model. Instead we aim at illustrating the possible
outcomes for the model solutions in examples.

\subsection{Simulation I: effect of size of mite population $M$}

In a first simulation experiment we investigate the effect of the
number of mites, $M$, on the fate of the bee colony. The analysis
of the autonomous case has shown that there is a critical number
of mites, $M_{\rm crit}$ [cf. \eqref{Mcrit}  above]. For mite
populations smaller than $M_{\rm crit}$ the virus can be fought off,
while for values greater than $M_{\rm crit}$ the disease will take its
toll and the bee population will die out.

In this simulation we test for mite populations $M$ somewhat
smaller than this critical value, $0.8 M_{\rm crit} \leq M \leq
M_{\rm crit}$. More specifically, we conduct the following 4
simulations
\begin{itemize}
\item $K:= 0.4 \mu/d_1$ (recall: $K< 0.5 \mu/d_1$ is required
      for $x^\ast_1$ to be stable),
\item $M=\delta M_{\rm crit}$ where $\delta=0.8, 0.85, 0.9, 0.95$.
\end{itemize}
Note that thus defined $K$ and $M$ are piecewise constant periodic
functions of time, since they are defined relative to the other model
parameters.

We assume that initially only a small number of mites carry the
virus. Moreover we assume that initially no sick bees are in the
colony. The number of healthy bees is chosen large enough to
exceed $K$  and  $x^\ast_1$ above.

\begin{itemize}
\item $m(0)=0.01 M$ (1\% of all mites carry virus),
\item $x(0)=x_0$ (10\% larger than $K$ and $x^\ast_1$ of the
       autonomous case),
\item $y(0)=0$ (no bees initially sick).
\end{itemize}
The simulations were run over a period of 6000 days or until the
colony dies out, whichever came first.
The results are plotted in Figure \ref{Mvary_fig}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.44\textwidth]{fig4a}
\includegraphics[width=0.44\textwidth]{fig4b}
 $\delta=0.80$ \hfil $\delta=0.85$ \\
\includegraphics[width=0.44\textwidth]{fig4c}
\includegraphics[width=0.44\textwidth]{fig4d}
 $\delta=0.90$ \hfil $\delta=0.95$ 
\end{center}
\caption{Simulation
of bee-mite-virus dynamics for varying mite population size
$M=\delta M_{\rm crit}$}\label{Mvary_fig}
\end{figure}

We observe that for the smallest value for $M$ in our survey,
$M=0.8 M_{\rm crit}$, the virus infestation is fought off. While for
the first few years virus carrying mites can be found in the
colony, their number decreases over the years. Like the bee
population itself, the number of virus carrying mites is highest
in summer and decreases in fall. The maximum number of $m$
decreases over the first 4 years, after which the virus has
virtually disappeared from the colony. The number of healthy bees
is periodic and seems to be little affected by the presence of the
virus.

As the number of mites $M$ in the population is increased to 85\%
of the critical mite load, the qualitative picture changes. While
the number of virus carrying mites still follows the seasonal
fluctuations, with highest numbers in summer, the maximum number
of virus carrying mites increases from year to year. Over the
first four years the healthy bee population appears to follow a
pattern similar to the previous case, although we do notice that
the maximum population size becomes slightly smaller each year.
After the fourth year of the infestation, the colony vanishes
rapidly, during the spring season. Over the whole period of four
years, the number of virus carrying bees is small compared to the
number of healthy bees. The colony dies because the number of
healthy bees drops below the levels required to maintain the
brood.

Increasing the number of mites $M$ further to 90\% and 95\% of
the critical value, we observe qualitatively the same picture, but
the decline of the bee population to levels that cannot sustain
the brood happens earlier, after two years and one year,
respectively.

The sudden decline of the bee population after several years is  a
consequence of gradual changes. Eventually the population size
falls below the critical threshold that is required for it to
recover, and then the colony dies rapidly.

\subsection{Simulation II: effect of brood maintenance parameter $K$}

In the second simulation experiment we fix the number of mites $M$
at 85\% of the critical value $M_{\rm crit}$ and vary the brood
maintenance coefficient $K$ between 40\% and 80\% of the value
that is necessary to allow for a nontrivial disease free
equilibrium. More specifically

\begin{itemize}
\item $K:= \kappa * 0.5 \mu/d_1$, where $\kappa=0.4, 0.6, 0.7,0.8$,
\item $M=0.85 M_{\rm crit}$.
\end{itemize}
Note that the case $\kappa=0.8$ is identical to $\delta=0.85$ in
the previous simulation experiment.
The initial data are chosen as in the previous simulation experiment:
\begin{itemize}
\item $m(0) =0.01 M$ (1\% of all mites carry virus),
\item $x(0)=x_0$ (10\% larger than $K$ and $x^\ast_1$),
\item $y(0)=0$ (no bees initially sick)
\end{itemize}
As before, the simulations were run over a period of 6000 days
or until the colony dies out, whichever came first.
The results are plotted in Figure \ref{Kvary_fig}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig5a}
\includegraphics[width=0.45\textwidth]{fig5b}
 $\kappa=0.40$ \hfil $\kappa=0.60$ \\
\includegraphics[width=0.45\textwidth]{fig5c}
\includegraphics[width=0.45\textwidth]{fig4b}  \\ 
% figure s85.eps was towice in the original tex file.
$\kappa=0.70$ \hfil $\kappa=0.80$
\end{center}
\caption{Simulation of bee-mite-virus dynamics for varying
$K=0.5 \mu/d_1 \kappa $}\label{Kvary_fig}
\end{figure}

The simulation results are easily summarized: if the number of
bees that is required to maintain the brood is small compared to
the birth/death rate ratio, i.e. if relatively few worker bees are
required for this task, an initially healthy population can fight
off a virus epidemic. Virus carrying mites can be found initially
but, after adjusting for seasonal fluctuations, their numbers
decline. Eventually the colony is disease free. The larger $K$ is,
the longer it will take to rid the colony from the virus. As the
parameter $K$ increases; i.e., if the number of worker bees that
is required to care for the brood becomes large enough, the colony
cannot fight off the disease and eventually dies out, -- 4 years
in our simulations. Again, the number of sick bees $y$ is always
small compared to the overall colony size.

\subsection{Simulation III: effect of initial number of virus
carrying mites $m(0)$}

In a third simulation experiment we fix the number of mites $M$
to be 80\% of the critical mite load, and the brood maintenance
parameter $K$ to be at 60\% of the value that is needed to
establish a stable bee population.
\begin{itemize}
\item $K:= 0.3 \mu/d_1$,
\item $M=0.80 M_{\rm crit}$.
\end{itemize}
In the previous simulations, where the number of mites carrying
the virus initially was small, at $m(0)$ being 1\% of all the
mites, it was found that these parameters are small enough for the
colony to fight of the diseases. In this simulation experiment, we
investigate whether this also holds true if at the starting point
of our simulation more mites carry the virus. We test the initial
data
\begin{itemize}
\item $m(0)=\nu M$ where $\nu=0.1, 0.11, 0.12$,
\item $x(0)=x_0$ (10\% larger than $K$ and $x^\ast_1$),
\item $y(0)=0$ (no bees initially sick).
\end{itemize}
The results are plotted in Figure \ref{m0vary_fig}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig6a}
\includegraphics[width=0.45\textwidth]{fig6b}
 $\nu=0.10$ \hfil $\nu=0.11$ \\
\includegraphics[width=0.45\textwidth]{fig6c} \\ 
 $\nu=0.12$  
\end{center}
\caption{Simulation of bee-mite-virus dynamics for varying
number of initially virus carrying mites $m_0$}\label{m0vary_fig}
\end{figure}

For the smallest value $m(0)$ of initially virus carrying mites
tested, the colony can fight off the disease. The temporal
patterns of the healthy population $x(t)$ and for the number of
mites that carry the virus $m(t)$ are qualitatively similar to the
ones observed before for small enough mite numbers $M$ and brood
maintenance coefficients $K$: $m(t)$ oscillates over the years,
but the peak value in summer decreases until $m$ virtually
disappears; $x(t)$ oscillates over the years as well, with
increasing maximum population sizes as the number of virus
carrying mites decreases.

As $m(0)$ is slightly changed from $m(0)=0.10 M$ to $m(0)=0.11 M$,
the long term behavior changes drastically. The number of virus
carrying mites, after adjusting for seasonal fluctuations,
increases slightly and the number of healthy bees, again after
adjusting for seasonal fluctuations, slightly decreases. After the
third cycle, the healthy bee population enters the spring too
small to be able to fight off the virus, and dies out. When
increasing $m(0)$ further to 12\% of the total number of mites,
this is accelerated and the colony disappears after only two
years. As before, the number of sick bees is small compared to the
overall colony size.

\subsection{Conclusion}

The mathematical model presented here is simple, yet it is too
complex for a complete rigorous analytical treatment. Therefore,
we studied it with a combination of analytical and computational
arguments. It is a straightforward extension of a model that was
previously formulated and analyzed in \cite{SM:2004}. We (i) focus
on the role of a brood maintenance term and how it affects the
disease dynamics of the Acute Paralysis Virus in a bee colony
qualitatively, and (ii) investigate how this affects the bee
colony long-term, over several years.

Our results indicate that models that do not account for the  fact
that a certain colony strength is required to maintain the brood
and thus to sustain  the colony development might underestimate
the maximum  number of mites that a bee population can tolerate.
Our mix of mathematical analysis and computer simulation indicates
that the long term behavior of a bee population can depend on a
variety of parameters, including the initial data. Our computer
simulations have shown that an infested colony might be able to
function seemingly like a healthy colony for years and then
suddenly  can collapse and disappear. In all our simulations, if
this occurs it happens in the spring season because the bee
population is too small to rear the brood. These preliminary
simulation results are in accordance with the recent study
\cite{Getal:2010} on wintering losses in Canada.

Many important processes that can affect the long term fate of  a
colony over the years are not yet included in this model. Examples
are gradual changes of the environment; the size of the mite
population as a dynamically varying variable; the possible loss of
the queen and its replacement by an emergency queen; or swarming
when a new queen emerges. Some such processes can be included in
this model in a rather straightforward manner, however, always at
the expense of introducing new unknown model parameters. Thus,
increasing biological complexity of the model is at the expense of
increasing also the mathematical complexity. Therefore we decided
to keep the model simple in our first study  of the topic.

\subsubsection*{Acknowledgements}
This study was supported in part by the Canadian Pollination
Initiative (CANPOLIN), a NSERC Strategic Network (PI: PGK). 
M. R. Frederick
acknowledges the support received through a NSERC Postgraduate
Scholarship.
H. J. Eberl acknowledges the support received from the
Canada Research Chairs program.

\begin{thebibliography}{99}

\bibitem{AO:2002} Amdam, G.V.; Omholt, S.W.;
  The Regulatory Anatomy of Honeybee Lifespan,
\emph{J. Theor. Biol.} 216(2):209-228, 2002.

\bibitem{Getal:2010} Guzm\'an-Novoa, E.; Eccles, Calvete Y.;
 McGowan, J.; Kelly, P. G.; Correa-Benitez, A.;
Varroa destructor is the main culprit for the death and reduced
populations of overwintered honey bee (Apis mellifera) colonies
in Ontario, Canada, \emph{Apidologie}, 41:443-450, 2010.

\bibitem{KHOG:2006} Kevan, P.G.; Hannan, M.A.; Ostiguy, N.;
Guzman-Novoa, E.;
A Summary of the Varroa-Virus Disease Complex in Honey Bees,
\emph{American Bee Journal}, 2006.

\bibitem{KGSvE:2007} Kevan, P.G.; Guzman, E.; Skinner, A.;
 van Engelsdorp, D.;
Colony collapse disorder in Canada: do we have a problem?
\emph{HiveLights}, May 2007 pp 14-16, 2007.

\bibitem{K:2010} Kevan, P.G.;  Bees,
\emph{Biology and Management},  Enviroquest, Cambridge, ON,
Canada, 2010.

\bibitem{M:2001} Martin, S.J.;
The role of Varroa and viral pathogens in the collapse
of honeybee colonies: a modeling approach,
\emph{J. Appl. Ecol.}, 38:1082-1093, 2001.

\bibitem{NRC:2007} National Research Council (NRC).
Committee on the Status of Pollinators in North America,
2007, Status of Pollinators in North America,
The National Academies Press, Washhington, DC, 2007.

\bibitem{O:1988} Omholt, S.;
Relationships between worker longevity and intracolonial
population dynamics of the honeybee, \emph{J. Theor. Biol.}, 130:275-284, 1988

\bibitem{S:2007} Shapley, D.;
 ``Famous Einstein Bee Quote Is Bogus",
\emph{The Daily Green}, June 22, 2007.

\bibitem{SM:2004} Sumpter, D.J.T.;  Martin, S.J.;
The dynamics of virus epidemics in \emph{Varroa}-infested honey
bee colonies, \emph{J. Animal Ecology} {bf 73}:51-63, 2004.

\bibitem{T:2008} Tautz, J.;
\emph{The Buzz about Bees. Biology of a superorganism}, Springer, 2008.

\bibitem{W:2000} Walter, W.; \emph{Gew\"ohnliche Differentialgleichungen},
 7th ed., Springer, 2000.

\bibitem{W:1987} Winston, M. L.; \emph{Biology of the Honey Bee},
Harvard University Press,  Cambridge, MA, USA, 1987.

\end{thebibliography}

\end{document}