Predicting Honeybee Colony Failure: Using the BEEHAVE Model to Simulate Colony Responses to Pesticides.

To simulate effects of pesticides on different honeybee (Apis mellifera L.) life stages, we used the BEEHAVE model to explore how increased mortalities of larvae, in-hive workers, and foragers, as well as reduced egg-laying rate, could impact colony dynamics over multiple years. Stresses were applied for 30 days, both as multiples of the modeled control mortality and as set percentage daily mortalities to assess the sensitivity of the modeled colony both to small fluctuations in mortality and periods of low to very high daily mortality. These stresses simulate stylized exposure of the different life stages to nectar and pollen contaminated with pesticide for 30 days. Increasing adult bee mortality had a much greater impact on colony survival than mortality of bee larvae or reduction in egg laying rate. Importantly, the seasonal timing of the imposed mortality affected the magnitude of the impact at colony level. In line with the LD50, we propose a new index of "lethal imposed stress": the LIS50 which indicates the level of stress on individuals that results in 50% colony mortality. This (or any LISx) is a comparative index for exploring the effects of different stressors at colony level in model simulations. While colony failure is not an acceptable protection goal, this index could be used to inform the setting of future regulatory protection goals.

Introduction

A number of stressors have been implicated in honeybee losses in many parts of the world[1], including habitat loss;[2] viral diseases;[3] parasites such as Varroa destructor,[4,5] which can be a disease vector;[6−8] and use of pesticides.[9] Because these stressors may interact, it is difficult to predict how they change the colony dynamics separately and in combination. Moreover, because of the many feedback mechanisms in honeybee colonies, understanding the relationship between the effects on individuals and the colony level effects is not straightforward. Ecological modeling enables us to disentangle these interactions and explore them both separately and in combination, in fully controlled simulations. An innate difficulty in studying the effect of pesticides on honeybee colonies is the level of replication needed to capture low-level effects at the field scale.[10] The European Food Safety Authority (EFSA) has described specific protection goals for honeybee colonies, stating, “The magnitude of effects on colonies should not exceed 7% reduction in colony size”.[11] To assess whether this level of impact is occurring a minimum of 60 pairs (control and treatment) of colonies and fields are needed for each study.[11] If multiple stressors are to be studied even higher numbers would be needed. Ecological models can help in designing and targeting empirical studies, generating specific hypotheses that later may be tested experimentally, and can be used to assess the risk of environmental chemicals to honeybees.[11]

There have been many laboratory, semifield, and field studies showing both acute and chronic effects of pesticides on adult honeybees[12−19] and bee larvae.[20] For example, pesticides have been shown to affect foraging via acute mortality,[21] or alternatively from sublethal effects.[22,18] Other effects, such as reduced learning acquisition,[23] decreased rate of learning from olfactory cues,[24] and reduced communication for recruitment to foraging[25] have also been shown to occur, but the realism of these exposures is unclear.[12] In this study, we will concentrate mainly on hypothetical direct lethal effects on different life stages and reduced egg laying rate that could result from sublethal effects on the queen.[26] This does not capture the complexity of real exposure events, but it is important to compare the sensitivity of the colony to mortality of different cohorts at different times of the year in a controlled way.

As pesticides can affect individuals in a number of ways, determining the colony level impact of an individual effect is difficult. Feedback loops may compensate for moderate stresses (e.g., earlier onset of foraging if food stores are low[27]) or exacerbate other processes; for instance, less comprehensive care of brood as a result of high in-hive worker mortalities. The BEEHAVE model[28] is a suitable tool to investigate this complexity because it integrates in-hive processes and foraging activities to simulate interactions between colony and environment. The model consists of four modules:[28] (1) a landscape module, allowing the user to define a landscape of nectar and pollen in food patches; (2) a colony module, an age-based cohort model including processes such as nursing and care of brood; (3) a foraging module, an individual-based model[29] calculating the foraging activities on a particular day and the quantity of both nectar and pollen brought back into the hive; and (4) a Varroa and virus module simulating the population dynamics of the varroa mite and the transmission of viruses. The large number of procedures and feedback loops allow a comprehensive view of the impacts of stressors on the honeybee colony.[30] Here, we report simulations using the BEEHAVE model[28] to explore the colony-level impact of altering the mortality of a number of honeybee life-stages and reducing the egg-laying rate of the queen at different times of the year.

