Age structure is critical to the population dynamics and survival of honeybee colonies.

Age structure is an important feature of the division of labour within honeybee colonies, but its effects on colony dynamics have rarely been explored. We present a model of a honeybee colony that incorporates this key feature, and use this model to explore the effects of both winter and disease on the fate of the colony. The model offers a novel explanation for the frequently observed phenomenon of 'spring dwindle', which emerges as a natural consequence of the age-structured dynamics. Furthermore, the results indicate that a model taking age structure into account markedly affects the predicted timing and severity of disease within a bee colony. The timing of the onset of disease with respect to the changing seasons may also have a substantial impact on the fate of a honeybee colony. Finally, simulations predict that an infection may persist in a honeybee colony over several years, with effects that compound over time. Thus, the ultimate collapse of the colony may be the result of events several years past.

Introduction

As honeybee populations continue to decline on a global scale [1], research efforts have been directed at identifying the underlying causes [2-5]. These efforts are necessitated by the ecological [6] and economical importance of honeybees worldwide [7-9]. To date, much of this research has been focused on the effects of pesticide and insecticide exposure on honeybee health and ultimately colony fitness [4,10]. Such hazards may cause injury or death to foraging bees, forcing surviving bees to begin foraging prematurely which will then disrupt the dynamics of the colony, ultimately leading to colony collapse [10-12]. Further research has focused on the effects of parasitism and disease [13-18], for example, the effects of Varroa destructor on honeybee colony dynamics [15-18], the effects of the microsporidian parasite Nosema ceranae [19] and the effects of communicable infections more generally [13].

Recent work has also explored the combined effects of two stressors such as disease, limited biodiversity and/or exposure to pesticides, which may create conditions detrimental to honeybee colony survival [20-25]. Simulation packages have been developed to model these effects using realistic parameter values [26,27]. For example, Pettis et al. [28] explored the interplay between environmental hazards and the subsequent susceptibility to Nosema ceranae. In previous work [13], we showed that the fate of a honeybee colony in the presence of disease is dependent on seasonality, particularly the onset of winter.

Previous mathematical models of honeybee dynamics have not considered the effects of age structure within the colony, yet the duties of female bees, which constitute the main work force of the hive, are determined primarily by their age [29]. Specifically, the structure of a colony hinges on both bee morphology and age [30-32]: caste polyethism differentiates a queen bee from a female worker bee [33], while age polyethism determines the functions of the worker bees [29].

Thus, the age distribution within a honeybee colony is an intrinsic property and this property is therefore important to colony survival, because it affects two key components in the dynamics of the colony, namely ongoing recruitment and death [29]. In this paper, we examine how the disease-free age distribution within a colony is altered in the face of a hazard, and how these changes affect previous predictions of colony survival.

We identify periods of increased vulnerability of the colony to the effects of winter, disease or a combination of the two, during which the colony would benefit from remedial actions (e.g. higher anti-microbial treatments, increased observation and management, etc.). In addition, we use the added dimension of age structure to explore the long-standing phenomenon of ‘spring dwindle’ whereby a drop in the number of bees within a colony occurs immediately after the end of winter [34]. While the phenomenon is generally suspected to be due to various stressors [35-37], the specific colony dynamics that lead to this phenomenon have not been established. The question of why the dwindle occurs after rather than during winter remains unanswered. The resolution of this puzzle is of particular value in the ongoing research efforts to understand and ultimately avert honeybee colony collapse, since spring dwindle leaves the colony in a particularly vulnerable state. Finally, we simulate colony dynamics in the presence of disease for multiple years to demonstrate the compounding effects of infection.

Model

We present a mathematical model that combines the disease-free demographics of a honeybee colony with the effects of seasonal changes and a disease that at first infects foragers, and then spreads to the rest of the colony. Earlier versions of this model were introduced by Khoury et al. in 2011 and 2013 [11,12] and were developed further in 2014 [13]. Analytical details of the model, including local and global stability of the disease-free equilibrium as well as a derivation of the basic reproduction number, R0, are presented in [14].

Briefly, the effects of the brood, guarding bees, as well as bees that work to repair the hive are neglected. The focus is solely on the remaining adult hive bees, H, which are responsible for ensuring the survival of the brood, and the foragers, F, which are responsible for bringing food, f, into the hive. The male honeybees, known as drones, are also neglected since they contribute only to reproduction [33]. In the presence of disease, the two classes of bees are further divided into the susceptible (disease-free) populations, HS and FS, and the infected populations, HI and FI.

Age structure

