# Behavioral Modulation of Infestation by Varroa destructor in Bee Colonies. Implications for Colony Stability

Joyce de Figueiró Santos, Flávio Codeço Coelho, Pierre-Alexandre Bliman

<figure><img src="image_1.png" /></figure>

## To cite this version:

Joyce de Figueiró Santos, Flávio Codeço Coelho, Pierre-Alexandre Bliman. Behavioral Modulation of Infestation by Varroa destructor in Bee Colonies. Implications for Colony Stability. PLoS ONE, 2016, 11 (9), pp.e0160465. 10.1371/journal.pone.0160465. hal-01361716

HAL Id: hal-01361716

https://hal.sorbonne-universite.fr/hal-01361716v1

Submitted on 7 Sep 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH ARTICLE

# Behavioral Modulation of Infestation by *Varroa destructor* in Bee Colonies. Implications for Colony Stability

Joyce de Figueiró Santos¹*, Flávio Codeço Coelho¹*,*, Pierre-Alexandre Bliman¹,²*

1 Applied Mathematics School, Getulio Vargas Foundation, Rio de Janeiro, RJ, Brasil, 2 Sorbonne Universités, Inria, UPMC Univ Paris 06, Lab. J.-L. Lions UMR CNRS 7598, Paris, France

© These authors contributed equally to this work.

* fccoelho@fgv.br

<figure><img src="image_3.png" /></figure>

CrossMark

click for updates

<figure><img src="image_4.png" /></figure>

## OPEN ACCESS

**Citation:** de Figueiró Santos J, Coelho FC, Bliman P-A (2016) Behavioral Modulation of Infestation by *Varroa destructor* in Bee Colonies. Implications for Colony Stability. PLoS ONE 11(9): e0160465. doi:10.1371/journal.pone.0160465

**Editor:** Olav Rueppell, University of North Carolina at Greensboro, UNITED STATES

**Received:** September 18, 2015

**Accepted:** July 19, 2016

**Published:** September 1, 2016

**Copyright:** © 2016 de Figueiró Santos et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper.

**Funding:** The authors received no specific funding for this work.

**Competing Interests:** The authors have declared that no competing interests exist.

## Abstract

Colony Collapse Disorder (CCD) has become a global problem for beekeepers and for the crops that depend on bee pollination. While many factors are known to increase the risk of colony collapse, the ectoparasitic mite *Varroa destructor* is considered to be the most serious one. Although this mite is unlikely to cause the collapse of hives itself, it is the vector for many viral diseases which are among the likely causes for Colony Collapse Disorder. The effects of *V. destructor* infestation differ from one part of the world to another, with greater morbidity and higher colony losses in European honey bees (EHB) in Europe, Asia and North America. Although this mite has been present in Brazil for many years, there have been no reports of colony losses amongst Africanized Honey Bees (AHB). Studies carried out in Mexico have highlighted different behavioral responses by the AHB to the presence of the mite, notably as far as grooming and hygienic behavior are concerned. Could these explain why the AHB are less susceptible to Colony Collapse Disorder? In order to answer this question, we have developed a mathematical model of the infestation dynamics to analyze the role of resistance behavior by bees in the overall health of the colony, and as a consequence, its ability to face epidemiological challenges.

## Introduction

In winter and spring of 2006/2007 American beekeepers started reporting heavier and widespread losses of bee colonies and so did Europeans beekeepers. This mysterious phenomenon was called “Colony Collapse Disorder” (CCD). Diseases, parasites, in-hive chemicals, agricultural insecticides, genetically modified crops, changed cultural practices and cool brood have all been suggested as possible causes for it [1] but nowadays the ectoparasitic mite *Varroa destructor* that parasitizes honey bees is considered the most likely cause. Although *V. destructor* has become a global problem its effects vary in different parts of the world. More intense losses have been reported in European honey bee colonies (EHB) in Europe, Asia and North America [2].
The mite's life cycle is tightly linked with that of the bees. Immature mites develop with immature bees, parasitizing them from an early stage. The mite's egg-laying behavior is coupled with that of the bees and thus depends on its reproductive cycle. In the northern hemisphere bees are much less active during the cold winter months. But since worker brood rearing (and thus *Varroa* reproduction) occurs all year round in tropical climates, one would expect that the impact of the parasite would be even worse in tropical regions. But even though *V. destructor* has been present in Brazil for more than 30 years, no colony collapses due to this mite, have been recorded [3]. It is worth noting that the dominant variety of bees in Brazil is the Africanized Honey Bee (AHB) which has spread throughout the entire country since its introduction in 1956 [4]. African bees and their hybrids are known to be more resistant to the mite *V. destructor* than the European bee subspecies [4, 5]. A review by Arechavaleta-Velasco et al. [6] in Mexico showed that EHB were twice as attractive to *V. destructor* as AHB.

## Resistance behaviors of the bee against the parasite

Both varieties of bees exhibit two types of resistance to the mite: firstly, grooming where workers use their legs and mandibles to remove the mite and then injure or kill it [7], and secondly hygienic behavior where workers destroy potentially infested brood cells [8]. Grooming behavior performed by adult bees, includes detecting and eliminating mites from their own body (auto-grooming) or from the body of another bee (allo-grooming). Hygienic behavior occurs when adult bees detect the presence of mite offspring still in the cells and in order to prevent the mites from spreading in the colony, the worker bees kill the infested brood. It has been demonstrated that the smell of the mite is capable of activating this behavior [9]. Hygienic behavior serves to combat other illnesses and parasites to which the brood is susceptible but it is not 100% accurate. Correa-Marques and De Jong [9] report that the majority (53%) of the uncapped cells display no apparent signs of parasitism or other abnormality which would justify killing of the brood.

AHB workers were more efficient in grooming mites from their bodies than EHB. AHB have been shown to be more effective in hygienic behavior than EHB. Vandame et al. [7] found in Mexico that the EHB are only able to remove 8% of infested brood whereas AHB removed up 32.5%. These types of behavior are important factors in keeping mite infestation low in the honey bee colonies but they come at a cost to the bees.

Our paper is not the first to model host-parasite systems; others exist in the literature and have recently been reviewed by Becher et al. [10]. In particular, Ratti et al. [11] modelled the population dynamics of bees and mites together with the acute bee paralysis virus. Here, we focus solely on the host-parasite interactions in order to understand the resilience of colonies in Brazil and the role of the more efficient resistance behaviors displayed by AHB to explain the lower infestation rates and the lower incidence of colony collapse [7].

The main goal of this paper is to propose a model capable of describing the dynamics of infestation by *V. destructor* in bee colonies taking into consideration bee's resistance mechanisms to mite infestation, grooming and hygienic behavior. In addition, by simulating the dynamics, we show how the resistance behaviors contribute to reducing infestation levels in the colony.

## Mathematical model

Vandame et al. [12] discuss the cost-benefit of resistance mechanism of bee against mite. The grooming behavior performed by adult bees, includes detecting and eliminating mites from their own body (auto-grooming) or from the body of another bee (allo-grooming). The hygienic behavior occurs when adult bees detect the presence of the mite offspring still in the
<figure><img src="image_5.png" /><figcaption>Fig 1. Diagram to describes the dynamics of the model.</figcaption></figure>

doi:10.1371/journal.pone.0160465.g001

cells and in order to prevent the mites from spreading in the colony, the worker bees kill the infested brood. Their study compared the results for two subspecies of bees—Africanized and European—to examine whether these two mechanisms could explain the observed low compatibility between Africanized bees and the mite *Varroa destructor*, in Mexico. The results showed that *grooming* and hygienic behavior appears most intense in Africanized bees than in European bees.

