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

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. 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 Europeans bees.
The model proposed is shown in the diagram of Fig 1, and detailed in the system of differential equations below: In the proposed model, I, I i , A and A i 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 π, δ 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]. μ is the mortality rate for adult bees, γ is the mortality rate induced by the presence of mites in the colony bees. The value used for γ (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. 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 i , 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 i , H-grooming, proper hygienic behavior and harmful hygienic behavior-were studied with respect to the existence of a stable infestation equilibrium.

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 nextgeneration matrix. Let x i , i = 1, 2, . . ., m be the number or proportion of individuals in the ith compartment. Then where F i ðxÞ is the rate of appearance of new infestations in compartment i and is the rate of transfer of individuals out of the ith compartment, and 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 F i and V i .
where x 0 is the disease free equilibrium. In our model, m = 2 and the infested compartments are: Now we write the matrices F and V, substituting the mite-free equilibrium values, A Ã ¼ dp mðdþHÞ and A Ã i ¼ 0.
Let the next-generation matrix G be the matrix product FV −1 . Then Now we can find the basic reproduction number, R 0 , which is the largest eigenvalue of the matrix G.

Well-Posed and Boundedness
For sake of simplicity, we denote in such a way that the System (1) rewrites  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. We assume that all the coefficients presented in Table 1 are all positive, that is: The previous system of equations is written 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 D the subset of those elements X ¼ ðI; A; I i ; A i Þ 2 R 4 þ such that A + A i 6 ¼ 0.  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. doi:10.1371/journal.pone.0160465.g003 Theorem 1 (Well-posedness and boundedness) If X 0 2 D, then there exists a unique solution of Eq (7) defined on [0, +1) such that X(0) = X 0 . Moreover, for any t > 0, XðtÞ 2 D, and ðAðtÞ þ A i ðtÞÞ dp ma min ð8bÞ Using the values for parameters π, δ, μ, H and γ from Table 1 The red region represent 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 R 0 < 1 which means that for these these combination of g and H i the mites will be eliminated. where by definition α min ¼ : AðtÞ ð9Þ and Define D 0 as the largest set included in D and fulfilling the inequalities of Theorem 1, that is: Theorem 1 shows that the compact set D 0 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 D 0 . In Theorem 1, the notations lim inf and lim sup 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, +1) is equal to inf t ! 0 sup τ ! t f (τ). It is the largest accumulation point of f at infinity.

Equilibria
Theorem 2 (Equilibria and asymptotic behavior) Define