We incorporate an age structure into the age-independent equations [13], using the standard approach of McKendrick [38] by writing where HS(a,t) is the number of susceptible hive bees of age a at time t and u(a) is the age-dependent rate of recruitment to foraging. Juvenile hormone III regulates the age at which honeybees begin foraging [29], making older bees more likely to be recruited to foraging duties. As well, it has been observed that there is a minimum age, aR, before which bees cannot be recruited [39]. Behavioural maturation of hive bees is further regulated by a pheromone, ethyl oleate, that is produced by foragers [40] and has the effect of delaying the age at which bees are recruited to foraging duties. This so-called ‘social inhibition’ process reduces recruitment when the number of foragers in the colony is high [40]. We account for these biological processes by defining u(a) as where α is the maximum rate of recruitment and Hv(a−aR) is the Heaviside function such that recruitment cannot begin before age aR. Thereafter, the recruitment rate we are using increases sigmoidally with age. At age k, recruitment will be half the maximum rate of recruitment. 1/σ is the maximum proportion of bees that are foraging at any given time and N is the total number of bees given by The emergence of new hive bees is modelled as the left boundary condition to equation (2.1). This, along with other boundary conditions are described at the end of this section.

The second term on the right-hand side of equation (2.1) governs the disease dynamics within the hive. We approximate the transmission of disease as a mass action process and assume that, on average, hive bees and foragers transmit infection between or within classes at rate β. The total number of infected bees is given by

The hive provides substantial safety for bees that are confined to it [11,30]. We therefore assume that the natural death rate of healthy hive bees is negligible compared to the rate of recruitment.

The equation governing the dynamics of infected hive bees is given by In contrast with their healthy peers, infected hive bees are at risk of dying due to disease at an age-dependent rate d(a).

Susceptible foragers are recruited from susceptible hive bees, and suffer age-dependent natural death at rate μ(a). Their dynamics are therefore governed by

Infected foragers can either be recruited from infected hive bees, or from susceptible foragers that have become infected. If we assume this class is subject to a disease-related death rate of d(a), then their dynamics are governed by

Food, f, is brought into the hive by both susceptible and infected foragers. Although it is likely that infected foragers would be less efficient at this task [41], for simplicity, we take an average rate of food intake, c (g/day/forager). Food is consumed by foragers and hive bees at an average rate γ. Therefore, the amount of food available at time t changes according to

The above system of equations governing the dynamics of the colony is subject to the following boundary conditions:

The first condition represents the emergence of new adult bees, where L is the daily egg-laying rate of the queen and S is a survivability function, which determines how many eggs survive to adulthood. The brood needs both sufficient food and sufficient care from the hive bees in order to survive [42]. Moreover, it has been shown that there is a range of ages within which hive bees will take on this care for the brood, a minimum age, a, and a maximum age, a between 11 and 16 days old [32]. After this age, hive bees tend to transition to foraging duties, or possibly security or hive maintenance [33]. Therefore, we define the survivability function S as Here b is the amount of food required for half the eggs to survive to adulthood provided the brood has sufficient care, w is the number of care providers required for half the eggs to survive provided the brood has sufficient food, and H(a,t)=HS(a,t)+HI(a,t).

The dynamics described in this section relate to the main active season of the bee colony, which is defined as the time interval between the end of one winter and the beginning of the next.

Winter

The winter season is assumed to last 155 days or roughly five months, roughly corresponding to a humid continental climate [43]. Over winter, no new hive bees emerge (although some eggs are still being laid by the queen) and foragers return to the hive [33]. Accordingly, in boundary conditions (2.9) we set L=0, and in equations (2.1), (2.5)–(2.8) we set u(a)=0 and c=0. Owing to an extended lifespan of bees over winter [33], the death rate of hive bees is no longer negligible and is set to an average rate of corresponding to an average lifespan over winter of six months [33]. All bees are performing the same function over winter (keeping the hive warm [44]); therefore, we set the natural death rate of foragers equal to that of the hive bees, μ(a)=μw.

As winter ends, bees resume their normal age-related duties, and bees that were foraging before winter resume their roles as foragers. We assume that a 21-day transition takes place during which parameter values change linearly from end-of-winter to new active season values as follows: where that is, during winter hive bees of all ages may contribute to brood care. Note that new hive bees that emerge during the early spring come from eggs that were laid over the winter months, after completing phases of development to adulthood.

Equations (2.1), (2.5)–(2.8) form a system of integro-partial differential equations which we solved simultaneously to predict the population dynamics of a honeybee colony. The parameter values used in the solution are given in the electronic supplementary material, table S1.