The model proposed is shown in the diagram of Fig 1, and detailed in the system of differential equations below:

$$ 
\begin{aligned}
\dot{I} &= \pi \frac{A}{A + A_i} - \delta I - HI \\
\dot{A} &= \delta I + gA_i - \mu A \\
\dot{I}_i &= \pi \frac{A_i}{A + A_i} - \delta I_i - H_i I_i \\
\dot{A}_i &= \delta I_i - gA_i - (\mu + \gamma) A_i
\end{aligned}
\quad (1) $$

In the proposed model, *I*, *I*<sub>i</sub>, *A* and *A*<sub>i</sub> represent the non-infested immature bees, infested immature bees, non-infested adult worker bees and infested adult worker bees, respectively.

Daily birth rate for bees is denoted by $\pi$, $\delta$ is the maturation rate, i.e., the inverse of number of days an immature bee requires to turn in adult, this rate is the same for both infested and non-infested immature bees. The infestation of immature bees is proportional to the fraction of infested adults because females mites initiate reproduction by entering the brood cell, before it is sealed [2]. $\mu$ is the mortality rate for adult bees, $\gamma$ is the mortality rate induced by the presence of mites in the colony bees. The value used for $\gamma$ (Table 1) is insignificant, but this parameter can be used in extensions of this model to represent additional mortality due to the impact of diseases transmitted by the mite. The parameters $H_i$, $H$ and $g$ are the rate of removal of infested pupae via hygienic behavior, the general hygienic rate (affecting uninfested pupae) and grooming rate, respectively.
<table><caption>Table 1. Parameters of the model.</caption><thead><tr><td>Parameters</td><td>Meaning</td><td>Value</td><td>Unit</td><td>Reference</td></tr></thead><tbody><tr><td>π</td><td>Bee daily birth rate</td><td>2500</td><td>bees × day<sup>-1</sup></td><td>[13]</td></tr><tr><td>δ</td><td>Maturation rate</td><td>0.05</td><td>day<sup>-1</sup></td><td>[13]</td></tr><tr><td>H</td><td>Generic hygienic behavior</td><td>-</td><td>day<sup>-1</sup></td><td>-</td></tr><tr><td>H<sub>i</sub></td><td>Hygienic behavior towards infested brood</td><td>-</td><td>day<sup>-1</sup></td><td>-</td></tr><tr><td>g</td><td>Grooming</td><td>-</td><td>day<sup>-1</sup></td><td>-</td></tr><tr><td>μ</td><td>Mortality rate</td><td>0.04</td><td>day<sup>-1</sup></td><td>[14]</td></tr><tr><td>γ</td><td>Mite induced mortality</td><td>10<sup>-7</sup></td><td>day<sup>-1</sup></td><td>[11]</td></tr></tbody></table>

doi:10.1371/journal.pone.0160465.t001

## Choosing parameters values

Some of the parameters associated with the bees life cycle, used for the simulations, can be found in the literature, as shown in Table 1. For the resistance behavior parameters, *g*, *H* and *H*<sub>i</sub>, very little information is available. Therefore we decided to study the variation of these parameters within ranges which allowed for the system to switch between a mite-free equilibrium to one of stable infestation. These ranges also reflected observations described in the literature (Table 2) [6, 12, 15].

The three unknown parameters representing resistance behaviors *g*, *H*<sub>i</sub>, *H*—grooming, proper hygienic behavior and harmful hygienic behavior—were studied with respect to the existence of a stable infestation equilibrium.

## Results

### Basic reproduction number $R_0$ of the infested bees

One way of looking for a boundary beyond which infestation by mites is possible, is to compute the basic reproduction number, $R_0$ of infestation. For our model, the basic reproduction number, or $R_0$ of infested bees, can be thought of as the number of new infestations that one infested bee when introduced into the colony generates on average over the course of its infestation period or before it is groomed, in an otherwise uninfested population.

**Deriving $R_0$ using the next-generation method.** To calculate the basic reproduction number of infested bees, we will use the next-generation matrix [16], where the whole population is divided into *n* compartments in which there are *m* < *n* infested compartments. The next-generation matrix defines the instantaneous rate of expansion of the infestation, right at the start.

In this method, $R_0$ is defined as the spectral radius, or the largest eigenvalue, of the next-generation matrix.

<table><caption>Table 2. Varying the parameters.</caption><thead><tr><td>Parameter</td><td>Minimum value</td><td>Maximum value</td></tr></thead><tbody><tr><td>g</td><td>0.01</td><td>0.01</td></tr><tr><td>H<sub>i</sub></td><td>0.08</td><td>0.4</td></tr><tr><td>H</td><td>0.04</td><td>0.2</td></tr></tbody></table>

doi:10.1371/journal.pone.0160465.t002
Let $x_i, i = 1, 2, \dots, m$ be the number or proportion of individuals in the $i$th compartment. Then

$$ \frac{dx_i}{dt} = \mathcal{F}_i(x) - \mathcal{V}_i(x) $$

where $\mathcal{F}_i(x)$ is the rate of appearance of new infestations in compartment $i$ and $\mathcal{V}_i(x) = \mathcal{V}_i^-(x) - \mathcal{V}_i^+(x)$. $\mathcal{V}_i^-$ is the rate of transfer of individuals out of the $i$th compartment, and $\mathcal{V}_i^+$ represents the rate of transfer of individuals into compartment $i$ by all other means.

The next-generation matrix is then defined by $FV^{-1}$, where $F$ and $V$ can be formed by the partial derivatives of $\mathcal{F}$ and $\mathcal{V}_i$.

$$ F = \left[ \frac{\partial \mathcal{F}_i(x_0)}{\partial x_j} \right] \quad \text{and} \quad V = \left[ \frac{\partial \mathcal{V}_i(x_0)}{\partial x_j} \right] $$

where $x_0$ is the disease free equilibrium.

In our model, $m = 2$ and the infested compartments are:

$$ \begin{aligned} \frac{dI_i}{dt} &= \pi \frac{A_i}{A + A_i} - \delta I_i - H_i I_i \\ \frac{dA_i}{dt} &= \delta I_i - gA_i - (\mu + \gamma)A_i \end{aligned} \qquad (2) $$

Now we write the matrices $F$ and $V$, substituting the mite-free equilibrium values, $A^* = \frac{\delta\pi}{\mu(\delta+H)}$ and $A_i^* = 0$.

$$ F = \begin{bmatrix} 0 & \frac{\mu(\delta+H)}{\delta} \\ 0 & 0 \end{bmatrix} $$

$$ V = \begin{bmatrix} \delta + H_i & 0 \\ -\delta & g + \gamma + \mu \end{bmatrix} $$

Let the next-generation matrix $G$ be the matrix product $FV^{-1}$. Then

$$ G = \begin{bmatrix} \frac{\mu(\delta+H)}{(\delta+H_i)(g+\gamma+\mu)} & \frac{\mu(\delta+H)}{\delta(g+\gamma+\mu)} \\ 0 & 0 \end{bmatrix} $$

Now we can find the basic reproduction number, $\mathcal{R}_0$, which is the largest eigenvalue of the matrix $G$.

$$ \mathcal{R}_0 = \frac{\mu(\delta + H)}{(\delta + H_i)(g + \gamma + \mu)} \qquad (3) $$

Figs 2, 3 and 4 show the boundary between mite-free (blue region, $\mathcal{R}_0 < 1$) and infestation equilibria (red region, $\mathcal{R}_0 > 1$).

## Well-Posed and Boundedness

For sake of simplicity, we denote