Because such simulations enable the user to examine a whole variety of stressors on individual bees, and the effects on the colony, finding a standard way of comparing the responses of the colony would be useful in risk assessment. For environmental chemicals the LD50 is the standard index used to describe the median lethal dose of a toxin i.e. that resulting in 50% subject mortality.[31] Here, we present an index to compare the impact of different imposed stresses on colony survival, the LIS50, describing the “lethal imposed stress” level resulting in a 50% colony mortality as predicted using the BEEHAVE model. We also present the LIS10, which predicts 10% colony failure from an imposed stress. We argue that these indices will be useful for comparing the impact of imposed stressors at colony level and could also inform the setting of pesticide protection goals in the future, once the indices have been applied to a wider variety of stressors and their variability has been quantified.

Methods

Model Parametrization

The BEEHAVE model[28] (BEEHAVE-Model Version 2014-03-04, free to download at www.beehave-model.net) was modified to increase the daily mortality of different life stages of bees in the colony from a defined day for a defined period to simulate potential effects of an exposure event (where exposure is defined as the period when toxic effects are imposed). We used the default setting as described in Becher et al.,[28] altering the landscape as explained below. Modifications to the model for the treatments below are outlined in the Supporting Information (Appendix 2).

The simulations started from first January with 10 000 worker bees in a colony and ran for 3 years, and each year had an identical annual weather cycle (based on maximal temperature and hours of sunlight at Rothamsted Research, Hertfordshire, U.K. in 2009). The colony was free from varroa and disease, as the purpose of the simulations was to look at effects of singular events increasing mortality in isolation. On the last day of each year, if fewer than 4000 adult bees were present, the colony is assumed to die due to winter mortality.[28] At the end of each three-year simulation the number of bees alive in the colony, or alternatively whether the colony had failed, was recorded.

Landscape

BEEHAVE allows users to define a dynamic landscape, giving values for the distance of each food patch to the colony, and the nectar quantity (L), nectar quality (sucrose concentration (mol/L)) and pollen quantity (kg) for each food patch on each day of the year. The simulations were set up in a very simplified and stylized modeled landscape: a single food patch was 1 km away from the hive offering 20 L of nectar and 1 kg of pollen each day of the year (although not representing the complexity of real landscapes, this enables tests of potential “exposure” in each month of the year in a controlled manner).

Imposed Stress

We ran simulations to contrast the effects of five different imposed stresses: reduced egg laying rate (ELR) of the queen, increased daily larval mortality, increased daily in-hive worker mortality, and increased forager mortality, applied daily or applied on each foraging trip. In reality, an exposure event may affect a combination of life stages over varying timeframes via different routes (nectar, pollen, honey, wax), but to specifically examine the sensitivity of different life stages, we chose a simplified set of simulations: examining increased mortality of individual life stages, during single exposure periods, when that exposure is assumed to be direct via consumption of nectar and pollen. We also ran the simulations for one combination of daily life stage mortalities (larvae, in-hive workers, foragers). Pupal mortality was not tested, as the pupae are in capped cells and do not receive food and therefore are unlikely to be exposed directly via nectar and pollen. In each simulation, a single stressor was applied for a continuous 30 day period each year. Duration of bloom of different crops differ widely as do persistence of different pesticides, but we here chose a 30 day exposure period as typical. Timing-dependent effects were investigated by running scenarios with the 30 day exposure period beginning on the first day of each month of the year. Imposed mortalities were applied as both multiples of the control value in the model and as set percent daily mortalities, while reduced egg-laying rate was only applied as a percent reduction. Testing a multiple of the control reflects the typical procedure of pesticide risk assessments.[11] We also ran simulations with set percent daily mortalities to determine the actual percentage of increased mortality that the colony could withstand.

Egg-Laying Rate

The daily egg-laying rate varies seasonally depending on the day of the year (for distribution, see Becher et al.[28]), and eggs are lost at a rate of ∼3% per day by default.[28] For the simulations, the egg laying rate (number of eggs produced on a particular day) was reduced by 25, 50, 75, and 90% for a period of 30 days with zero reduction applied as a control.

