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.

Copyright: © 2016 de Figueiró Santos et al. This is an open access article distributed under the terms of the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

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.

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

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.

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.

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.

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

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

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.

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

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

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.

Initial conditions: I = 5000, I i = 0, A = 20000, A i = 0 and parameters g = 0.01, H i = 0.1, μ = 0.04, δ = 0.05, γ = 10 −7 and H = 0.19. On time t = 100 days, a single infested adult bee is introduced into the colony. For this simulation, β = 0.375 and .

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 Fig 6 .

When the parameter α is greater than α 0 , coexistence equilibrium (I i > 0) exists. When α < α 0 , only the mite-free equilibrium exists. Blue dots correspond to the equilibrium values of I i .

The point , that is β = 0, is the point of a transcritical bifurcation, that appears when 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 , and conjecture that the same property holds for β in the interval . Using α as bifurcation parameter, the bifurcation appears for , after substituting the parameter values ( Fig 5 ).

• If , then there exists two equilibrium points in , namely X MF and a infestation equilibrium defined by (14) Moreover, for all initial conditions in except in a zero measure set, the trajectories tend towards X CO .

• If β ≤ 0, then there exists a unique equilibrium point of System (7) in , that corresponds to a mite-free situation. It is globally asymptotically stable, and given by (13)

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, +∞) is equal to inf t ≥ 0 sup τ ≥ t f(τ). It is the largest accumulation point of f at infinity.

Define as the largest set included in and fulfilling the inequalities of Theorem 1, that is: (11) Theorem 1 shows that the compact set 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 .

Theorem 1 (Well-posedness and boundedness) If , then there exists a unique solution of Eq (7) defined on [0, +∞) such that X(0) = X 0 . Moreover, for any t > 0, , and (8a) (8b) where by definition α min ≐ min{α; α i }, α max ≐ max{α; α i }. Also, (9) and (10)

For sake of simplicity, we denote (4) in such a way that the System (1) rewrites (5a) (5b) (5c) (5d) We assume that all the coefficients presented in Table 1 are all positive, that is: (6) The previous system of equations is written (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 the subset of those elements such that A + A i ≠ 0.

Using the values for parameters π, δ, μ, H and γ from Table 1 The red region represent 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 which means that for these these combination of g and H i the mites will be eliminated.

g = 0.01 and other parameters as given in Table 1 . The region in red (bottom-right) corresponds to , the black line to 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.

H i = 0.01 and remaining parameters set as described in Table 1 . The region in red (top-left) corresponds to , the black line to 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.

Now we can find the basic reproduction number, , which is the largest eigenvalue of the matrix G. (3)

The next-generation matrix is then defined by FV −1 , where F and V can be formed by the partial derivatives of and . where x 0 is the disease free equilibrium.

Let x i , i = 1, 2, …, m be the number or proportion of individuals in the ith compartment. Then where 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 represents the rate of transfer of individuals into compartment i by all other means.

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

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.

One way of looking for a boundary beyond which infestation by mites is possible, is to compute the basic reproduction number, of infestation. For our model, the basic reproduction number, or 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.

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 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 ≠ 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 : (16) Integrating this differential inequality between any two points X(0) = X 0 and X(t) of a trajectory for which , τ ∈ [0; t], one gets (17) where the right-hand side is in any case positive for any t > 0.

Similarly, one has (18) and therefore (19) This proves in particular that the inequalities in Eq (8a) hold for any portion of trajectory remaining inside .

We now consider the evolution of A, A i . Similarly to what was done for I, I i , one has (20) Therefore, (21) Integrating the lower bound of I + I i extracted from Eq (17) yields the conclusion that any solution departing from indeed remains in as long as it is defined. On the other hand, we saw previously that trajectories remaining in could be extended on the whole semi-axis [0, +∞). Therefore, any trajectory departing from a point in can be extended to [0, +∞), and remains in for any t > 0. In particular, Eq (8a) holds for any trajectory departing inside .

Let us now achieve the proof by bounding A + A i from above. One has (22) and thus (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 (24a) (24b) The auxiliary linear time-invariant system (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 (26) and thus Hurwitz because α(μ + g) − μα min = (α − α min )μ + αg > 0. Therefore, all trajectories of Eq (25) tend towards the unique equilibrium: (27) 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): (28) This gives Eq (9).

• One finally proves Eq (10). Using Eq (8b), identity Eq (5c) implies (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.

We first write System (5) under the form of an I/O system, namely (30a) (30b) (30c) (30d) (30e) (30f) i , A i ) for Eqs (30a)–(30d), and the corresponding values of y as given by Eq (30e). 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. One then shows that the I/O system u ↦ 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]. 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 (31a) (31b) (31c) (31d) (31e) Summing up the first and third identities gives (32) and thus necessarily: (33)

• The case λ = 0 yields I = 0, and then A = 0 by Eq (31a), and therefore u has to be zero from Eq (31b). Also, , by Eq (31d), and then . in Eq (11) and should be discarded. obtained point is located outside 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: (34) which yield (35) (The map u ↦ 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 (36) Using Eq (33) on the one hand and summing the two identities Eqs (31b)–(31d) on the other hand, yields (37) This permits to express A as a function of λ, namely: (38) Using this formula together with Eqs (33), (31d) and (36) now allows to find an equation involving only the unknown λ, namely: (39) Simplifying (as λ ≠ 0, 1) gives: (40) The previous condition is clearly affine in λ. It writes (41) which, after developing and simplifying, can be expressed as: (42) and finally (43) For u ≥ 0, this equation admits a solution in (0;1) if and only if (44) and the latter is given as (45) The state and output values may then be expressed explicitly as functions of u. In particular, one has (46) value

•Eq (31) admits exactly one solution in for any u ≥ 0; admits a supplementary solution in for any u ∈ [0; u*). Figs 8 and 9 summarize the number of solutions of Eq (31) for all nonnegative values of u. (The map u ↦ 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 λ = 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 (47) This corresponds to the equilibrium denoted X MF in the statement.

• Let us consider now the case of λ ∈ (0;1). For this case to be considered, it is necessary that β > 0, that is . The value of u should be such that (see Eq (46)) (48) that is (49) or again (50) after replacing β by its value defined in Eq (12). The corresponding value of (51) 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: (52a) (52b) (52c) and corresponds to X CO defined in the statement.

diagonal that comes from the loop closing.

3. Let be the cone in defined as the product of orthants . We endow the state space with this order. In other words, for any X = (I, A, I i , A i ) and in , means: (53) With this structure, one may verify that the System (30a)–(30e) has the following monotonicity properties [20, 21]

For any function (54) 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.

, u) denotes the value at time t of the point in the trajectory departing at time 0 from X and subject to input u. The Jacobian matrix of the I/O system is (55) i ≠ 0. The system is therefore strongly monotone in

≠ 0. The system is therefore strongly monotone in The input-to-state map is monotone, that is: for any inputs (56)

The state-to-output map is anti-monotone, that is: for any (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 . As shown by Theorem 1, all trajectories are precompact.

• When β ≤ 0, it has been previously established that for any there exists at most one equilibrium in 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 (58) This matrix is block triangular, with diagonal blocks (59) The first of them is clearly Hurwitz, while the second, whose characteristic polynomial is (60) (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 ∈ [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 ∈ [0; u*], for almost any , one has (61) and, from the monotony properties, for any scalar-valued continuous function u, for almost any : (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, (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 (64) and it thus smaller than 1 if and only if , which is an hypothesis of the statement.