Modeling the Influence of Mites on Honey Bee Populations

The Varroa destructor mite has been associated with the recent decline in honey bee populations. While experimental data are crucial in understanding declines, insights can be gained from models of honey bee populations. We add the influence of the V. destructor mite to our existing honey bee model in order to better understand the impact of mites on honey bee colonies. Our model is based on differential equations which track the number of bees in each day in the life of the bee and accounts for differences in the survival rates of different bee castes. The model shows that colony survival is sensitive to the hive grooming rate and reproductive rate of mites, which is enhanced in drone capped cells.


Introduction
The decline of honey bees (Apis mellifera) has been associated with many factors, which include the Varroa destructor and Varroa jacobsoni mite and the viruses that they transmit, the microsporidia fungus Nosema ceranae, pesticides, including neonicotinoids, and bee management practices [1][2][3]. Here, we study the effect of the V. destructor mite with a mathematical model. The life cycle of the V. destructor mite is composed of a phoretic phase and a reproductive phase. In the phoretic phase, mites remain attached to adult bees where they feed on adult bee fat body tissue [4]. In the reproductive phase, the foundress female mite enters a worker brood or drone cell just before it is capped and lays eggs. Mites increase the mortality rate of bees, especially if they transmit a virus (e.g., the deformed wing virus and the acute paralysis virus) [5].
Models of honey bees have been developed by several authors [6][7][8][9][10][11][12]. Other authors have developed models to describe the growth of mites within honey bee colonies. Calis et al. [13] extends the model of Fries et al. [14] to study the growth of V. jacobsoni in colonies. Mites use available brood cells for reproduction but do not affect the bee population. Mite invasion rates are based on data from Boot et al. [15]. The model accounts for differences in mite mortality in the summer and winter as well as mite fertility in drone and worker cells. Wilkinson and Smith [16] use difference equations for adult bees and phoretic mites to study the growth rate of mites and perform sensitivity studies to study the effect of seasonal variation, the amount of drone brood, and post-capping times. Their model bases the growth of adult bees on the egg laying rate of the queen and similar to Calis et al. [13] accounts for differences in the invasion and reproductive rates of mites in worker and drone brood.
Many honey bee models include the effects of the V. destructor mite and the viruses that they transmit on the honey bee colony. DeGrandi-Hoffman and Curry [17] create a validated mathematical model that includes the effects of miticides. Martin [18] constructs a model of viral infection which uses mites as vectors. The model modifies the model developed by DeGrandi-Hoffman [8] and integrates meteorological conditions. Survival rates used are based on data from Fukuda and Sakagami [19]. Sumpter and Martin [20] create a model based on differential equations which tracks healthy hive bees, hive bees which acquire a virus, mites, and virus carrying mites. Hive bees are defined to be the adult worker bees that live inside the hive. Kang et al. [21] propose a model that includes parasitism, virus transmission terms, and allee effects. Ratti et al. [22] assemble a model of four differential equations, which tracks healthy and virus infested mites, and healthy and virus infested worker bees. Dénes and Ibrahim [23] use differential equations to model forager bees and three compartments of hive bees (susceptible, infested with virus free mites, and infested with virus-infected mites).
Torres et al. [24] develop a model based on differential equations that tracks each day in the life of the bee and uses different survival rates for each of the different bee castes. Th survival rates are affected by the relative number of bees in each caste. For example, the brood survival rate is decreased if there is insufficient hive bees. The model is unique in that it accounts for the brood and ethyl oleate pheromone by slowing or accelerating the maturation of hive bees into foraging bees. In this article, we add the V. destructor species to our existing model. This modification entails adding differential equations for the mite population, the drone population, and classes of infested pupae, drones, hive, and forager bees.