Stage-Specific Daily Mortalities

The daily mortalities of the larvae, in-hive workers, and foragers were altered in two ways: (1) Control daily background mortality was multiplied by a factor of 1, 1.5, 2, or 3. (2) Daily background mortality was set to a set percentage during the treatment period: one of 1, 5, 10, 25, and 50% each day.

The control daily background mortality is typically low: ∼0.4% for the in-hive workers[28] and foragers and ∼1% for the larvae.[28] For the larvae, this mortality does not include the chance of dying from lack of food or brood care.

Forager Mortality During A Foraging Trip

Increased foraging trip mortality was simulated in two ways, similar to the stage-specific daily mortality simulations: (1) control value was multiplied by 1 (control), 1.5, 2 or 3; and (2) set values of 1, 2.5, 5, 7.5, and 10% mortality per trip. Lower values than for the daily mortalities were used because foragers take multiple trips on a single day, so the majority of foragers may die if a forager had a 25 or 50% chance of dying on each trip, reducing the impact of higher mortalities as each bee can only die once. These settings were used to simulate pesticide exposure at levels high enough to cause death in the foragers through either immediate acute mortality, gradual weakening during the return flight, or through behavioral changes leading to impaired orientation and consequent homing failure.[18] Foraging mortality in the model depends on the duration of a foraging trip and is applied before an individual forager returns to the colony. For the single food patch present in the simulations, the mortality is ∼1.5% for nectar foragers and ∼0.9% for pollen foragers under control conditions (values taken from BEEHAVE model during control simulations). Although these values can vary during a day (as handling time of a food patch is increased when the patch is depleted), enough nectar and pollen are provided at the patch that this variation is negligible.

Simultaneous Daily Mortality of Larvae, In-Hive Workers and Foragers

To simulate an event in which several life stages are affected, the mortality of each of the larvae, in-hive workers, and foragers were all modified simultaneously in two ways: (1) Control daily background mortalities of larvae, in-hive bees and foragers were multiplied by a factor of 1, 1.5, 2, or 3. (2) Daily background mortalities of larvae, in-hive bees, and foragers were set to a set percentage during the treatment period: one of 1, 5, 10, 25, and 50% each day.

While it is unrealistic that different life stages are affected with identical effect levels, this scenario demonstrated the colony’s sensitivity to multiple effect types.

Each scenario was run for 30 replicates, with the mean number of live bees at the end of three years as output. For each combination of mortality type and 30 day exposure period, a linear regression was carried out between either the factor increase of the control or the percentage imposed mortality per bee and the mean number of bees alive per colony at the end of the 3 year simulation. The slopes of these regressions were plotted (Figure ), showing how sensitive the colony is to increased mortality of each life stage.

Figure 1

Calculated sensitivity of colonies to each stage mortality imposed for 30 days, calculated as the slope of the linear regression for the simulation data shown in Figures S1–S3. For each combination of imposed stress and treatment month, a linear regression was performed with the colony population at the end of the third year against the magnitude of the imposed stress. The graphs show the reduction in colony size (A) per percent decrease in egg-laying rate (ELR); (B) per multiple of the control background daily mortality for larvae, in-hive bees, and foragers; (C) per multiple of control background per-trip mortality; (D) per percent daily mortality of larvae, in-hive bees, and foragers; and (E) per percent daily per-trip mortality. (†) In these months, all levels of the combined mortality except 1% mortality lead to all colonies being lost; therefore, it was not possible to fit a linear regression. (‡) In these months, all levels of both the forager and the combined mortalities except 1% mortality lead to all colonies being lost; therefore, it was not possible to fit a linear regression.

LIS50: An Index for Comparing Lethal Imposed Stresses at Colony Level