$$ \alpha \doteq \delta + H, \quad \alpha_i \doteq \delta + H_i, \quad \mu_i \doteq \mu + \gamma \qquad (4) $$
<figure><img src="image_6.png" /><figcaption>Fig 2. Plot of values of $R_0$ for a range of values of g and H. $H_i = 0.01$ and remaining parameters set as described in Table 1. The region in red (top-left) corresponds to $R_0 > 1$, the black line to $R_0 = 0$ and the blue region (bottom-right) otherwise. This figure shows a slightly narrower range of the parameters as described in Table 2, for a better visualization of the threshold.</figcaption></figure>

doi:10.1371/journal.pone.0160465.g002

in such a way that the System (1) rewrites

$$ \dot{I} = \pi \frac{A}{A + A_i} - \alpha I \qquad (5a) $$

$$ \dot{A} = \delta I - \mu A + gA_i \qquad (5b) $$

$$ \dot{I}_i = \pi \frac{A_i}{A + A_i} - \alpha_i I_i \qquad (5c) $$

$$ \dot{A}_i = \delta I_i - (\mu_i + g) A_i \qquad (5d) $$
<figure><img src="image_7.png" /><figcaption>Fig 3. Values of $R_0$ for various combinations of $H_i$ and $H$. $g = 0.01$ and other parameters as given in Table 1. The region in red (bottom-right) corresponds to $R_0 > 1$, the black line to $R_0 = 0$ and the blue region (top-left) otherwise. This figure illustrates one of the conditions for infestation(given other parameters values fixed as in Table 1) that $H$ must be larger than $H_i$. This figure shows a slightly narrower range of the parameters as described in Table 2, for a better visualization of the threshold.</figcaption></figure>

doi:10.1371/journal.pone.0160465.g003

We assume that all the coefficients presented in Table 1 are all positive, that is:

$$ \pi, \delta, \mu > 0, \quad \alpha, \alpha_i > \delta, \quad \mu_i > \mu. \qquad (6) $$

The previous system of equations is written

$$ \dot{X} = f(X), \quad X = (I, A, I_i, A_i) \qquad (7) $$

The right-hand side of Eq (7) is not properly defined in the points where $A + A_i = 0$. However, the following result demonstrates that this has no consequence on the solutions, as the latter stays away from this part of the subspace. For subsequent use, we denote $\mathcal{D}$ the subset of those elements $X = (I, A, I_i, A_i) \in \mathbb{R}_+^4$ such that $A + A_i \neq 0$.
<figure><img src="image_8.png" /><figcaption>**Fig 4.** Implicit plot for $\mathcal{R}_0$ letting $g$ and $H_i$ vary. Using the values for parameters $\pi, \delta, \mu, H$ and $\gamma$ from Table 1 The red region represents $\mathcal{R}_0 > 1$ which means that for these combination of $g$ and $H_i$ the mite will stay in the colony. On the other hand, the blue region represents $\mathcal{R}_0 < 1$ which means that for these these combination of $g$ and $H_i$ the mites will be eliminated.</figcaption></figure>

doi:10.1371/journal.pone.0160465.g004

**Theorem 1** (Well-posedness and boundedness) If $X_0 \in \mathcal{D}$, then there exists a unique solution of Eq (7) defined on $[0, +\infty)$ such that $X(0) = X_0$. Moreover, for any $t > 0$, $X(t) \in \mathcal{D}$, and

$$ \frac{\pi}{\alpha_{\max}} \leq \liminf_{t \to +\infty} (I(t) + I_i(t)) \leq \limsup_{t \to +\infty} (I(t) + I_i(t)) \leq \frac{\pi}{\alpha_{\min}} \quad (8a) $$

$$ \frac{\delta\pi}{\mu_i\alpha_{\max}} \leq \liminf_{t \to +\infty} (A(t) + A_i(t)) \leq \limsup_{t \to +\infty} (A(t) + A_i(t)) \leq \frac{\delta\pi}{\mu\alpha_{\min}} \quad (8b) $$
where by definition $\alpha_{\min} \doteq \min\{\alpha; \alpha_i\}$, $\alpha_{\max} \doteq \max\{\alpha; \alpha_i\}$. Also,

$$ \frac{1}{(\alpha - \alpha_{\min})\mu + \alpha g} \frac{\pi g \mu \alpha_{\min}}{\mu_i \alpha_{\max}} \leq \liminf_{t \to +\infty} I(t), \qquad \frac{1}{(\alpha - \alpha_{\min})\mu + \alpha g} \frac{\delta \pi g \alpha}{\mu_i \alpha_{\max}} \leq \liminf_{t \to +\infty} A(t) \quad (9) $$

and

$$ (I_i(0), A_i(0)) \neq (0, 0) \Rightarrow \forall t \geq 0, I_i(t) > 0, A_i(t) > 0 \quad (10) $$

Define $\mathcal{D}'$ as the largest set included in $\mathcal{D}$ and fulfilling the inequalities of Theorem 1, that is:

$$ \mathcal{D}' = \left\{ (I, A, I_i, A_i) \in \mathbb{R}_+^4 : \begin{array}{l} \frac{\pi g \mu \alpha_{\min}}{\mu_i \alpha_{\max}} \leq I, \quad \frac{\delta \pi g \alpha}{\mu_i \alpha_{\max}} \leq A, \quad \frac{\pi}{\alpha_{\max}} \leq I + I_i \leq \frac{\pi}{\alpha_{\min}}, \\ \\ \frac{\delta \pi}{\mu_i \alpha_{\max}} \leq A + A_i \leq \frac{\delta \pi}{\mu \alpha_{\min}} \end{array} \right\}. \quad (11) $$

Theorem 1 shows that the compact set $\mathcal{D}'$ is positively invariant and attracts all the trajectories. Therefore, in order to study the asymptotics of System (5), it is sufficient to consider the trajectories of Eq (5) that are in $\mathcal{D}'$.

In Theorem 1, the notations $\liminf$ and $\limsup$ correspond respectively to the limit inferior and limit superior of a function (or lower limit and upper limit). We recall e.g. that the limit superior at infinity of a real-valued function $f$ defined on $[0, +\infty)$ is equal to $\inf_t \geq 0 \sup_{\tau \geq t} f(\tau)$. It is the largest accumulation point of $f$ at infinity.

## Equilibria

**Theorem 2 (Equilibria and asymptotic behavior) Define**

$$ \beta \doteq \frac{\mu}{\alpha_i} - \frac{\mu_i + g}{\alpha} \quad (12) $$

* If $\beta \leq 0$, then there exists a unique equilibrium point of System (7) in $\mathcal{D}'$, that corresponds to a mite-free situation. It is globally asymptotically stable, and given by

$$ X_{MF} = \frac{\pi}{\alpha} \begin{pmatrix} 1 \\ \frac{\delta}{\mu} \\ 0 \\ 0 \end{pmatrix}. \quad (13) $$

* If $\beta > \frac{1}{\alpha_i}$, then there exists two equilibrium points in $\mathcal{D}'$, namely $X_{MF}$ and a infestation equilibrium defined by

$$ X_{CO} = \frac{\delta \pi g}{\alpha_i(\mu_i + g)} \frac{\alpha \mu - \alpha_i(\mu_i + g)}{\alpha(\mu + g) - \alpha_i(\mu_i + g)} \begin{pmatrix} \frac{1}{\delta} \frac{\alpha_i(\mu_i + g)}{\alpha \mu - \alpha_i(\mu_i + g)} \\ \frac{\alpha}{\alpha \mu - \alpha_i(\mu_i + g)} \\ \frac{\mu_i + g}{\delta g} \\ \frac{1}{g} \end{pmatrix}. \quad (14) $$
<figure><img src="image_9.png" /><figcaption>Fig 5. Bifurcation diagram showing the transcritical bifurcation with bifurcation point $\alpha_0 \approx 0.125$ ($\beta = 0, R_0 = 1$). When the parameter $\alpha$ is greater than $\alpha_0$, coexistence equilibrium ($I_i > 0$) exists. When $\alpha < \alpha_0$, only the mite-free equilibrium exists. Blue dots correspond to the equilibrium values of $I_i$.</figcaption></figure>