Mathematical Model
We begin with the fundamental Equation (1) from our transient model [24] where i refers to the age of the bee in days, S i is the daily survival rate of a bee that is i days old, and B i is the number of bees that are i days old. The term a i accelerates or decelerates the maturation of hive bees into foraging bees due to the presence of pheromones. Equation (1) is solved using a forward Euler difference equation, B n+1 i = B n i + tS n i−1 B n i−1 a n i − tB n i a n i (2) where the subscript n refers to the time n t and t is the time step. While there are more accurate schemes, (2) benefits from simple conservation properties. The term − tB n i accounts for the loss of tB n i bees from the number of bees that are i days old due to natural aging in time t. Similarly, the term tS n i−1 B n i−1 accounts for the movement of tB n i−1 bees that are i − 1 days old into the number of bees that are i days old. However, the survival factor S n i−1 reduces the number of bees tB n i−1 that become i days old. Table 1 shows how an entire bee lifespan is divided into castes based on Schmickl et al. [12] and defines the total number of days a bee spends as an egg E t (i), brood or larvae L t (i), pupae P t (i), hive bee H t (i), and foraging bee F t (i). Each different caste has a different associated survival rate S i .
We also add the castes for drones B D i , as shown in Table 2. The days spent in each drone caste are reported by DeGrandi-Hoffman and Curry [17] and used in the table.
Adult Drone (D) 21 Table 3 shows the daily mortality (m) and survival rates (S = 1 − m) based on data from Schmickl et al. [12] who base their rates on experimental data from Sakagami and Fukuda [25]. Hive bees also have a higher survival rate in the winter months, which we assume to be 0.9947 using the winter mortality rate that was provided by Sumpter and Martin [20]. If there are insufficient hive bees to tend the larvae, the survival rate of the larvae S reduce i , 4 ≤ i ≤ 8 is reduced according to the equation where S i = 0.99 is the normal survival rate of larvae provided in Table 3, R H L is the ratio of hive bees to larvae and (R H L ) ideal is the ideal ratio, which is assumed to be 2 [12]. The parameter α = 1 5 prevents the survival rate from declining too steeply if r is less than but close to 1. Martin [18] also modifies the survival rate of egg, larva, and pupae, but bases the change in survival on the total number of hive bees. The summer season is determined by the beginning S begin and end days S end . In the summer, hive bees mature into foraging bees, while in the winter they remain hive bees and survive at a higher daily survival rate of 0.9947. The daily egg laying rate is determined by the parameters S b < S egg < S begin < S peak < S end < S e . We found that the equation approximately matches the growth rate of the egg laying curve of Schmickl and Crailsheim [12] and can be used with adjustable beginning and ending summer dates. The parameters S b , S egg , and S e allow the queen bee to begin and end the egg laying rate at a non-zero value in early spring and late summer. S peak is the day when the egg laying rate peaks for the summer. L max represents the maximum daily egg laying rate during the summer season.

Adding Mites to the Mathematical Model
In order to account for mites, we add equations to model the aging and survival of mites within the colony and assume mites survive an average of 27 days [21] in the summer when larvae brood is present Here, M i represents the number of mites that are i days old, and S M i−1 refers to the daily survival rate of mites. The survival rate of the mites is increased during the winter months. Martin provides the daily mortality rates of mites [26], which are 0.006 (corresponding to a survival rate of S M = 0.994) in the summer and 0.002 (corresponding to a survival rate of 0.998) in the winter.
Phoretic mites are assumed to be equally distributed among the infested hive and drone bees. We assume that, if a bee is groomed, the mites on the bee die. Mites also die in a capped cell if the pupae dies within the cell. These reductions in the mite population are incorporated in the term D i in (6). The product of the number of infested hive and drone bees times a maximum mite per bee value ξ determines the maximum number of phoretic mites that can be sustained by the colony.