To compare between the effects of these imposed stresses at colony level, we calculated a new index analogous to the LD50: the LIS50 was calculated as the level of imposed stress at the individual level that led, statistically, to 50% of the colonies dying in the BEEHAVE simulations within 3 years (using a threshold for survival of at least 4000 bees alive on the last day of each year).[28] In these simulations, the level of imposed stress was the percentage stage-specific daily mortality or percentage chance of dying during a foraging trip. The LIS50 was calculated using the “dose.p” function in R’s MASS library on a generalized linear model (GLM) built using data on the number of colonies alive after increased mortalities were imposed. For each of the imposed stresses in question, the mortality was applied from 0 to 100% (in 5% increments) for each month of the year (separately) with 50 replicates. For foraging mortality per trip, preliminary runs showed that colony death occurred when foraging mortality was 40% for all tested months, so higher mortalities were not tested.

The LIS50 was chosen for its theoretical parallel to the LD50, but a LIS could be calculated for any percentage of colony failure (x) that is of interest, for example, LIS10 figures are also presented, predicting the level of stress resulting in 10% colony deaths.

Results

Reducing Egg-Laying Rate (ELR)

Reducing the ELR for 30 days had only a moderate impact with none of the colonies dying in any of the simulations (Figure S1). A reduction of the daily egg laying rate by 90% (i.e., to 10% of the control) in June led to the average colony size at the end of three years being reduced by 35% of the initial population (Figure S1). Between April and August, each percent reduction in ELR led to only 50 fewer bees per colony after three years (Figure A). Nevertheless, colony dynamics were affected to varying degrees depending on the season and the reduction in the egg laying rate (Figure S1).

Effect of Increasing Mortality As Multiple of the Control

The colony was not highly sensitive to an increase in daily mortality of the larvae or the in-hive workers within the tested range (Figure B; Figure S2A,B). For the larvae, the control background mortality was already low (∼1%) and the majority of larvae that died in the control simulations did so from a lack of resources (food or brood care). Therefore, small larval losses from increased background mortality could be compensated in the model by allowing resources to be spread among remaining larvae, reducing mortality from a lack of these resources. As with background larval mortality, the control value of daily in-hive worker mortality in the model is small (∼0.4%), such that trebling it equates to 1.2% daily mortality and does not result in large losses over the course of the month. The impact of increasing the control daily forager mortality was low when imposed in January to August, but the colony was sensitive to increased mortality imposed in September to December (Figure B; Figure S2C). The critical threshold for colony survival in the BEEHAVE model was applied on the last day of December; therefore the colony had the whole year to recover from increased individual mortality applied in January, before winter survival was calculated. When the increased foraging mortality was applied to foragers at the food patch on each successful foraging trip, the impact on the colony was much larger (colony reduced down to almost 4,000 bees in June) (Figure C; Figure S2 D), as the mortality was applied many times per day and background mortality is higher than for in-hive life stages. This impact on the colony was likely due to the decreased food stores in the colony. These effects were particularly strong if the stress was imposed during the summer months when foragers were making the most foraging trips. For simulations multiplying the stage specific control mortalities, per-trip foraging mortality was the only single imposed stress to lead to colony failure with 3x mortality in June leading to 77% colony survival. When the mortality of larvae, in-hive workers and foragers were applied simultaneously as a multiple of the control mortalities, the impact on the colony was similar to the worst case equivalent single mortality (Figure B, Figure S4), the single life-stage daily mortality to which the colony is most sensitive when applied as a multiple of default.

Effect of Increasing Mortality by Set Percentage

The largest impact on the colony from larval mortality came when the effect was applied in one of the months between April and August (Figure D,E; Figure S3A). During this period, the colony has a lot of larvae as it is building to peak numbers, and increased mortality reduced or delayed this peak (Figure S5). Very high larval mortality in May and June led to the colony population being reduced to between 2000 and 3000 individuals, and winter mortality was high (Figure S3A). When larval mortality was increased (illustrated in Figure S5 for a level of 25% daily mortality), the resulting loss of larval numbers had the effect of reducing deaths due to lack of food or care during the treatment period, as these became more readily available for the surviving larvae. This feedback allowed the colony to compensate for moderate increases in larval mortality. However, high larval mortalities during summer led to a reduction in the worker population, which in turn led to a further peak in larval mortality 1 week after the end of the treatment period (Figure S5).