doi:10.1371/journal.pone.0160465.g005

Moreover, for all initial conditions in $D'$ except in a zero measure set, the trajectories tend towards $X_{CO}$.

Recall that $R_0 = \frac{z_\mu}{\alpha_i(\mu_i+g)}$, in such a way that

$$ \beta > 0 \Leftrightarrow R_0 > 1. \qquad (15) $$

The point $R_0 = 1$, that is $\beta = 0$, is the point of a transcritical bifurcation, that appears when $R_0$ gets larger than 1. For larger values, two equilibria are found analytically, a mite-free one, that is unstable, and a infestation equilibrium which is stable. We've shown (Theorem 2) that the latter is globally asymptotically stable if $\beta > \frac{1}{\alpha_i}$, and conjecture that the same property holds for $\beta$ in the interval $(0, \frac{1}{\alpha_i}]$. Using $\alpha$ as bifurcation parameter, the bifurcation appears for $\alpha = \frac{z_i(\mu_i+g)}{\mu} \approx 0.125$, after substituting the parameter values (Fig 5).

If we solve numerically the system from Eq (5), we confirm the existence of two equilibria when $\alpha$ crosses the bifurcation value of 0.125. The instability and stability of the mite-free and infestation equilibria, respectively is shown in the simulation of Fig 6.
<figure><img src="image_10.png" /><figcaption>**Fig 6.** Simulation showing the infestation of a colony, by a single infested adult bee, with parameters giving $R_0 \approx 1.33$. Initial conditions: $I = 5000$, $I_i = 0$, $A = 20000$, $A_i = 0$ and parameters $g = 0.01$, $H_i = 0.1$, $\mu = 0.04$, $\delta = 0.05$, $\gamma = 10^{-7}$ and $H = 0.19$. On time $t = 100$ days, a single infested adult bee is introduced into the colony. For this simulation, $\beta = 0.375$ and $R_0 \approx 3.199$.</figcaption></figure>

doi:10.1371/journal.pone.0160465.g006

Figs 6 and 7 show simulations representing the infestation and mite-free equilibria, respectively. The time range of simulations is between 2 and 3 years, with daily time steps, which is enough for the dynamic to converge to the equilibria.

## Proofs of the theorems

**Proof of Theorem 1.** • Clearly, the right-hand side of the system of equations is globally Lipschitz on any subset of $\mathcal{D}$ where $A + A_i$ is bounded away from zero. The existence and uniqueness of the solution of System (5) is then obtained for each trajectory staying at finite distance of this boundary. We will show that the two formulas provided in the statement are valid for each trajectory departing initially from a point where $A + A_i \neq 0$. As a consequence, the fact that all trajectories are defined on infinite horizon will ensue.
<figure><img src="image_11.png" /><figcaption>Fig 7. Simulation showing the elimination of the mites from a colony, by a single infested adult bee, when R₀ < 1. Initial conditions: *I* = 15000, *I*ᵢ = 5000, *A* = 20000, *A*ᵢ = 6000 and parameters *g* = 0.01, *H*ᵢ = 0.1, *μ* = 0.05, *δ* = 0.05, *γ* = 10⁻⁷ and *H* = 0.1.</figcaption></figure>

doi:10.1371/journal.pone.0160465.g007

• Summing up the first two equations in Eq (5) yields, for any point inside D:

$$ \dot{I} + \dot{I}_i = \pi - \alpha I - \alpha_i I_i \geq \pi - \alpha_{\max}(I + I_i). \quad (16) $$

Integrating this differential inequality between any two points X(0) = X₀ and X(t) of a trajectory for which X(τ) ∈ D, τ ∈ [0; t], one gets

$$ I(t) + I_i(t) \geq \frac{\pi}{\alpha_{\max}} (1 - e^{-\alpha_{\max}t}) + (I(0) + I_i(0))e^{-\alpha_{\max}t}, \quad (17) $$

where the right-hand side is in any case positive for any *t* > 0.

Similarly, one has

$$ \dot{I} + \dot{I}_i \leq \pi - \alpha_{\min}(I + I_i), \quad (18) $$
and therefore

$$I(t) + I_i(t) \leq \frac{\pi}{\alpha_{\min}} (1 - e^{-z_{\min}t}) + (I(0) + I_i(0))e^{-z_{\min}t}. \quad (19)$$

This proves in particular that the inequalities in Eq (8a) hold for any portion of trajectory remaining inside $D$.

We now consider the evolution of $A, A_i$. Similarly to what was done for $I, I_i$, one has

$$\dot{A} + \dot{A}_i = \delta(I + I_i) - \mu A - \mu_i A_i \geq \delta(I + I_i) - \mu_i(A + A_i) \quad (20)$$

Therefore,

$$A(t) + A_i(t) \geq (A(0) + A_i(0))e^{-\mu_i t} + \delta \int_0^t (I(\tau) + I_i(\tau))e^{-\mu_i(t-\tau)}. d\tau. \quad (21)$$

Integrating the lower bound of $I + I_i$ extracted from Eq (17) yields the conclusion that any solution departing from $D$ indeed remains in $D$ as long as it is defined. On the other hand, we saw previously that trajectories remaining in $D$ could be extended on the whole semi-axis $[0, +\infty)$. Therefore, any trajectory departing from a point in $D$ can be extended to $[0, +\infty)$, and remains in $D$ for any $t > 0$. In particular, Eq (8a) holds for any trajectory departing inside $D$.

Let us now achieve the proof by bounding $A + A_i$ from above. One has

$$\dot{A} + \dot{A}_i \leq \delta(I + I_i) - \mu(A + A_i) \quad (22)$$

and thus

$$A(t) + A_i(t) \leq (A(0) + A_i(0))e^{-\mu_i t} + \delta \int_0^t (I(\tau) + I_i(\tau))e^{-\mu_i(t-\tau)}. d\tau. \quad (23)$$

Using Eq (19) then permits to achieve the proof of Eq (8b), and finally the proof of Eq (8).

* Let us now prove Eq (9). One deduces from Eqs (5a) and (5b) and the bounds established earlier the differential inequalities

$$\dot{I} \geq \frac{\pi}{\limsup(A + A_i)} A - \alpha I \geq \frac{\mu \alpha_{\min}}{\delta} A - \alpha I, \quad (24a)$$

$$\dot{A} \geq \delta I - \mu A + g(\liminf(A + A_i) - A) \geq \delta I - (\mu + g)A + \frac{\delta \pi g}{\mu_i \alpha_{\max}} \quad (24b)$$

The auxiliary linear time-invariant system

$$\frac{d}{dt} \begin{pmatrix} I' \\ A' \end{pmatrix} = \begin{pmatrix} -\alpha & \frac{\mu \alpha_{\min}}{\delta} \\ \delta & -(\mu + g) \end{pmatrix} \begin{pmatrix} I' \\ A' \end{pmatrix} + \begin{pmatrix} 0 \\ \frac{\delta \pi g}{\mu_i \alpha_{\max}} \end{pmatrix} \quad (25)$$

is monotone, as the state matrix involved is a Metzler matrix [17]. Moreover, it is asymptotically stable, as the associated characteristic polynomial is equal to