In [45], the authors found that the age distribution in a colony is as follows: 41% of bees were one week old, 23% were two weeks old, 17% were three weeks old, 11% were four weeks old and 8% were five weeks old. Accordingly, the death rate μex(a) in our model was constructed to approximately match the experimental findings. To account for the whole range of ages and remain consistent with the five weeks of study in [45], we break the age range into five cohorts of 10 days each. The comparison between experimental data and our model is summarized in table 1. The death rate that best matched the observed data is defined by where parameter values were set such that the percentage of bees in their first week of life matches the experimental data. This death rate is relatively high for very young (approx. <10 days) or very old bees (approx. >40 days), as shown in figure 1. We see that the drop in population density after week 1, as observed in [45], can be explained by the high death rate of young foragers.

Table 1.

Comparison of experimental results with model results. Standard deviations for experimental results are given; z-scores computed from experimental means and standard deviations are given in parentheses next to model values.

cohort	experimental values [45]	model values, μex(a)	model values, μ(a)
1	41±6.4%	41% (0)	41% (0)
2	23±5.1%	25% (0.39)	30% (1.37)
3	17±3.7%	18% (0.27)	16% (0.27)
4	11±2.2%	11% (0)	7% (1.81)
5	8±1.3%	5% (2.31)	6% (1.54)

Figure 1.

Death rate distribution which provides the best fit to experimental results observed in [45]; parameter Cex=0.42.

We performed sensitivity analysis (see the electronic supplementary material, figure S1) and found that our results are insensitive to the details of the death rate distribution. We therefore use a simplified version of the natural death rate distribution to make the mathematics more tractable. This death rate is plotted in figure 2, and also compared with the experimental results in table 1, showing good agreement.

Figure 2.

The natural death rate, μ(a), and an example of disease-related death rate, d(a), motivated by deformed wing virus, which predominantly affects young bees [46]. During the active season, honeybees live an average of five to seven weeks [33], thus we define our death rates until the age of 60 days. Beyond this age few, if any, bees are still alive and we set their natural death rate as .

The disease-related death rate is motivated by deformed wing virus [46] that disproportionately affects young bees. The constants C, K1 and K2 are chosen such that both the natural and disease-related death rates are non-negative and have an average of 0.14 per day as in [13]. In other cases, we assume the disease affects all bees uniformly, i.e. d(a)=d. The effects of different death rate distributions are shown in the electronic supplementary material, figures S4–S12. The results are qualitatively the same.

Electronic supplementary material, figure S1, shows the age distribution within a disease-free colony under the natural death rate μ(a) shown in figure 2. This is used as a baseline distribution for all simulations that begin from equilibrium. We predict from this distribution that the mean forager age of a healthy colony is approximately 25 days, in agreement with [33,47].

Results

Figure 3 shows the effects of changing seasons on the number of bees in a colony in the absence of disease. The seasonal shift back to normal operating conditions is seen to bring with it a drop in the total number of bees as the colony emerges from winter. It is obvious from these results that a longer winter may have a deleterious effect on the colony since the effects of winter are largely negative. This is because the surviving bees are now much older and are subject to a higher natural death rate, while a new generation of younger bees has not yet emerged to replace them. This seasonal drop in the bee population, which has been referred to as ‘spring dwindle’ [34-37], is discussed further in the next section. The colony dynamics that produce the phenomenon have not been fully understood in the past because these dynamics emerge only when the age distribution within the colony is considered.

Figure 3.

Time course of total bee population in a disease-free colony. Shaded area denotes winter, during which no brood is produced. We see that the aged bees that restart the colony post-winter create an added risk to the colony, because of their relatively high death rate, until new brood is able to hatch and contribute to population rebound. The post-winter dip seen in the figure has been referred to in the literature as ‘spring dwindle’ but its dynamics have not been previously described.

Our next objective is to compare the dynamics of the age-structured honeybee colony model with the dynamics of the age-independent model such as that presented in [13]. To give the comparison a measure of equivalency, the same value of the basic reproduction number is used in both cases, namely R0=1.43, calculated using the formulae derived in [14]. Figure 4 shows that an age-independent model makes substantially different predictions about the timing and the severity of disease within a bee colony. This is further explored in the electronic supplementary material, figure S3, where we demonstrate that an age-dependent recruitment function both lowers the severity of a disease and delays the timing of an epidemic, while age-dependent death rates, μ(a) and d(a), trade off the two effects.

Figure 4.

The effect of age distribution on the dynamics of disease within the bee colony. The blue line represents the total number of susceptible bees and the green dashed line represents the total infected bees. The shaded area denotes winter. Disease begins with a single infected forager at t=112. We use the death rates given in figure 2. The two time courses show that the two models have significant differences in their predictions when simulating disease dynamics in the colony.

The vulnerability of a bee colony in early spring is explored further in figure 5 which illustrates that an infection which is endemic in a colony during one active season and one winter may pose the greatest risk to the colony during the small window in which the hive is recovering from the winter months. Moreover, the colony may suffer much greater losses after its second, third and fourth winters with disease (figure inset).