The modeled colony was sensitive to losses of adult bees in most months (Figure D,E; Figures S3B–D, S6, and S7). With respect to daily mortalities, the colony was more sensitive to losses of the younger in-hive workers than to the older foragers. During the period of April–September, the same period in which brood mortalities had a noticeable impact (Figure D; Figure S3), a daily in-hive worker mortality of over 25% led to all colonies being lost between May and July (Table S2) and a 5% daily mortality led to up to 40 000 more in-hive worker deaths over the course of the month (Figure S7). Loss of in-hive bees led to a large increase in brood loss from lack of care or food over the rest of the year (Figure S6C) and reduced the honey stores in the hive (Figure S6E). This combined stress was very damaging to the colony. Outside spring and summer, high daily forager mortality was devastating to the colony (Figure D; Figure S3C). This is because during the autumn and winter only few eggs are laid, so, the colony consists primarily of older bees still termed “foragers” (even though they rarely exited the colony).

The results of these simulations also highlighted the potential sensitivity of the colony to patch-specific forager mortalities, experienced on each foraging trip (Figure E). There was little effect at the very beginning or end of the year due to the lack of foraging activity at these times. Between May and October, there was a large impact on the colony from increasing this foraging mortality (June was the most sensitive month, as there was more time to forage in June than other months). A 5% mortality at the food patch applied in June led to an average colony size of ∼1000 bees at the end of the 3rd year (Figure S3D) with only 5 of 30 replicate colonies surviving (17% Table S2).

When the mortalities were applied as a set percent to several life stages simultaneously, this had a consistently higher impact than the worst case individual daily life stage mortality in each particular month (forager mortality in winter and in-hive mortality in summer) (Figure D). When colonies were subjected to combined daily mortalities of over 10% at any time of the year, then no colonies survived (Table S2).

LIS50: An Index for Comparing the Lethality of Different Stressors on Colonies

The LIS50 values represent the statistical likelihood that a certain imposed stress will lead to 50% colony failure for the specific control conditions used in the model (in this case calculated after 3 years). With agreement on appropriate control settings, the LIS50 could be standardized for use over any number of months or years, depending on the sensitivity required.

Table contains LIS50 values for four stressors imposed for 30 day exposures in 4 different months. The months were chosen to be those when foragers are active and crops flower and the colony is therefore most likely to be exposed. A low value indicates that a low daily percentage mortality imposed on individuals led to high colony failure and, therefore, identifies stressors to which the colonies are most sensitive. No values for the reduction in egg-laying rate were given, as no 30 day reduction of egg-laying rate led to colonies dying for any of the chosen months.

Table 1

Daily Percentage Mortality of Specific Honeybee Life Stages Required over 30 Days to Statistically Kill 50 and 10% (LIS50 and LIS10) of Colonies over 3 Years in an Otherwise Beneficial Environment (Ample Food and No Pathogens) ± Standard Errora

	larval mortality per day (%)		in-hive mortality per day (%)		forager loss per day (%)		forager loss per trip (%)
	LIS50	LIS10	LIS50	LIS10	LIS50	LIS10	LIS50	LIS10
April	> 100	76 ± 3	> 100	76 ± 5	> 100	68 ± 3	23 ± 0.4	18 ± 0.5
May	47 ± 1	33 ± 1	11 ± 0.2	9 ± 0.3	68 ± 1	49 ± 2	16 ± 0.4	8 ± 0.7
June	63 ± 1	31 ± 2	7 ± 0.1	6 ± 0.2	20 ± 0	16 ± 0.7	4 ± 0.1	3.6 ± 0.1
July	81 ± 1	38 ± 5	11 ± 0.2	9 ± 0.2	25 ± 1	18 ± 0.9	6 ± 0.2	5 ± 0.3

Values of >100% imply that in all of the simulations, 50% colony loss was not reached. No simulations imposing reduced egg-laying rate lead to colony loss.

The colonies were most sensitive to in-hive worker daily mortality and to per-trip foraging mortality (Table , Figure ). The in-hive worker mortality started to become very influential in the period between April and May, when the respective LIS50 went from over 100% in April to 11%, in May and 7% in June (Table , Figure ) as the brood nest was growing exponentially, requiring a large workforce of nursing bees, and the colony structure moved from mostly foragers (the overwinter bees) to more younger in-hive workers. The LIS50 for the forager mortality per trip has a relatively low value during the summer, with a value of just 4% in June. Such mortality was applied many times a day to foragers, especially those that were particularly active within the treated patch so the cumulative daily mortality was higher. The colonies were also sensitive, but to a lesser extent, to daily forager and larval mortality. Daily forager mortality had the largest impact in June and July, larval mortality in May and June.