$$\left| \begin{array}{cc} s + \alpha & -\frac{\mu \alpha_{\min}}{\delta} \\ -\delta & s + \mu + g \end{array} \right| = s^2 + (\alpha + \mu + g)s + \alpha(\mu + g) - \mu \alpha_{\min}, \quad (26)$$

and thus Hurwitz because $\alpha(\mu + g) - \mu \alpha_{\min} = (\alpha - \alpha_{\min})\mu + \alpha g > 0$. Therefore, all trajectories of
Eq (25) tend towards the unique equilibrium:

$$ 
\begin{aligned}
\lim_{t \to +\infty} \begin{pmatrix} I'(t) \\ A'(t) \end{pmatrix} &= - \begin{pmatrix} -\alpha & \frac{\mu z_{\min}}{\delta} \\ \delta & -(\mu + g) \end{pmatrix}^{-1} \begin{pmatrix} 0 \\ \frac{\delta \pi g}{\mu_i z_{\max}} \end{pmatrix} \\
&= \frac{1}{(\alpha - \alpha_{\min})\mu + \alpha g} \begin{pmatrix} \mu + g & \frac{\mu z_{\min}}{\delta} \\ \delta & \alpha \end{pmatrix} \begin{pmatrix} 0 \\ \frac{\delta \pi g}{\mu_i z_{\max}} \end{pmatrix} \quad (27) \\
&= \frac{1}{(\alpha - \alpha_{\min})\mu + \alpha g} \begin{pmatrix} \frac{\pi g \mu z_{\min}}{\mu_i z_{\max}} \\ \frac{\delta \pi g \alpha}{\mu_i z_{\max}} \end{pmatrix}.
\end{aligned}
$$

Invoking Kamke's Theorem, see e.g. ([18] Theorem 10, p. 29), one deduces from Eq (24) and the monotony of Eq (25) the following comparison result, that holds for all trajectories of Eq (31):

$$ \liminf_{t \to +\infty} \begin{pmatrix} I(t) \\ A(t) \end{pmatrix} \ge \frac{1}{(\alpha - \alpha_{\min})\mu + \alpha g} \begin{pmatrix} \frac{\pi g \mu z_{\min}}{\mu_i z_{\max}} \\ \frac{\delta \pi g \alpha}{\mu_i z_{\max}} \end{pmatrix}. \quad (28) $$

This gives Eq (9).

* One finally proves Eq (10). Using Eq (8b), identity Eq (5c) implies

$$ \dot{I}_i \ge \frac{\pi}{\limsup(A+A_i)} A_i - \alpha_i I_i \ge \frac{\mu z_{\min}}{\delta} A_i - \alpha_i I_i \quad (29) $$

Joining this with Eq (5d) and using Kamke's result as before, ones deduces that both $I_i$ and $A_i$ have positive values when at least one of their two initial values are positive. This achieves the proof of Theorem 1.

**Proof of Theorem 2.** The proof is organized as follows.

1. We first write System (5) under the form of an I/O system, namely

$$ \dot{I} = \pi \frac{A}{A + A_i} - \alpha I \quad (30a) $$

$$ \dot{A} = \delta I - \mu A + u \quad (30b) $$

$$ \dot{I}_i = \pi \frac{A_i}{A + A_i} - \alpha_i I_i \quad (30c) $$

$$ \dot{A}_i = \delta I_i - (\mu_i + g) A_i \quad (30d) $$

$$ y = g A_i \quad (30e) $$

where $u$, resp. $y$, is the input, resp. the output, closed by the unitary feedback

$$ u = y. \quad (30f) $$

For subsequent use of the theory of monotone systems, one determines, for any (nonnegative) constant value of $u$, the equilibrium values of $(I, A, I_i, A_i)$ for Eqs (30a)–(30d), and the corresponding values of $y$ as given by Eq (30e).
2. The equilibrium points of System (5) are then exactly (and easily) obtained by solving the fixed point problem $u = y$ among the solutions of the previous problem. unique equilibrium points when $\beta \le 0$, and there exist exactly two equilibrium points when $\beta > 0$. equilibrium points.

3. One then shows that the I/O system $u \mapsto y$ defined by Eqs (30a)-(30e) is anti-monotone with respect to certain order relation, and the study of the stability of these equilibria shows that it admits single-valued I/S and I/O characteristics, as in [19].

4. Using this properties, the stability of the equilibria of the system obtained by closing the loop Eqs (30a)-(30e) by Eq (30f) is then established using arguments similar to Angeli and Sontag [17].

1. For fixed $u > 0$, the equilibrium equations of the I/O system (30a)-(30e) are given by

$$ \pi \frac{A}{A + A_i} - \alpha I = 0 \qquad (31a) $$

$$ \delta I - \mu A + u = 0 \qquad (31b) $$

$$ \pi \frac{A_i}{A + A_i} - \alpha_i I_i = 0 \qquad (31c) $$

$$ \delta I_i - (\mu_i + g) A_i = 0 \qquad (31d) $$

$$ y = gA_i \qquad (31e) $$

Summing up the first and third identities gives

$$ \pi = \alpha I + \alpha_i I_i, \qquad (32) $$

and thus necessarily:

$$ \exists \lambda \in [0; 1], \quad I = \lambda \frac{\pi}{\alpha}, \quad I_i = (1 - \lambda) \frac{\pi}{\alpha_i}. \qquad (33) $$

* The case $\lambda = 0$ yields $I = 0$, and then $A = 0$ by Eq (31a), and therefore $u$ has to be zero from Eq (31b). Also, $I_i = \frac{\pi}{\alpha_i}$, $A_i = \frac{\delta\pi}{\alpha_i(\mu_i + g)}$ by Eq (31d), and then $y = gA_i = \frac{g\delta\pi}{\alpha_i(\mu_i + g)}$. in Eq (11) and should be discarded. obtained point is located outside $\mathcal{D}$ and has to be discarded; or

* The case $\lambda = 1$ yields $I_i = 0$, and then $A_i = 0$ by Eqs (31d) or (31c), and $y = 0$. There remains the two following conditions:

$$ \pi = \alpha I, \quad \delta I = \mu A - u \qquad (34) $$

which yield

$$ I = \frac{\pi}{\alpha}, \quad A = \frac{\delta\pi}{\alpha\mu} + \frac{u}{\mu} \qquad (35) $$

(The map $u \mapsto y(u)$ is therefore multivalued.) Notice that these solutions do not systematically correspond to equilibrium points for the closed-loop System (30). unconditionally.
• Let us now look for possible values of $\lambda$ in (0;1). From Eqs (33) and (31a)–(31c), one deduces

$$ \frac{A}{A_i} = \frac{\alpha I}{\alpha_i I_i} = \frac{\lambda}{1 - \lambda}. \qquad (36) $$

Using Eq (33) on the one hand and summing the two identities Eqs (31b)–(31d) on the other hand, yields

$$ \delta(I + I_i) = \delta\pi\left(\frac{\lambda}{\alpha} + \frac{1-\lambda}{\alpha_i}\right) = \mu A + (\mu_i + g)A_i - u = A\left(\mu + (\mu_i + g)\frac{1-\lambda}{\lambda}\right) - u. \quad (37) $$

This permits to express A as a function of $\lambda$, namely:

$$ A = \frac{\lambda}{\lambda\mu + (1 - \lambda)(\mu_i + g)} \left[ \delta\pi\left(\frac{\lambda}{\alpha} + \frac{1 - \lambda}{\alpha_i}\right) + u \right]. \quad (38) $$

Using this formula together with Eqs (33), (31d) and (36) now allows to find an equation involving only the unknown $\lambda$, namely:

$$ \begin{aligned} \delta I_i = \frac{\delta\pi}{\alpha_i}(1 - \lambda) &= (\mu_i + g)A_i = (\mu_i + g)\frac{A_i}{A} \\ &= (\mu_i + g)\frac{1 - \lambda}{\lambda}\frac{\lambda}{\lambda\mu + (1 - \lambda)(\mu_i + g)}\left[\delta\pi\left(\frac{\lambda}{\alpha} + \frac{1 - \lambda}{\alpha_i}\right) + u\right]. \end{aligned} \quad (39) $$

Simplifying (as $\lambda \neq 0, 1$) gives:

$$ \frac{\delta\pi}{\alpha_i} = \frac{\mu_i + g}{\lambda\mu + (1 - \lambda)(\mu_i + g)} \left[ \delta\pi \left( \frac{\lambda}{\alpha} + \frac{1 - \lambda}{\alpha_i} \right) + u \right]. \quad (40) $$

The previous condition is clearly affine in $\lambda$. It writes

$$ (\lambda\mu + (1 - \lambda)(\mu_i + g))\frac{\delta\pi}{\alpha_i} = (\mu_i + g)\left(\delta\pi\left(\frac{\lambda}{\alpha} + \frac{1 - \lambda}{\alpha_i}\right) + u\right) \quad (41) $$

which, after developing and simplifying, can be expressed as:

$$ \lambda\mu\frac{\delta\pi}{\alpha_i} = (\mu_i + g)\left(\delta\pi\frac{\lambda}{\alpha} + u\right) \quad (42) $$

and finally

$$ (\mu_i + g)u = \delta\pi\left(\frac{\mu}{\alpha_i} - \frac{\mu_i + g}{\alpha}\right)\lambda = \delta\pi\beta\lambda. \quad (43) $$

For $u \ge 0$, this equation admits a solution in (0;1) if and only if

$$ \beta > 0 \quad \text{and} \quad u < u^* = \frac{\delta\pi\beta}{\mu_i + g}, \quad (44) $$

and the latter is given as

$$ \lambda = \frac{\mu_i + g}{\delta\pi\beta} u. \quad (45) $$

The state and output values may then be expressed explicitly as functions of u. In particular,
<table><caption>**Fig 8.** $\mathcal{R}_0 \le 1$ (i.e. $\beta \le 0$). See details in the text.</caption><thead><tr><th>Values of u ≥ 0</th><th>Number of distinct solutions y(u) of (31)</th></tr></thead><tbody><tr><td>u = 0</td><td>2</td></tr><tr><td>0 < u</td><td>1</td></tr></tbody></table>

doi:10.1371/journal.pone.0160465.g008

one has

$$y(u) = gA_i = \frac{\delta g}{\mu_i + g} I_i = \frac{\delta \pi g}{\alpha_i(\mu_i + g)} (1 - \lambda) = \frac{\delta \pi g}{\alpha_i(\mu_i + g)} \left(1 - \frac{\mu_i + g}{\delta \pi \beta} u\right). \quad (46)$$

value

• Eq (31) admits exactly one solution in $D'$ for any $u \ge 0$; admits a supplementary solution in $D'$ for any $u \in [0; u^*)$. Figs 8 and 9 summarize the number of solutions of Eq (31) for all non-negative values of $u$. (The map $u \mapsto y(u)$ is therefore multivalued.) Notice that these solutions do not systematically correspond to equilibrium points for the closed-loop System (30).

2. The equilibrium points of System (5) are exactly those points for which $u = y(u)$ for some nonnegative scalar $u$, where $y(u)$ is one of the output values corresponding to $u$ previously computed. We now examine in more details the solutions of this equation.

• For the value $\lambda = 0$ in the previous computations, one should have $u = 0$, due to Eq (45); but on the other hand $y > 0$ for $u = 0$, due to Eq (46). Therefore this point does not correspond to an equilibrium point of System (31).

• The value $\lambda = 1$ yields a unique equilibrium point. Indeed, $y = 0$, so $u$ should be zero too, and the unique solution is given by

$$I = \frac{\pi}{\alpha}, \quad A = \frac{\delta \pi}{\alpha \mu}, \quad I_i = 0, \quad A_i = 0, \quad y = 0. \quad (47)$$

This corresponds to the equilibrium denoted $X_{MF}$ in the statement.

• Let us consider now the case of $\lambda \in (0;1)$. For this case to be considered, it is necessary that $\beta > 0$, that is $\mathcal{R}_0 > 1$. The value of $u$ should be such that (see Eq (46))

$$y = \frac{\delta \pi g}{\alpha_i(\mu_i + g)} - \frac{g}{\alpha_i \beta} u = u, \quad (48)$$

<table><caption>**Fig 9.** $\mathcal{R}_0 > 1$ (i.e. $\beta > 0$). See details in the text.</caption><thead><tr><th>Values of u ≥ 0</th><th>Number of distinct solutions y(u) of (31)</th></tr></thead><tbody><tr><td>u = 0</td><td>3</td></tr><tr><td>0 < u < u*</td><td>2</td></tr><tr><td>u* ≤ u</td><td>1</td></tr></tbody></table>

doi:10.1371/journal.pone.0160465.g009
that is

$$ \left(1 + \frac{g}{\alpha_i \beta}\right) u = \frac{\delta \pi g}{\alpha_i (\mu_i + g)}, \qquad (49) $$

or again

$$ u = \frac{\delta \pi \beta g}{(\alpha_i \beta + g)(\mu_i + g)} = \frac{\delta \pi g}{\alpha_i (\mu_i + g)} \frac{\alpha \mu - \alpha_i (\mu_i + g)}{\alpha (\mu + g) - \alpha_i (\mu_i + g)}, \qquad (50) $$

after replacing $\beta$ by its value defined in Eq (12). The corresponding value of

$$ \lambda = \frac{\mu_i + g}{\delta \pi \beta} u = \frac{g}{\alpha_i \beta + g}, \qquad (51) $$

given by Eq (45), is clearly contained in (0;1) when $\beta > 0$. Therefore, when $\beta > 0$, there also exists a second equilibrium. The latter is given by:

$$ I = \lambda \frac{\pi}{\alpha} = \frac{\mu_i + g}{\alpha \delta \beta} u = \frac{1}{\delta} \frac{\alpha_i (\mu_i + g)}{\alpha \mu - \alpha_i (\mu_i + g)} u, \quad A_i = \frac{u}{g}, \qquad (52a) $$

$$ I_i = \frac{\mu_i + g}{\delta} A_i = \frac{\mu_i + g}{\delta g} u \qquad (52b) $$

$$ A = \frac{1}{\mu}(\delta I + u) = \frac{1}{\mu}\left(\frac{\alpha_i(\mu_i + g)}{\alpha\mu - \alpha_i(\mu_i + g)} + 1\right)u = \frac{\alpha}{\alpha\mu - \alpha_i(\mu_i + g)}u, \qquad (52c) $$

and corresponds to $X_{CO}$ defined in the statement.

diagonal that comes from the loop closing.

3. Let $\mathcal{K}$ be the cone in $\mathbb{R}_+^4$ defined as the product of orthants $\mathbb{R}_+ \times \mathbb{R}_+ \times \mathbb{R}_- \times \mathbb{R}_-$. We endow the state space with this order. In other words, for any $X = (I, A, I_i, A_i)$ and $X' = (I', A', I'_i, A'_i)$ in $\mathbb{R}_+^4$, $X \leq_{\mathcal{K}} X'$ means:

$$ I \leq I', A \leq A', I_i \geq I'_i, A_i \geq A'_i. \qquad (53) $$

With this structure, one may verify that the System (30a)-(30e) has the following monotonicity properties [20, 21]