Figure 5.

Time course of the bee population over 4 years. (a) The total number of susceptible bees, HS+FS and (b) the total number of infected bees, HI+FI. An infection is introduced into a colony at equilibrium through a single-infected forager at time t=0. The inset shows the minimum population size of susceptible bees during the spring dwindle following each winter. There is a clear decrease in this minimum year after year, thus the risk of colony collapse increases year after year. The ultimate collapse of a bee colony may therefore be the consequence of events long past.

The time interval between the beginning of an infection and the onset of winter is an important determinant of the health of the colony at the end of winter, as illustrated in figure 6. The figure shows that the number of bees surviving to the end of winter, NW, depends on when the infection begins relative to winter. In the figure, the infection begins from a single-infected forager Δt days before winter. The time interval between the end of one winter and the beginning of the next is taken to be approximately 200 days (approx. seven months). The timing of the infection is seen to produce a point of highest vulnerability (local minimum in population size) which depends on the severity of the disease as represented by the basic reproduction number.

Figure 6.

Number of bees at the end of winter as a function of the time interval, Δt between the onset of disease and the beginning of winter, based on various values of the basic reproduction number, R0. The minima show times during which the colony is most vulnerable due to the combined effects of winter and disease. The colony is started from the equilibrium age distribution shown in the electronic supplementary material, figure S1. Infection is introduced at time t=200−Δt via one infected forager. The total number of bees is calculated at the end of winter. R0 values are calculated using the expression derived in [14].

The timing and depth of the points of highest vulnerability seen in figure 6 are tightly coupled to the time at which an infection peaks within the colony, which in turn is related to R0. Figure 7 shows this relationship, estimating the time of greatest risk to the colony, Δt*, with respect to the onset of winter, based on the basic reproduction number. Of course, for high values of R0 the colony is likely to suffer substantial losses regardless of when a disease occurs. In figure 6, we see, for example, that for R0=2.0, Δt* is unique, but when R0=2.2 there are two Δt* values producing the observed bifurcation in figure 7. To measure the relative significance of these minima, we use the metric where the maximum and minimum are computed for each value of R0, over the range of possible Δt. The larger the value of this metric, the greater the significance of the time of onset of disease, Δt.

Figure 7.

The timing and depth of points of highest vulnerability as a function of the basic reproduction number, R0. The colour bar reflects the depth of the minima seen in figure 6 (see text for details). Red indicates the most severe combinations of R0 and Δt*, where Δt* is the value of Δt in figure 6 at which a local minimum occurs.

Discussion

Our results demonstrate that the age structure that is inherent to all honeybee colonies has critical effects on colony survival that are not captured by previously explored age-independent models. Most notably, an age-dependent recruitment rate has a profound impact on the dynamics, severity and spread of an infection in a bee colony. In fact, our results predict that age polyethism in a honeybee colony has the beneficial effect of impeding the spread of infection within the colony. Thus, in practice, any disturbance in the natural age structure of a bee colony may also increase the severity of an infection within the hive. Results in the electronic supplementary material, figure S3, indicate that the natural age distribution within a colony creates a delay in the speed of spread of an infection. Again, this suggests that any intervention that would alter the natural age distribution may have the adverse effect of increasing the spread of infection within a hive.

Our explicit treatment of the colony age distribution also revealed the dynamics underlying the frequently observed phenomenon of ‘spring dwindle’ in which the total number of bees within a colony weakens rapidly immediately after the end of winter. While the phenomenon is generally suspected to be due to various stressors that occur during winter [35-37], our results suggest that declining spring populations are a natural consequence of honeybee colony dynamics, leading to seasonal vulnerability during which the colony may be particularly susceptible to other hazards. This post-winter dip is not captured by age-independent models because the phenomenon is due predominantly to the ageing of the bees over the winter months and the time required to replace them by younger bees in the early spring.

Simulations of colony dynamics over multiple winter seasons demonstrated that endemic diseases may weaken the colony year after year, increasing the risk of colony collapse during spring. Thus, the eventual collapse of a bee colony may be due to the compounded effects of an infection over several years.

The time interval between the onset of an infection and the beginning of an approaching winter is a critical determinant of the ultimate course and consequence of the disease within a bee colony. As illustrated in figures 6 and 7, this effect is most important for infections with a basic reproduction number close to 2. In this case, the colony is particularly vulnerable in the early spring, approximately five months before the onset of winter.

In conclusion, age structure is an inherent property of honeybee colonies, determining the division of labour and defining how age-related events such as disease and death unfold within the colony. Our results suggest that age structure should be a key consideration in future studies of honeybee population dynamics, and that accurate models of age structure will be necessary if we are to understand, and ultimately reverse, the increases in colony collapse observed in recent years.