Adding Drones and Infected Bee Populations to the Mathematical Model
An infested (denoted by an asterisk *) hive or forager bee is defined to be a bee with attached mites. In the case of pupae, an infested pupae harbors mites within its capped cell. We add additional equations for pupae infested with mites (8), and noninfested or "healthy" drone larvae and pupae B D i , 4 ≤ i ≤ 10, 11 ≤ i ≤ 24 (9), Infested pupae and foragers : Infested drone pupae : Healthy drone larvae and pupae : For the healthy drone castes, we assume the same survival rates as their hive bee counterparts (see Table 3). However, we use different infested survival rates S * i = 1 − m * i based on an increased mortality rate m * i for the infested castes. The equations pertaining to noninfested or healthy B i bees (10) and infested hive bees B * i (11), 21 ≤ i ≤ 41, and healthy B D i (12), and infested drone bees B D * i (13), 25 ≤ i ≤ 45, require additional terms to account for the movement of healthy bees B i , B D i into the subpopulation of infested bees B * i , B D * i using the transmission rate β and the migration of infested bees back to the healthy subpopulation using the grooming rate γ Healthy hive bees : Infested hive bees : Healthy drones : Infested drones : where is the ratio of the number of infested hive bees to the total number of hive bees. The term R is consistent with the SIR (Susceptible, Infected, Recovered) model of infectious disease [27]. These equations allow for a hive or drone bee to be infested with mites either through infested pupae maturation or as adults through the transmission term. Because we assume that infested hive bees are less effective in tending to larvae, we modify the ideal hive to larvae ratio (R H L ) ideal in (3) using the weighted mean in the presence of infested hive bees, where (R H L ) * ideal = 3.0.

Proportion of Drone Eggs
The proportion Z of eggs that are selected to be drone eggs is determined by the equation that was provided by Martin [18] where D L is the day length which is determined by the day of the year and the latitude. We find that the equation produces fairly similar rates when compared to the rates of drone egg production used by Sumpter and Martin [20] in the spring 1%, summer 3.3%, and autumn 1%. The day length is computed using the equations that were provided by the CBM model in Forsythe et al. [28] and Brock [29].

Mite Reproduction
While the female foundress mite lays up to 5-6 eggs within the capped cell, surviving daughter mite offspring are produced at the rate of 1.3-1.45 for worker brood and 2.2-2.6 for drone brood [5]. The reproductive rate is higher in drone brood because drone brood remain capped for a longer period [5]. Martin [26] uses 1.01 and 2.91 for the reproductive rate for mites emerging from hive worker brood and drone brood, respectively, to account for a multitude of factors including mite infertility. DeGrandi-Hoffman et al. [30] use 1.5 and 2.6 for the number of surviving daughter mite offspring in worker and drone brood and we adopt these values in our model. Sumpter and Martin [20] assume the percentage of time mites spend in capped pupae cells is 75%. While Kang et al. [21] state that the phoretic period may last 4.5-11 days when larvae are present, Martin notes that the mean number of days a mite spends as a phoretic mite is 4-6 [26], which slightly underestimates the 75% percentage if the mite lifespan is 27 days. Our model distributes 75% of total mites to available larvae cells on the day before they are capped. Let us refer to the mites that need to be distributed to available larvae as invading mites. Invading mites remain in capped cells for 12 days in worker pupae cells and 14 days in drone pupae cells. We assume that 85% of invading mites invade drone brood first up to four mites per cell due to the high propensity (5.5 to 12.1 times) of mites to invade drone cells compared to worker cells [14]. If there are insufficient drone brood to hold the invading drone mites, the remaining mites that did not invade a drone cell and 15% of the initial pool of invading mites are then distributed to available worker brood up to four mites per cell. The reduction in mite offspring due to multiple infestations of foundress mites is determined by using the values provided by Martin [26]. Specifically the number of offspring is reduced by 0.91, 0.86, and 0.60 in cells with two, three, and four foundress mites per cell in drone cells and 0.84, 0.65, and 0.66 in cells with two, three, and four foundress mites per cells in worker cells. The survival rates of the worker and drone pupae with multiple mite infestations are also reduced according to DeGrandi-Hoffman and Curry [17]. The survival rates are reduced 10%, 20%, and 40% for pupae with two, three, or four mites per cell respectively. A female mite usually has 2-3 reproductive cycles during her lifetime [5]. Martin [26] assumes the number of reproductive cycles for a mite to be 2.4. DeGrandi-Hoffman and Curry [17] only allow 60% of mites to invade a cell a second time. We only consider mites that are 16 days old or younger to be eligible for brood invasion, which allows for a maximum of two invasions.