* For any function $u \in \mathcal{U} = \{u : [0; +\infty) \to \mathbb{R}, \text{ locally integrable and taking on positive values almost everywhere}\}$, for any $X_0, X'_0 \in \mathbb{R}_+^4$,

$$ X_0 \leq_{\mathcal{K}} X'_0 \Rightarrow \forall t \geq 0, X(t; X_0, u) \leq_{\mathcal{K}} X(t; X'_0, u) \qquad (54) $$

where by definition $X(t; X_0, u)$ denotes the value at time $t$ of the point in the trajectory departing at time 0 from $X_0$ and subject to input $u$.

* The Jacobian matrix of the I/O system is

$$ \begin{pmatrix} -\alpha & \pi \frac{A_i}{(A+A_i)^2} & 0 & -\pi \frac{A}{(A+A_i)^2} \\ \delta & -\mu & 0 & 0 \\ 0 & -\pi \frac{A_i}{(A+A_i)^2} & -\alpha_i & \pi \frac{A}{(A+A_i)^2} \\ 0 & 0 & \delta & -(\mu_i + g) \end{pmatrix}, \qquad (55) $$
which is irreducible when $A \neq 0$ and $A_i \neq 0$. The system is therefore strongly monotone in $D' \setminus \{X : A_i = 0\}$ (notice that $D'$ does not contain points for which $A = 0$), and also on the invariant subset $D' \cap \{X : I_i = 0, A_i = 0, \}$.

* The input-to-state map is monotone, that is: for any inputs $u, u' \in U$, for any $X_0 \in \mathbb{R}_+^4$, 
  $u(t) \le u'(t)$ a.e. $\Rightarrow \forall t \ge 0, X(t; X_0, u) \le X(t; X'_0, u)$. (56)

* The state-to-output map is anti-monotone, that is: for any $X, X' \in \mathbb{R}_+^4$, 
  $X \le X' \Rightarrow \forall t \ge 0, gA_i \ge gA'_i$ (57)

monotone (due to the irreducibility of the Jacobian matrix) for any constant value of $u$.

* In order to construct I/S and I/O characteristics for System (31), we now examine the stability of the equilibria of System (31) for any fixed value of $u \in \mathbb{R}_+$. As shown by Theorem 1, all trajectories are precompact.

* When $\beta \le 0$, it has been previously established that for any $u \in \mathbb{R}$ there exists at most one equilibrium in $D'$ to the I/O System (31). The strong monotonicity property of this system depicted above then implies that this equilibrium is globally attractive ([20] Theorem 10.3). Therefore, System (31) possesses I/S and I/O characteristics. As for any value of $u$, this equilibrium corresponds to zero output, the I/O characteristics is zero. Applying the results of Angeli and Sontag [19], one gets that the closed-loop system equilibrium $X_{MF}$ is an almost globally attracting equilibrium for System (5).

* Let us now consider the case where $\beta > 0$. We first show that the equilibrium point with $I_i = 0, A_i = 0$ and Eq (34) is locally unstable. Notice that this point is located on a branch of solution parametrized by $u$ and departing from $X_{MF}$ for $u = 0$. The Jacobian matrix Eq (55) taken at this point is

$$ \begin{pmatrix} -\alpha & 0 & 0 & -\frac{\mu\alpha\pi}{\delta\pi+\alpha u} \\ \delta & -\mu & 0 & 0 \\ 0 & 0 & -\alpha_i & \frac{\mu\alpha\pi}{\delta\pi+\alpha u} \\ 0 & 0 & \delta & -(\mu_i+g) \end{pmatrix}. \quad (58) $$

This matrix is block triangular, with diagonal blocks

$$ \begin{pmatrix} -\alpha & 0 \\ \delta & -\mu \end{pmatrix} \quad \text{and} \quad \begin{pmatrix} -\alpha_i & \frac{\mu\alpha\pi}{\delta\pi+\alpha u} \\ \delta & -(\mu_i+g) \end{pmatrix}. \quad (59) $$

The first of them is clearly Hurwitz, while the second, whose characteristic polynomial is

$$ \begin{align} s^2 + (\alpha_i + \mu_i + g)s + \alpha_i(\mu_i + g) - \frac{\mu\alpha\delta\pi}{\delta\pi + \alpha u} &= s^2 + (\alpha_i + \mu_i + g)s - \alpha\alpha_i(\beta - u(\mu_i + g)) \tag{60} \\ &= s^2 + (\alpha_i + \mu_i + g)s - \alpha\alpha_i(\mu_i + g)(u^* - u) \end{align} $$

(where $u^*$ is defined in Eq (44)) is not Hurwitz when $\beta > 0$ and $0 \le u \le u^*$, and has a positive root for $0 < u < u^*$. Therefore, the corresponding equilibrium of the I/O System (30) is unstable for these values of $u$.
The other solution, given as a function of u by Eq (52), is located on a branch of solution parametrized by u and departing from $X_{CO}$ for $u = 0$. As the other solution is unstable for $0 < u < u^*$, one can deduce from Hirsch [20] that these solutions are locally asymptotically stable.

• One may now associate to any $u \in [0; u^*]$ the corresponding unique locally asymptotically stable equilibrium point, and the corresponding output value, defining therefore respectively an I/S characteristic $k_X$ and an I/O characteristic k for System (30).

For any scalar $u \in [0; u^*]$, for almost any $X_0 \in D'$, one has

$$ \lim_{t \to +\infty} X(t; X_0, u) = k_X(u), \quad \lim_{t \to +\infty} y(t; X_0, u) = k(u), \qquad (61) $$

and, from the monotony properties, for any scalar-valued continuous function u, for almost any $X_0 \in D'$:

$$ k(\limsup_{t \to +\infty} u(t)) \le \liminf_{t \to +\infty} y(t; X_0, u) \le \limsup_{t \to +\infty} y(t; X_0, u) \le k(\liminf_{t \to +\infty} u(t)). \qquad (62) $$

Using the fact that k is anti-monotone and that $u = y$ for the closed-loop system, one deduces, as e.g. in Gouzé [22] that, for the solutions of the latter,

$$ k^{2l}(\liminf_{t \to +\infty} y(t; X_0, u)) \le \liminf_{t \to +\infty} y(t; X_0, u) \le \limsup_{t \to +\infty} y(t; X_0, u) \le k^{2l}(\limsup_{t \to +\infty} y(t; X_0, u)). \qquad (63) $$

Here $k(u)$, defined by Eq (46), is a linear decreasing map. When its slope is smaller than 1, then the sequences in the left and right of Eq (63) tend towards the fixed point that corresponds to the output value at $X = X_{CO}$, see Eq (50).

This slope value, see Eq (46), is equal to

$$ \frac{\delta \pi g}{\alpha_i(\mu_i + g)} \frac{\mu_i + g}{\delta \pi \beta} = \frac{1}{\alpha_i \beta}, \qquad (64) $$

and it thus smaller than 1 if and only if $\beta > \frac{1}{\alpha_i}$, which is an hypothesis of the statement.

Under these assumptions, one then obtains that the $\liminf$ and $\limsup$ in Eq (63) are equal, and thus that $y$, and thus $u$, possesses limit for $t \to +\infty$. Moreover, the state itself converges towards the equilibrium $X_{CO}$ when $t \to +\infty$ for almost every initial conditions $X(0)$. This achieves the proof of Theorem 2.

## Discussion

The parasitism of bees by Varroa mites in nature is an undeniable fact. However, this parasitic relationship is fraught with dangers for the bees, since Varroa mites can be vectors of lethal viral diseases. These deleterious effects for the health of the individual workers and the whole colony, has led to the evolution of resistance behaviors such as the hygienic behavior and grooming.