Figure 2

Imposed stress response curves for the four individual mortalities investigated: (A) daily larval mortality, (B) daily in-hive worker mortality, (C) daily forager mortality, and (D) forager mortality per foraging trip. These show the % survival of 30 colonies for each of the varying mortalities, at the end of three years. Different colored lines are shown for these mortalities applied for 30 days in (yellow) April, (orange) May, (red) June, or (black) July. The intercept between a response curve and the solid horizontal line indicates its LIS50 (Table ), while the intercept with the dashed horizontal indicates its LIS10.

Importantly, the different imposed stresses had their greatest impact (smallest LIS50) at different times of the year. Specifically, larval mortality had the lowest LIS50 in May, whereas the two forager mortalities and the in-hive mortalities were at their lowest in June. The reason was that increased larval mortality led to a reduction of the adult in-hive population later in the year (Figure S5A) therefore reducing larval population in May could lead to reduced in-hive worker population in June (and all impacts which arise from that).

Alternative thresholds of colony failure may be explored with LIS: so LIS10 values have also been included in Table showing the level of stress causing 10% colony failure. Although the LIS10 is likely to be quite variable, it gives useful information when used in conjunction with LIS50: for certain life stages in certain months (in-hive bees and foraging trip mortality), LIS10 and LIS50 were remarkably similar (Figure B,D) where the stress response curve was so steep that there were effectively tipping points when any increase in daily mortality rates (imposed for this period of 30 days) led to all colonies failing.

Discussion

Using a set of simplified scenarios, our results showed a large variation in the impacts of imposed stress on the honeybee colony depending on both the demographic stage targeted by the imposed stress and the time of year in which the stress is applied. Imposed stress on all stages, except for daily forager mortality, led to highest colony losses from April to August. Imposed stress to the adult workers was, most often, more damaging to the colony than effects on the brood; and imposed stress applied to in-hive workers had a larger impact in late summer than the rest of the year. This was partly a consequence of the large proportion of in-hive workers, and that individuals were in-hive bees for a relatively long time such that cumulative stage mortality was higher, even if daily mortality was the same as for other life stages. In-hive bees (Table S1) would have consumed a lot of resources in their development but as yet not started generating them for the colony by foraging for nectar and pollen.

When the stress was imposed as a multiple of the default, there was little impact on the colony, except in the case of mortality calculated per foraging trip. This suggests that similar fluctuations in daily mortality are insufficient to cause colony loss in isolation, whereas foraging mortality calculated on a per trip basis carries a higher risk to the colony because a forager can perform many trips per day. When mortality of a life stage is a set percentage, > 10% daily (or per-trip) mortality for 30 days is very damaging in all cases.

The impact of decreased egg laying rate was sizable but not lethal to the colony and other models show similar results. Bromenshenk et al.[32] (using the PC BEEPOP model) find that eggs are the least damaging of the life stages to lose and adults are the most. Similarly, Schmickl and Crailsheim[33] (using the HoPoMo model) find that reducing the ELR to 60% reduced the number of bees in the colony on day 360 down by 40% from control. In real colonies, that may experience swarming events, there may be periods of around 3 weeks with no egg laying. Even when egg production is high, as during the summer, these require no food and only a small amount of nursing, so the investment lost with an egg is minimal. Large increases of larval mortality at particular times of year (Figure ) can have a significant impact on the colony. Individuals only spend 6 days as larvae so, as the stress was applied daily, an individual’s chance of being affected during their larval period is lower for the same daily mortality than for other life stages with longer developmental periods (Table S1). However, the colony invests honey and pollen in feeding larvae and has invested care from in-hive workers, and so, losing larvae will represent a net loss of effort to the colony. Figure S5 shows the feedback effect of high larval mortalities, at a sensitive time of year, leading to far fewer bees in the colony as a result of reduced potential worker population.