Code and Computational Times
The code is written in MATLAB but also runs with Octave. Octave can be downloaded for free at gnu.org/software/octave. The free code can be downloaded at the github link: https://github.com/ davytorres/beecode-with-mites/. Computational times vary on different computers and versions of MATLAB. A two-year simulation with a time step of 0.1 days runs in approximately 0.07 s on a Dell Ultrabook with a 2.6 GHz Intel Core i5 processor with 8 GB of RAM with MATLAB R2018b, while the same simulation runs in approximately 0.12 s on a MacBook Air with a 1.6 GHz Intel Core i5 processor with 8 GB of memory with MATLAB R2017a. Computational times can increase by orders of magnitude with earlier versions of MATLAB or Octave. The same simulation runs in approximately 4 s on a Dell Inspiron 3671 with a 2.9 GHz Intel Core i5 processor with 8 GB of installed memory with MATLAB R2014a and in 83.7 s with Octave 5.2.0 on the same machine. Despite the differences in computational times, the algorithm is fairly efficient and sensitivity studies can be conducted with two interacting variables. Figures 1 and 2 show a two-year simulation of a honey bee colony with mites with separate plots for eggs E t (i), larvae B t (i), total pupae P t (i), infested pupae P * t (i), total hive bees H t (i), infested hive bees H * t (i), total foraging bees F t (i), and infested foraging bees F * t (i) (see Table 1). Note that the total caste subpopulation (e.g., H t (i)) includes the infected caste (e.g., H * t (i)) subpopulation. The survival rates in Table 3 and the simulation parameters in Tables 4 and 5 are used. In the simulation, the mites do not affect the survival rates of any of the bee castes, and the hive is provided with an infinite supply of food. Differences in the number of bees in each caste and the height of each caste graph are mostly determined by the number of days bees spend in each caste, but are also determined to a lesser degree on the survival rate of each caste and pheromones in the hive caste. The number of drones D t (i) is also shown, but multiplied by 10 for plotting purposes, since the number of drones is relatively small. The location of the peaks are staggered beginning with the egg caste, and continuing with the larvae, pupae, hive, and foraging castes and reflect the maturation time spent in each caste. We note that the mite population peaks later than the peak of the foraging bee population. DeGrandi-Hoffman et al. [30] also show that mite populations are the largest in September through November. Figure 3 shows the percentage of eggs that become drones using (15) at latitude 35 • . The graph has an average percent value of 2.5% over the two years during the summer. Figure 4 shows the positive values of the percent daily increase in the mite population,