Those behaviors are not entirely without cost to the bees, exacerbated hygienic behavior— when both $H$ and $H_i$ are intensified—can exert a substantial toll on the fitness of the queen. So it is safe to say that this parasitic relationship has evolved within a very narrow range of parameters. Even if the mite-free equilibrium is advantageous to the colony, maintaining it may be too expensive to the bees.

On the other hand, in the absence of viral diseases, mite parasitism seems to be fairly harmless. If we look at the expression for the $R_0$ of infestation Eq (3), we can see that the mite-induced bee mortality, $\gamma$, must be kept low or risk de-stabilizing the colony.

Africanized Honey Bees, having evolved more effective resistance behaviors, are more resistant to colony collapse through this ability to keep infestation levels lower when compared to
their European counterparts [23, 24]. Unfortunately, the lack of more detailed experiments measuring the rates of grooming and higienic behaviors in both groups (EHB and AHB), makes it hard to position them accurately in the parameter space of the model presented.

In this model, we chose to leave seasonal effects out, for simplicity, even though it is known that colonies in temperate climates suffer substantial losses during the winter. Such effects can be added to this model through the use of a time-varying mortality and birth rates. Nevertheless, we are convinced this simple model still applies to tropical colonies, and our observations about infestation levels and colony vulnerability remain relevant regardless of external morbidity factors such as hard winters.

Finally, we hope that the model presented here along with its demonstrated dynamical properties will serve as a solid foundation for the development of other models including viral dynamics and other aspects of bee colony health.

## Acknowledgments

The authors are grateful for valuable comments by Moacyr A. H. Silva, Max O. Souza and Jair Koiler on an early version of the manuscript.

## Author Contributions

**Conceived and designed the experiments:** JFS FCC PAB.

**Performed the experiments:** JFS FCC PAB.

**Analyzed the data:** JFS FCC PAB.

**Contributed reagents/materials/analysis tools:** JFS FCC PAB.

**Wrote the paper:** JFS FCC PAB.

## References

1. Oldroyd BP. What's Killing American Honey Bees? PLoS Biol. 2007; 5(6):e168. doi: 10.1371/journal.pbio.0050168 PMID: 17564497

2. Calderón RA, Veen JWv, Sommeijer MJ, Sanchez LA. Reproductive biology of Varroa destructor in Africanized honey bees (Apis mellifera). Experimental and Applied Acarology. 2010; 50(4):281–297. doi: 10.1007/s10493-009-9325-4 PMID: 19851876

3. Carneiro FE, Torres RR, Strapazzon R, Ramírez SA, Guerra JCV Jr, Koling DF, et al. Changes in the reproductive ability of the mite Varroa destructor (Anderson e Trueman) in africanized honey bees (Apis mellifera L.) (Hymenoptera: Apidae) colonies in southern Brazil. Neotropical Entomology. 2007; 36(6):949–952. doi: 10.1590/S1519-566X2007000600018 PMID: 18246271

4. Pinto FA, Puker A, Barreto LMRC, Message D. The ectoparasite mite Varroa destructor Anderson and Trueman in southeastern Brazil apiaries: effects of the hygienic behavior of Africanized honey bees on infestation rates. Arquivo Brasileiro de Medicina Veterinária e Zootecnia. 2012; 64(5):1194–1199. doi: 10.1590/S0102-09352012000500017

5. Medina LM, Martin SJ. A Comparative Study of Varroa Jacobsoni Reproduction in Worker Cells of Honey Bees (Apis Mellifera) in England and Africanized Bees in Yucatan, Mexico. Experimental & Applied Acarology. 1999; 23(8):659–667. doi: 10.1023/A:1006275525463

6. Arechavaleta-Velasco ME, Guzman-Novoa E. Relative effect of four characteristics that restrain the population growth of the mite *Varroa destructor* in honey bee (*Apis mellifera*) colonies. Apidologie. 2001; 32(2):157–174. doi: 10.1051/apido:2001121

7. Vandame R, Colin ME, Morand S, Otero-Colina G. Levels of compatibility in a new host-parasite association: Apis mellifera/Varroa jacobsoni. Canadian Journal of Zoology. 2000; 78(11):2037–2044. doi: 10.1139/z00-109

8. Spivak M. Honey bee hygienic behavior and defense against Varroa jacobsoni. Apidologie. 1996; 27:245–260. doi: 10.1051/apido:19960407
9. Corrêa-Marques MH, De Jong D. Uncapping of worker bee brood, a component of the hygienic behavior of Africanized honey bees against the mite Varroa jacobsoni Oudemans. Apidologie. 1998; 29 (3):283-289. doi: 10.1051/apido:19980307

10. Becher MA, Osborne JL, Thorbek P, Kennedy PJ, Grimm V. Review: towards a systems approach for understanding honeybee decline: a stocktaking and synthesis of existing models. Journal of Applied Ecology. 2013; 50(4):868-880. doi: 10.1111/1365-2664.12112 PMID: 24223431

11. Ratti V, Kevan PG, Eberl HJ. A mathematical model for population dynamics in honeybee colonies infested with Varroa destructor and the Acute Bee Paralysis Virus. Canadian Applied Mathematics Quarterly. 2012;

12. Vandame R, Morand S, Colin ME, Belzunces LP. Parasitism in the social bee Apis mellifera: quantifying costs and benefits of behavioral resistance to Varroa destructor mites. Apidologie. 2002; 33(5):433-445. doi: 10.1051/apido:2002025

13. Pereira FdM, Lopes MTR, Camargo RCR, Vilela SLO. Organização Social e Desenvolvimento das abelhas Apis mellifera; 2002. Available from: http://sistemasdeproducao.cnptia.embrapa.br/FontesHTML/Mel/SPMel/organizacao.htm

14. Khoury DS, Myerscough MR, Barron AB. A Quantitative Model of Honey Bee Colony Population Dynamics. PLoS ONE. 2011; 6(4):e18491. doi: 10.1371/journal.pone.0018491 PMID: 21533156

15. Mondragón L, Spivak M, Vandame R. A multifactorial study of the resistance of honeybees Apis mellifera to the mite Varroa destructor over one year in Mexico. Apidologie. 2005; 36(3):345-358. doi: 10.1051/apido:2005022

16. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences. 2002; 180(1):29-48. doi: 10.1016/S0025-5564(02)00108-6 PMID: 12387915

17. Angeli D, Sontag ED. Monotone control systems. Automatic Control, IEEE Transactions on. 2003; 48 (10):1684-1698. doi: 10.1109/TAC.2003.817920

18. Coppel WA. Stability and asymptotic behavior of differential equations. vol. 11. Heath Boston; 1965.

19. Angeli D, Sontag E. Interconnections of monotone systems with steady-state characteristics. In: Optimal control, stabilization and nonsmooth analysis. Springer; 2004. p. 135-154.

20. Hirsch MW. Stability and convergence in strongly monotone dynamical systems. J reine angew Math. 1988; 383(1):53.

21. Smith HL. Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems. vol. 41. American Mathematical Soc.; 2008.

22. Gouzé JL. A criterion of global convergence to equilibrium for differential systems. Application to Lotka-Volterra systems; 1988. RR-0894. Available from: https://hal.inria.fr/inria-00075661

23. Moretto G, Gonçalves LS, De Jong D, Bichuette MZ, others. The effects of climate and bee race on Varroa jacobsoni Oud infestations in Brazil. Apidologie. 1991; 22(3):197-203. doi: 10.1051/apido:19910303

24. Moretto G, Gonçalves LS, De Jong D. Heritability of Africanized and European honey bee defensive behavior against the mite Varroa jacobsoni. Revista Brasileira de Genetica. 1993; 16:71-71.