b¼
: • If β 0, then there exists a unique equilibrium point of System (7) in D 0 , that corresponds to a mite-free situation. It is globally asymptotically stable, and given by • If b > 1 a i , then there exists two equilibrium points in D 0 , namely X MF and a infestation equilibrium defined by Moreover, for all initial conditions in D 0 except in a zero measure set, the trajectories tend towards X CO .
The point R 0 ¼ 1, that is β = 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 b > 1 a i , and conjecture that the same property holds for β in the interval 0; 1 . Using α as bifurcation parameter, the bifurcation appears for a ¼ a i ðm i þgÞ m % 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 α 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

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 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 6 ¼ 0. As a consequence, the fact that all trajectories are defined on infinite horizon will ensue. • Summing up the first two equations in Eq (5) yields, for any point inside D: Integrating this differential inequality between any two points X(0) = X 0 and X(t) of a trajectory for which XðtÞ 2 D, τ 2 [0; t], one gets where the right-hand side is in any case positive for any t > 0. Similarly, one has and therefore IðtÞ þ I i ðtÞ p a min 1 À e Àa min t ð Þþð Ið0Þ þ I i ð0ÞÞe Àa min t : ð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 Therefore, 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, +1). Therefore, any trajectory departing from a point in D can be extended to [0, +1), 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 and thus 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 The auxiliary linear time-invariant system d dt 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 and thus Hurwitz because α(μ + g) − μα min = (α − α min )μ + αg > 0. Therefore, all trajectories of Eq (25) tend towards the unique equilibrium: 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): This gives Eq (9).
• One finally proves Eq (10). Using Eq (8b), identity Eq (5c) implies 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 where u, resp. y, is the input, resp. the output, closed by the unitary feedback 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 β 0, and there exist exactly two equilibrium points when β > 0. equilibrium points.
3. One then shows that the I/O system u 7 ! 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 Summing up the first and third identities gives and thus necessarily: • The case λ = 0 yields I = 0, and then A = 0 by Eq (31a), and therefore u has to be zero from Eq (31b). Also, I i ¼ p a i , A i ¼ dp a i ðm i þgÞ by Eq (31d), and then y ¼ gA i ¼ gdp a i ðm i þgÞ . in Eq (11) and should be discarded. obtained point is located outside D and has to be discarded; or • The case λ = 1 yields I i = 0, and then A i = 0 by Eqs (31d) or (31c), and y = 0. There remains the two following conditions: which yield (The map u 7 ! 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 λ in (0;1). From Eqs (33) and (31a)-(31c), one deduces Using Eq (33) on the one hand and summing the two identities Eqs (31b)-(31d) on the other hand, yields This permits to express A as a function of λ, namely: Using this formula together with Eqs (33), (31d) and (36) now allows to find an equation involving only the unknown λ, namely: Simplifying (as λ 6 ¼ 0, 1) gives: The previous condition is clearly affine in λ. It writes lm þ ð1 À lÞðm i þ gÞ ð Þ dp a i ¼ ðm i þ gÞ dp l a which, after developing and simplifying, can be expressed as: and finally For u ! 0, this equation admits a solution in (0;1) if and only if and the latter is given as The state and output values may then be expressed explicitly as functions of u. In particular, one has 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 λ = 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 λ = 1 yields a unique equilibrium point. Indeed, y = 0, so u should be zero too, and the unique solution is given by This corresponds to the equilibrium denoted X MF in the statement.
• Let us consider now the case of λ 2 (0;1). For this case to be considered, it is necessary that β > 0, that is R 0 > 1. The value of u should be such that (see Eq (46)) that is or again after replacing β by its value defined in Eq (12). The corresponding value of given by Eq (45), is clearly contained in (0;1) when β > 0. Therefore, when β > 0, there also exists a second equilibrium. The latter is given by: and corresponds to X CO defined in the statement. diagonal that comes from the loop closing.
3. Let K be the cone in R 4 þ defined as the product of orthants R þ Â R þ Â R À Â R À . We endow the state space with this order. In other words, for any X = (I, A, I i , A i ) and X 0 ¼ ðI 0 ; A 0 ; I 0 i ; A 0 i Þ in R 4 þ , X K X 0 means: With this structure, one may verify that the System (30a)-(30e) has the following monotonicity properties [20,21] • For any function u 2 U¼ : fu : ½0; þ1Þ ! R; locally integrable and taking on positive values almost everywhere}, for any X 0 ; X 0 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 which is irreducible when A 6 ¼ 0 and A i 6 ¼ 0. The system is therefore strongly monotone in D 0 n fX : A i ¼ 0g (notice that D 0 does not contain points for which A = 0), and also on the invariant subset D 0 \ fX : • The input-to-state map is monotone, that is: for any inputs u; u 0 2 U, for any X 0 2 R 4 þ , uðtÞ u 0 ðtÞ a:e: ) 8t ! 0; Xðt; X 0 ; uÞ K Xðt; X 0 0 ; uÞ: ð56Þ • The state-to-output map is anti-monotone, that is: for any X ; X 0 2 R 4 þ , 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 2 R þ . As shown by Theorem 1, all trajectories are precompact.
• When β 0, it has been previously established that for any u 2 R there exists at most one equilibrium in D 0 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 β > 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 This matrix is block triangular, with diagonal blocks The first of them is clearly Hurwitz, while the second, whose characteristic polynomial is (where u Ã is defined in Eq (44)) is not Hurwitz when β > 0 and 0 u 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 2 [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 2 [0;u Ã ], for almost any X 0 2 D 0 , one has lim t!þ1 Xðt; X 0 ; uÞ ¼ k X ðuÞ; lim t!þ1 yðt; X 0 ; uÞ ¼ kðuÞ; and, from the monotony properties, for any scalar-valued continuous function u, for almost any X 0 2 D 0 : 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, yðt; X 0 ; uÞÞ: ð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 dpg a i ðm i þ gÞ and it thus smaller than 1 if and only if b > 1 a i , which is an hypothesis of the statement. Under these assumptions, one then obtains that the lim inf and lim sup in Eq (63) are equal, and thus that y, and thus u, possesses limit for t ! +1. Moreover, the state itself converges towards the equilibrium X CO when t ! +1 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 behaviorwhen 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 miteinduced bee mortality, γ, must be kept low or risk de-stabilizing the colony.
Africanized Honey Bees, having evolved more effective resistance behaviors, are more resistant to colony colapse 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.