For newly emerged workers, a large amount of resources have been used to raise them to adulthood and they have not yet contributed to the colony. Depending on the state of the colony and the time of year, in-hive bees spend around 20 days before becoming a forager, meaning there is a long period over which the imposed stress can have an effect. Daily mortalities will build up quickly; a 5% daily in-hive worker mortality in August can lead to overall stage mortality of 65% and approximately 40 000 more in-hive workers dying (Figure S7) and further impacts from a 5% daily in-hive worker mortality (e.g., increase in larval loss and reduced food stores) are shown in Figure S6.

Figure shows that the colony was more sensitive to in-hive worker losses than forager losses during summer, and the colony was more sensitive to forager mortality toward the end of the year. At the beginning and end of the year, the modeled colony contains mostly foragers; no new workers were emerging, there was little or no foraging taking place, and in BEEHAVE, existing workers are classified as foragers once reaching a certain age. Therefore, at these times of year, impacts from forager mortality should be seen as impacts from general adult mortality. Hence, high daily adult mortalities can heavily reduce either the colony’s ability to survive over winter or the colony’s ability to build resources early the following year. High foraging mortality will also trigger in-hive workers to become foragers at an earlier age, reducing the age of first foraging. In reality, such precocious foragers may not be as successful as older foragers, resulting in further stress to the colony.[36]

The colony is very sensitive to high percentages of this combined mortality throughout the year (Figure D) but, in many months, not much more sensitive than the worst-case daily life-stage mortality at any one time period. An explanation could be that the loss of certain life stages can lead to the loss of other life stages, and hence, removing, for example, both in-hive bees and larvae will not necessarily cause more damage than only removing in-hive bees, as the larvae would have died anyway due to a lack of brood care. At many time points in the year, there is one life stage in the model that the colony is highly sensitive to losing, but at other times, the colony is more sensitive to losing multiple life stages (e.g., April and September), and these dynamics need further investigation.

Setting the BEEHAVE Simulations in the Context of Empirical Evidence

The set of simulations described here use a precisely defined exposure period, effects on single life stages, and a stylized landscape. In reality, the heterogeneity of the cropped landscape over time and space and the relative toxicity and persistence of different pesticides in the landscape, and in the hive, may lead to a diverse range of sublethal and lethal impacts on individual bees at different life stages. The next steps in using BEEHAVE to examine more realistic scenarios will involve using detailed empirical evidence to capture those exposure routes and timeframes, for specific chemicals in precise locations, and a specific module for this is in development. There are many empirical studies showing how stressors affect individuals or, in some cases, the colony in the short term. Long-term, multiyear studies are available,[19,34] yet uncommon,[17] so the impact of imposed stresses over multiple years is not fully understood.[17] Sandrock et al.[14] find that 1.5 month exposure to two neonicotinoid insecticides through pollen patties starting in May leads to a 28% reduction in worker population in the following April, along with effects on brood size and food stores. Dively et al.[19] also find effects on colony strength and overwintering success after 12 weeks (May–August) exposure to diet patties with high (20–100 ppb) levels of imidacloprid. In contrast, three studies of honeybee colony growth and survival in the field, when exposed via natural foraging on flowering crops treated with neonicotinoids, have shown no significant impact of the pesticide exposure on the colonies (for clothianidin[16,35] and for thiamethoxam[34]). Carreck and Ratnieks[12] suggest that the levels of pesticide encountered by foraging honeybees are lower in the field than used in many lab experiments. In large-scale field studies, in which the bees are placed near treated crops to forage[16,34,35] the bees may have lower and more variable pesticide exposure than in studies where bees are fed with an artificial feed, with pesticide added at “field realistic” levels[14,18,19] and this may explain why the former studies often find less damaging effects. In addition field studies offering the colony a known amount of pesticide (such as Sandrock et al.[14] and Dively et al.[19]) find that the impact upon the hive from the pesticide can appear sometime after exposure. The simulations presented here show how this can occur within the model: Figure S5B shows how the number of larvae in the colony is affected by a 25% larval mortality in June. It is clear that there is an additional delayed impact likely due to a reduction in workers providing brood care.

