http://orcid.org/0000-0003-3868-4391
Figures
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.
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. https://doi.org/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, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: 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.
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 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:
(1)
In the proposed model, I, Ii, A and Ai 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 Hi, 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 Hi, 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, Hi, 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 
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, 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.
Deriving 
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, is defined as the spectral radius, or the largest eigenvalue, of the next-generation matrix.
Let xi, 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.
The next-generation matrix is then defined by FV−1, where F and V can be formed by the partial derivatives of and
.
where x0 is the disease free equilibrium.
In our model, m = 2 and the infested compartments are:
(2)
Now we write the matrices F and V, substituting the mite-free equilibrium values, and
.
Let the next-generation matrix G be the matrix product FV−1. Then
Now we can find the basic reproduction number, , which is the largest eigenvalue of the matrix G.
(3)
Figs 2, 3 and 4 show the boundary between mite-free (blue region, ) and infestation equilibria (red region,
).
Hi = 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.
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 Hi. This figure shows a slightly narrower range of the parameters as described in Table 2, for a better visualization of the threshold.
Using the values for parameters π, δ, μ, H and γ from Table 1 The red region represent which means that for these combination of g and Hi 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 Hi the mites will be eliminated.
Well-Posed and Boundedness
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 + Ai = 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 + Ai ≠ 0.
Theorem 1 (Well-posedness and boundedness) If , then there exists a unique solution of Eq (7) defined on [0, +∞) such that X(0) = X0. Moreover, for any t > 0,
, and
(8a)
(8b)
where by definition αmin ≐ min{α; αi}, αmax ≐ max{α; αi}. Also,
(9)
and
(10)
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
.
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 inft ≥ 0supτ ≥ t f(τ). It is the largest accumulation point of f at infinity.
Equilibria
Theorem 2 (Equilibria and asymptotic behavior) Define
(12)
• 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)
• If , then there exists two equilibrium points in
, namely XMF and a infestation equilibrium defined by
(14)
Moreover, for all initial conditions in
except in a zero measure set, the trajectories tend towards XCO.
Recall that , in such a way that
(15)
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).
When the parameter α is greater than α0, coexistence equilibrium (Ii > 0) exists. When α < α0, only the mite-free equilibrium exists. Blue dots correspond to the equilibrium values of Ii.
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.
Initial conditions: I = 5000, Ii = 0, A = 20000, Ai = 0 and parameters g = 0.01, Hi = 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 .
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 = 15000, Ii = 5000, A = 20000, Ai = 6000 and parameters g = 0.01, Hi = 0.1, μ = 0.05, δ = 0.05, γ = 10−7 and H = 0.1.
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 + Ai 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 + Ai ≠ 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) = X0 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, Ai. Similarly to what was done for I, Ii, one has
(20)
Therefore,
(21)
Integrating the lower bound of I + Ii 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 + Ai 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 Ii and Ai 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) where u, resp. y, is the input, resp. the output, closed by the unitary feedback
(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, Ii, Ai) 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 Ii = 0, and then Ai = 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).
See details in the text.
See details in the text.
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 XMF 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 XCO 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, Ii, Ai) 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
locally integrable and taking on positive values almost everywhere}, for any
,
(54) where by definition X(t; X0, u) denotes the value at time t of the point in the trajectory departing at time 0 from X0 and subject to input u.
- The Jacobian matrix of the I/O system is
(55) which is irreducible when A ≠ 0 and Ai ≠ 0. The system is therefore strongly monotone in
(notice that
does not contain points for which A = 0), and also on the invariant subset
.
- The input-to-state map is monotone, that is: for any inputs
, for any
,
(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 XMF 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 Ii = 0, Ai = 0 and Eq (34) is locally unstable. Notice that this point is located on a branch of solution parametrized by u and departing from XMF 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 XCO 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 kX 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 = XCO, 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.
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 → +∞. Moreover, the state itself converges towards the equilibrium XCO when t → +∞ 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 Hi 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 of infestation Eq (3), we can see that the mite-induced 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.
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. 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. 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. 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.
- 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.
- 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.
- 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.
- 8. Spivak M. Honey bee hygienic behavior and defense against Varroa jacobsoni. Apidologie. 1996;27:245–260.
- 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.
- 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. 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.
- 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. 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.
- 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. pmid:12387915
- 17. Angeli D, Sontag ED. Monotone control systems. Automatic Control, IEEE Transactions on. 2003;48(10):1684–1698.
- 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.
- 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.