Results
corresponding to the mite population curve in Figure 2. After a period of oscillation, the rate relaxes. The average rate in the first year is 2.6% and the average rate in the second year is 1.3%. The reported percent rate of growth of the mite population in virus-free colonies ranges from 2.1% [31,32] to 2.2-2.5% [33]. Figure 5 shows the percentage of mites that spend their time in capped cells during the summer season. We see that approximately 70% of the mites reside within capped cells. Martin [26] predicts that 50% to 70% of mites reside within sealed brood with his model.     (16). After a period of oscillation, the rate relaxes. The average rate in the first year is 2.6% and the average rate in the second year is 1.3%.   In Figures 6 and 7, a simulation is performed with mites under the same conditions as Figures 1 and 2, except the mites increase the mortality rates of the infested pupae, hive, drone, and forager castes shown in Table 3 by four times. We also apply (14), where (R H L ) * ideal = 3.0 in order to require additional infested hive bees to effectively attend to larvae. This requirement is not enforced in Figures 1 and 2 The differences in the bee castes are plotted. For example, in regard to hive bees, the difference is plotted in Figure 7. Figure 6 shows that the simulation with unchanged mortality rates (Figures 1  and 2) has larger numbers of larvae, pupae, and drones when compared to the simulation with increased bee mortality. However, the differences are less than 200 in each category. Figure 7 shows that more mites exist in the simulation with unchanged bee mortality, because the mites are negatively affected by increased bee mortality. In addition, the simulation with unchanged bee mortality has significantly larger numbers of hive and foraging bees. The number of hive bees that persist late into the summer season is significantly reduced with increased bee mortality and the reduced number of hive bees does not allow the colony to survive past the first year. Figure 8 shows how the maximum number of hive bees varies in the second year when 32 different values of the grooming rate ranging from 0 to 0.15 and 32 different values of transmission rate ranging from 0 to 0.2 are used for a total of 1024 simulations. In this sensitivity study, mites increase the mortality rates of the infested pupae, hive, drone, and forager castes shown in Table 3 by four times. We also apply (14), where (R H L ) * ideal = 3.0. With the exceptions noted above, the simulation parameters in Tables 3-5 are used. A sharp line separates the conditions under which the colony can be sustained (yellow) from the conditions under which the colony perishes (blue) beyond the first year. Figure 8 shows the importance of the grooming rate in sustaining the colony. Figure 9 shows how the maximum number of hive bees varies in the second year when 32 different grooming rates ranging from 0 to 0.1 and 32 different reproductive rates ranging from 1.5 to 3.0 for the drone foundress mite are used. As with Figure 8, in this sensitivity study, mites increase the mortality rates of the infested pupae, hive, drone, and forager castes by four times. A transmission rate of β = 0.2 is used. We also apply (14) where (R H L ) * ideal = 3.0. With the exceptions noted above, the simulation parameters in Tables 3-5 are used. Again, a sharp line separates the conditions under which the colony can be sustained. Figure 9 shows the importance of the higher reproductive level of mites in drone cells. Mite populations are significantly reduced in simulations of colonies without drones, despite the relatively small number of drones when compared to hive bees. We note that beekeepers do remove drone brood in practice to control mites [34].  Tables 4 and 5. Mites increase the mortality rates in Table 3 of infested pupae, hive, drone, and foragers by four times. Only the first year is shown since the colony does not survive past the first year. In this figure, the difference in the bee subpopulations (larvae, total pupae, infested pupae, total drones, infested drones) with unaffected mortality (Figures 1 and 2) and the bee subpopulations with increased mortality is shown.  Tables 4 and 5. Mites increase the mortality rates in Table 3 of infested pupae, hive, drone, and foragers by four times. Only the first year is shown since the colony does not survive past the first year. In this figure, the difference between the number of bee subpopulations (total hive bees, infested hive bees, total foraging bees, and infested foraging bees) with unchanged mortality (Figures 1 and 2) and the bee subpopulations with increased mortality is plotted. Figure 8. Effect of grooming rate γ and transmission rate β (10)-(13) on the second year maximum hive bee population. The transmission rate models the rate at which uninfested bees become infested. The grooming rate models the rate at which infested bees becomes uninfested. A sharp line divides the conditions under which the colony can be sustained (yellow) versus when the colony perishes (blue). Figure 9. Effect of grooming rate γ (10)- (13) and drone reproductive rate on the second year maximum hive bee population. The grooming rate models the rate at which infested bees becomes uninfested. The transmission rate β = 0.2. A sharp line divides the conditions under which the colony can be sustained (yellow) versus when the colony perishes (blue). The figure shows that drones play an important role in sustaining the mite population due to the higher reproductive rate of mites in drone capped cells.

Discussion
We have developed a mathematical model of a honey bee colony with mites, which tracks healthy and infested pupae, hive bees, drones, and foragers. The model is based on differential equations that solve for the number of bees in each day in the life of the bee. Model simulations shed light on the effect of different parameters on colony sustainability. Specifically, the grooming rate, transmission rate, and reproductive rate of the foundress mites are all important parameters that determine whether a colony will survive or perish.
Our model shows that the mite population is significantly reduced in the absence of drones despite the relatively low number of drones compared to hive bees. The importance of drone brood to mite population growth has also been reported in other models [13,17,35]. Our model also shows that grooming is an effective mechanism for controlling the growth of mites. The absence of grooming in the model can lead to very high populations of mites and the importance of grooming to bee survival has been noted by other authors [13,36,37].

Conclusions
The observations from our mathematical model could lead to continued or new experimental investigations in order to test the impact drone reduction [23] and increased grooming rate have on colony sustainability and ultimately help guide bee management processes (e.g., removal of drone brood and selective breeding) to improve colony health.
Future research efforts would be aimed at adding the influence of viruses and the increased mortality that they bring to specific castes. A computational comparison of available models of mite growth would also be useful to quantitatively compare the differences and similarities of facets within each model.