These contrasts also highlight the difficulties of scaling from individual level effects to those at the colony level. The BEEHAVE model contains a large number of feedback loops, allowing in-depth investigation into how multiple stressors can disturb the colony dynamics in terms of mechanism, and which particular stresses are more damaging to the homeostasis of the colony. For example, with high forager mortality, worker bees become foragers earlier to compensate; this in turn may reduce the nursing force, increasing larval mortality.[36] A small increase in larval mortality can reduce the mortality of the surviving larvae from other causes such as lack of food or brood care, and reduces further losses, i.e. show a compensatory effect (Figure S6).

How do BEEHAVE Simulations Compare to Those of Other Models?

The BEEHAVE model is a useful tool in the risk assessment of stressors to bees as many potential stressors can be assessed simultaneously, and testable hypotheses can be developed. Indeed EFSA have recently published a review to suggest, with further development, BEEHAVE could be the model of choice for regulatory pesticide risk assessment.[37] Several models[32,38,39] have been used to explore the impact of pesticides on bee colonies. However, the models of Khoury et al.[38] and Bryden et al.[39] focus on limited portions of colony dynamics, and lack key processes required to accurately predict how a bee colony reacts to numerous stressors. One major feature lacking in these models is seasonality. We have shown that the time at which a stress is applied greatly affects the colonies’ response. A stress imposed in April has little effect, while the same stress imposed in June will devastate the colony. PCBEEPOP[32] is a model that includes colony dynamics similar to BEEHAVE; although, BEEHAVE also includes a number of factors, such as the landscape and foraging dynamics (including the flow of energy in the form of honey stores) integrated with the colony module, which may be key to understanding how pesticides can impact the colony.[30] The BEEHAVE model is the only tool to date that also includes a dynamic landscape module and weather providing the potential for climate or location-specific simulations, as well as integration with the foraging and varroa and virus modules, to apply many stressors to the colony at any one time, as would be happening to real colonies in the field. Further development of BEEHAVE, with a “pesticide module”, to ensure correct implementation of exposure routes from flower, via forager, into the colony is underway.

Can Such Simulations Be Used in Future Risk Assessment?

To quantify the impact of a pesticide on a hive, the EFSA guidance[11] classifies the magnitude of an effect by the % change in colony size. EFSA considers a change between 3.5 and 7% negligible; and a change larger than 35% to be large. EFSA use the Khoury et al. 2011 model to estimate what forager loss would be permissible for “negligible” change and find that forager losses of 1.5 × control for 6 days; or 2 × control for 3 days or 3 × control for 2 days would be permissible. We have shown that effects on fecundity or brood mortality are not as impactful on the colony as adult loss, so worker loss is a conservative measure of the damage possible from a pesticide to the colony. The modeled colony has a certain capacity for compensation, which varies with the life stage affected and time of year and durations, but once the compensation threshold is exceeded the colony is likely to fail. Levels of “background” mortality in the absence of pesticide exposure, depending on weather, forage quality and other stressors present, are likely to influence compensation capacity. Due to this compensation capacity and how it may vary with the health of the colony, a percent reduction of bees in the colony could have highly variable results on the health of the colony. Modeling, such as with BEEHAVE, could, therefore, help supplement the risk assessment procedure by teasing apart such dynamics.

LIS50 and LIS10 provide a tool to compare the effects of a variety of imposed stresses on the colony using BEEHAVE, treating the colony as an individual “super” organism by using the percent chance of colony mortality as a measure of sensitivity to imposed stress. Colony failure as a result of an introduced chemical is not an acceptable end point, but these indices (calculated from simulations) could be used to provide theoretical comparisons of the effects of different stressors on the colony, which may be informative in discussions of future regulatory risk assessment procedures, and protection goals. This study was designed with impacts of pesticides in mind, but impacts on the colony driven by Varroa destructor, related diseases, Nosema sp. or lack of forage sources could also be compared with LIS. Also, comparing LIS10 and LIS50 provides hypothetical evidence of which stressors, at which levels, may lead to colony tipping points, with the caveat that the tipping point will depend on the control scenario (e.g., forage availability and weather will all affect the colony’s capacity for compensation). The BEEHAVE model, together with the use of LIS50, allows consistent investigation into the impact of multiple stressors on the honeybee colony, and could be key for future risk evaluation.