BEEHAVE: a systems model of honeybee colony dynamics and foraging to explore multifactorial causes of colony failure

Summary A notable increase in failure of managed European honeybee Apis mellifera L. colonies has been reported in various regions in recent years. Although the underlying causes remain unclear, it is likely that a combination of stressors act together, particularly varroa mites and other pathogens, forage availability and potentially pesticides. It is experimentally challenging to address causality at the colony scale when multiple factors interact. In silico experiments offer a fast and cost‐effective way to begin to address these challenges and inform experiments. However, none of the published bee models combine colony dynamics with foraging patterns and varroa dynamics. We have developed a honeybee model, BEEHAVE, which integrates colony dynamics, population dynamics of the varroa mite, epidemiology of varroa‐transmitted viruses and allows foragers in an agent‐based foraging model to collect food from a representation of a spatially explicit landscape. We describe the model, which is freely available online (www.beehave-model.net). Extensive sensitivity analyses and tests illustrate the model's robustness and realism. Simulation experiments with various combinations of stressors demonstrate, in simplified landscape settings, the model's potential: predicting colony dynamics and potential losses with and without varroa mites under different foraging conditions and under pesticide application. We also show how mitigation measures can be tested. Synthesis and applications. BEEHAVE offers a valuable tool for researchers to design and focus field experiments, for regulators to explore the relative importance of stressors to devise management and policy advice and for beekeepers to understand and predict varroa dynamics and effects of management interventions. We expect that scientists and stakeholders will find a variety of applications for BEEHAVE, stimulating further model development and the possible inclusion of other stressors of potential importance to honeybee colony dynamics.


Setup & Running options
The interface tab of BEEHAVE provides several setup and running buttons that initiate and start the simulation run as well as offering the option to choose pre-defined scenarios.
The Setup button initiates the model settings called by the parameterization procedure whereas the Run button starts the model run over an indefinite period (or until the death of the honeybee colony). Alternatively you can run the model for defined time intervals from one day to many years by using the 1 Day, 1 Month, 1 Year and run X days buttons.
Furthermore, you can start pre-defined scenarios by the setup options of the BEEHAVE interface shown above. You may choose between the Default, 2 patches, 1 feeder, RRes, varroa and the beekeeping scenario.
When the Default scenario is initialized, all input variables will be reset to their default values, except for RAND_SEED, DotDensity, and INPUT_FILE. These default settings include 2 food patches located at distances of 1500 and 500 m from the hive, seasonal food flow as well as local weather data of Rothamsted from 2009, but do not include any beekeeping options or infestation by DWV-transmitting (deformed wing virus) Varroa mites.
The other pre-defined scenarios change some input variables. Thus, the beekeeping scenario implements a virtual beekeeper that feeds bees, harvests honey, merges weak colonies and performs Varroa treatments.
In addition, 1 feeder offers the opportunity to implement only one active feeder, constant food flow and a short, constant handling time for nectar and pollen collection.
The varroa scenario adds 20 Varroa mites, 50 percent of which are infected with the deformed wing virus (DWV), to the default settings.
The RRes scenario integrates the locations and sizes of real patches of forage into the model run from a landscape map of the area around Rothamsted, and includes seasonal food flow and varying flower handling times for nectar and pollen foraging. This gives a more realistic foraging landscape for the bees. These data are read in from an input file.
The dance slope of a dancing bee can be set to the mean (0.51) or maximum value (1.16) in order to calculate the number of circuits a bee dances for a patch depending on the patch quality (energetic efficiency) Default: 1.16 The dance slope and intercept of a dancing bee can be set to the mean (0.51 and -11.1) or maximum value (1.16 and -17.7) in order to calculate the number of circuits a bee dances for a patch depending on its energetic efficiency: e.g.
distance to the hive Default: dance slope (1.16) and intercept (0) add 1 kg pollen or 1 kg fondant to the colony's pollen store at any time during the model run runs the model for 1825 days to detect, if the official BEEHAVE code was changed The advanced button activityList offers the option to write all foraging activities of all active forager squadrons on that day to the 'Command Center'. So, the ID of the respective forager squadron, its flight distance travelled as well as its foraging activities per foraging round on that day are shown. Table 1.2 gives an overview of all possible foraging activities. frN/ frP a recruit finds a nectar / pollen patch and becomes a successful nectar / pollen forager eSn/ eSp a recruit found a patch that provides no nectar / pollen, and becomes a scout mSn/ mSp a recruit didn't find the advertised nectar/ pollen patch and searches a new patch Dn/ Dp a successful nectar / pollen forager dances and communicates EEF and ID of the patch rFnNF/ rFpPF a resting and unexperienced forager is recruited to an advertised flower patch by a dancing forager and becomes a nectar / pollen forager rFnxNF/ rFpxPF if the advertised nectar patch has higher EFF or the advertised pollen patch has a shorter trip duration than the known patch, the dance follower becomes a nectar / pollen forager and searches for the advertised patch Xnr/ Xpr an experienced forager recruited by following a dance, but start foraging at a known nectar / pollen patch with the same or higher EEF / trip duration than the advertised patch An/ Ap resting foragers that abandoned their nectar / pollen patch still rest AnSn/ ApSp an active forager that abandoned its nectar / pollen patch searches for a new one AfR an active forager that was not recruited in the foraging round before, abandons foraging and becomes a resting bee Rx a pollen or nectar forager that switches to resting N/ P collect nectar / pollen at a flower patch uS an unsuccessful searching bee finds no flower patch End end of the today's foraging period of the respective forager squadron the viewing of updates. Use the speed slider to control the speed of the model run, and use the chooser to adjust the frequency of viewing updates (from continuous to tick-based). Thus, continuous updates means that the Netlogo view will be updated many times a second, whereas the tick-based updates results in a Netlogo view only updated when the tick counter proceeds.
The ticks count daily as the model run proceeds. If you click the Settings button on the toolbar, then you will see the 'Model Settings' window (below) with a box to choose whether to 'Show tick counter' or not.

Input options
The BEEHAVE interface tab offers a large number of input options using 'Input fields', 'Switch' and 'Chooser' buttons (blue). These define the parameter settings, input data and a variety of adjustment options of colony and foraging processes, Varroa infestation and beekeeping practices for simulation experiments in order to investigate the effects of combined stressors on the performance of honey bee colonies. Thus, the following tables in this section give a review of the various input options in order to clarify these complex model Here, important input data and settings are specified. You can define specific feedback loops in dependence of actual colony composition and stores, food flow conditions and foraging properties for simulation experiments using "on" and "off" switches (table 1.4).
The Default scenario sets the constant handling time for nectar and pollen collection (ConstantHandlingTime) "off" and sets a season-dependent pattern in food availability (SeasonalFoodFlow) "on". Hence, the maximum amount of either nectar or pollen available at the "red" and "green" patches depends on the season (if ReadInFile = false). If ConstantHandlingTime is off, then higher depletion of pollen and nectar at the flower patch during the day results in longer handling times to collect nectar or pollen loads.
Moreover, the effect of queen's age on queen's egg-laying rate (QueenAgeing) is turned off.
This implies that the seasonal egg-laying is not reduced with increasing age of the queen. But the egg-laying rate of the queen is affected by the number of available in-hive bees that feed the queen as well as prepare and nurse the brood (EggLaying_IH).
All beekeeping options such as feeding bees (FeedBees) and harvesting honey (HoneyHarvesting) are switched off. and plots give a review of all actually integrated output variables.

Note:
The range of output monitors and plots that display important variables of colony status as well as foraging and Varroa dynamics can be extended as you like, by programming just a few additional code lines and adding advanced Netlogo output plots and monitors to the interface tab.

'Output monitor'
The integrated monitors of the BEEHAVE user interface display the most important variables of the actual colony composition and stores, foraging trips and infestation by Varroa mites on that day (table 1.7). At the beginning of the model run starting on the 1 st January with an initial colony size of 10,000 worker bees, all these bees are defined as foragers with an average age of 130 days. Hence, the total number of workers (# workers) is identical to # foragers, until the first cohort of in-hive bees emerges ( fig. 1.30). Note: The abrupt, distinct decrease in the average age (mean age foragers) illustrates the extinction of the long-living foragers (surviving winter bees) and thus the emergence of newly foragers attaining the age of first foraging (aff).  In addition, forager squadrons are divided into healthy, infected as adults and infected as pupae, whereby all these squadrons are moved to the right when they ageing. Moreover, successful foragers assume the colour of their last visited flower patch, while scouts and inactive bees are shown in grey.
The current time step of the simulation run is shown by a tick counter. Hence, the model proceeds in daily time steps represented by ticks.
However, the BEEHAVE world provides several additional symbols that illustrate the foraging opportunities, the current state of the brood, the infestation of immature stages by mites and also the colony stores as well as applied beekeeping treatments.
Thus, the positioning of the 'red' and 'green flower' to the left and right of the hive represents the distances of the nectar-and pollen-providing flower patches to the colony. If the mapped flowers are depicted as withered no further nectar and pollen is supplied by these patches.
In addition, local weather conditions strongly influence flight activities of scouting and foraging bees on that day. Hence, the following weather symbols indicate the duration of the foraging period on that day.
'Cloud' indicates that the local weather conditions do not support foraging trips on that day. The 'exclamation mark' is shown on days, when then foraging probability is set to zero, which may happen even if the weather is suitable for foraging.
Death of the brood due to lack of pollen or nurse bees is a further important indicator of the current colony state and is displayed as a symbol.
Moreover, beekeeping options like honey harvesting, feeding bees, stocking up colony stores, merging weak colonies and treatment of invaded Varroa mites (not included in the Default setting) may also be very important drivers of colony dynamics.
Thus, a 'Beekeeper' represents the onset of at least one of the following beekeeping options. Now, honey or pollen stores are added to the colony store ('food' or 'pollen grain').
Additionally, the beekeeping symbol 'Honey' indicates the harvest of honey during the predefined harvesting period regarding the defined threshold of honey stores and the amount of honey [kg] harvested.
Bees can be added to a weak colony in autumn if the current colony size is below a defined threshold value ('2 hives'). Furthermore, treatments of Varroa mites take place in autumn ('Varroa').
In the end, if the colony size falls below the threshold of the critical colony size during winter, the honey stores are below zero or all bees have died, then the death of the honeybee colony is displayed by the 'skull'.
In this section we will demonstrate your first pre-defined BEEHAVE scenarios with step-bystep instructions.
2. Type in the input field 'RAND_SEED' 8 and press the Default button.
However, the food sources ('red' and 'green flower') providing nectar and pollen have different distances to the hive. The withered flowers indicate that the existing patches offer no nectar and pollen to foraging bees at the moment whereas the arrows display the initial stores of pollen (left: 100 g) and honey (right: 25 kg) of the colony. 10,000 long-living winter bees, here referred to as foragers, constitute the initial colony size (figure 2.1).

Note:
The choice of this pre-defined scenario resets all input variables to their default values (table 2.1) and requires 2 feeders ("red" and "green" patch), seasonal food flow as well as local weather data at Rothamsted from 2009. In the Default scenario beekeeping options and Varroa infestation are turned off.  5. Now, let's take a closer look at those 'output monitors' that offer an easy way to see how the most important properties of the system change as the model runs. Firstly, in order to recall the definition of the several displayed variables click the right mouse button on any 'output monitor' and then click Edit. The dialog window of the monitor appears; here, the monitors of four outputs, referred to their input field "Display Name", are shown: The monitor windows show the displayed name of the output variable and the instruction for computing their daily value (called 'Reporter' in NetLogo). Thus, the total daily number of workers (# workers) sums up the total number of in-hive bees and foragers (totalIHbees + totalForagers) on that day. The total amount of brood (# brood) is composed of the total number of worker and drone brood cells (totalWorkerAndDroneBrood) on the current day.
The daily workload of the colony is a function of the total amount of brood cells and the total number of working bees in association with the contribution of foraging to nursing bees and the ratio of brood cells to nursing bees (work load).
6. Now, start the model run by clicking the 1 Day button. On the top left of the interface the in-hive cohorts and forager squadrons, and several 'output monitors' illustrating the current colony status are displayed. Moreover, four worker eggs are laid by the queen (blue bar), but no age-cohorts of larvae (yellow bar) and pupae (brown bar) are developed during the current day. Thus, the daily total colony's amount of immature stages is only composed of four brood cells (# brood). The actual relative amount of brood to maximum nursing force (work load) amounts to zero on that day. Compared to the day before the stores of honey and pollen decreases (red arrows: left: pollen; right: honey). The current pollen store only exceeds 34.1 % of the ideal pollen store (% pollen store). In addition, the actual weather conditions allow no foraging trips as represented by the 'Cloud' (figure 2.2).
Nevertheless, one scouting forager squadron leaves the hive (grey bee). But only one patch 'red flower' provides nectar on that day. In contrast, the withered 'green flower' symbolizes no pollen and nectar provision by the green patch.
Thus, we will check the food availability of the existing patches and also the current visits at the flower patches by clicking the advanced inputs (bottom right of interface): and .
Now, click of the 'Command Center' to enlarge its window. The following information of the available patches is printed out by the 'Command Center': Hence, only the 'red flower' patch offers very low amounts of nectar. Moreover, neither flower patch is visited by any forager squadron.  8. Click the 1 Day button several times on the trot until the 'output monitor' Julian Day displays day '7' (tick counter). Now, press the advanced button active patches.
The egg-laying rate ('Generic plot 1') and also the number of brood cells (# brood) increases.
Furthermore, the first age-cohort of the larval developmental stage (yellow bars) hatches, but no pupae (brown bars) and in-hive bee cohorts (orange bars) emerge. At this time point no drone brood cells have developed (bars of age-cohorts on the left side of the panel). The 'output monitors' display the increase in the mean age of foragers, the decrease in the age of first foraging (aff) and a lack of pollen stores (% pollen store). In addition, the amount of nectar provided by the 'red flower' patch increases, but the 'green flower' patch still does not offer nectar rewards ('Command Center', 'Generic plot 7'). Insufficiently good weather conditions symbolized by the 'Cloud' inhibit foraging trips on that day. Therefore, honey and pollen stores decrease compared to the day before (red arrows).  10. Now, we want to inspect several age-cohorts of the developmental stages of the inhive bees by using the agent monitors of Netlogo. Click the right mouse button on the first egg age-cohort and pupae age-cohort (the circle with a label on it; it is next to the bar indicating the size of this cohort). Next, choose egg-cohort 649 (and pupaecohort 653) and click inspect egg-cohort 649 (and inspect pupae-cohort 653).
The window of the egg-cohort appears and displays the cohort ID (who), the ploidy level and age of the egg-cohort as well as the number of individuals in this age-cohort (number). Now close these monitors.  fig. 2.7). Nevertheless, 'Generic plot 8' also illustrates the loads of returning foragers from day 1 to day 127 of simulation run; thereby crops of pollen and nectar collected by active forager squadrons were brought to the hive within the last days.  Check the colony dynamics (panel, 'Generic plot 2') and stores as well as foraging activities.

Click the
Note, the number of foragers' (# foragers) increases, but at the same time the mean age of foragers decreases. The current foragers emerge from the newly developed in-hive bees attaining their age of first foraging (aff). Thus, forager squadrons are not mainly constituted of long-living foragers (winter bees) any more. Moreover, no foraging activities and thus no nectar and pollen visits are performed on that day due to bad weather conditions as displayed by the 'Could' and 'output monitor' trips/ h ( fig. 2.8). 13. Before continuing our simulation example, we select the option 'All visits' using the 'Foraging map' chooser in order to visualize all foraging trips of all active forager squadrons at the available patches. Now, click the 1 Month button.
Also, we want to take a look at the offered pollen and nectar amounts as well as the number of visits at the available flower patches on that day. Thus, click visited patches.
Check the development of age-based cohorts and the current stores of the colony.
Furthermore, compare nectar and pollen availability, and also the number of visits at the existing patches ( fig. 2.9).  14. Click 1 Month. Now, the 'output monitor' Julian Day shows day '217`of model run.
Again, consider the dynamics of the age-cohorts and the seasonal progress of egglaying ( fig. 2.10).
You will observe that the actual intra-colonial dynamics are characterized by a distinct decrease in the queen's egg-laying activity ('Generic plot 1') as well as in the amount of brood ('Generic plot 2'), and a considerably increase in the number of in-hive bees emerging from the previous immature stages ('Generic plot 2'). Therefore, the workload resulting from the relative amount of brood to the maximum nursing force decreases as displayed by 'Generic plot 3' of figure 2.10. As illustrated by 'Generic plot 5' the age of first foraging increases and attains an actual value of 24 days (aff). Moreover, the pollen ('Generic plot 6') and nectar ('Generic plot 7') availability of the 'red' patch providing nectar and pollen early in the season strongly decreases, while the 'green' patch offers large amounts of rewards. In total 6746 foraging trips per hour (trips/h) due to suitable weather conditions, that allow a foraging duration of 5.9 hours ('Sun'), are performed on that day. 15. Next, click the 1 Month button twice. Well, we should attain day '277' of our simulation run. Apply the advanced buttons show Patches and visited patches in order to recheck the provided nectar rewards, energetic efficiency (EEF) and forager visits of the available patches at the current day ( fig. 2.11).
When you take a look at the intra-colonial dynamics you see a distinct decrease in the individual number of all age-cohorts as illustrated by 'Generic Plot 2' of figure 2.11.
Although the current weather conditions are suitable for a foraging period of 4.9 hours on that day ('Sun'), no foraging trips are performed (monitors: Nectar visits, Pollen visits and trips/h, 'Command Center'), because the foraging probability is set to zero as illustrated by the 'exclamation mark'. Furthermore, only the 'green' patch flowering later in the season actually offers nectar and pollen rewards ('Command Center').     Note: You can also write 1095 in the 'input field' X_Days and click run X days. Figure 2.14: Colony and foraging dynamics of three consecutive years.

The Varroa Scenario
1. Type in the input field 'RAND_SEED' 8 and press the varroa button in order to initialize the pre-defined varroa scenario. 4. When the pre-defined varroa scenario has been initialized, several 'output monitors' display specific variables of initial mite infestation ( fig. 2.15). Thus, in total 20 phoretic mites (total mites) initially infest the honeybee colony, 50 percent of these are infected with the deformed wing virus (rate healthy mites).
right mouse button on any mite-specific 'output monitor' and then click on Edit. The dialog window of the monitor appears. Here, the monitors of five outputs, referred to by the input field "Display Name", are shown. Now, take a closer look at the reporters by which the values of mite-specific variables are daily computed: 5. Next, click Run.
Now, this pre-defined simulation runs until the death of the colony illustrated by the 'skull'.
Subsequently, the day of colony death, colony sizes and honey stores of the several years are written to the 'output box'. Furthermore, for the day of colony's death the number of healthy and infected phoretic mites (phoretic mites healthy, phoretic mites INFECTED) is printed out (figure 2.16). Furthermore, the number of mites invading the honeybee colony and proportion of infected mites increases over time, but you can also see that the mite dynamics strongly depend on the honeybee colony dynamics ('Generic plot 4 & 7').

Figure 2.17:
The temporal progress of egg-laying, structure of immature and adult stages of worker bees and drones, healthy and DWV-transmitting mites and mite drops, honey and pollen stores of the hive as well as the age of first foraging and lifespan of worker bees is presented over several years until the death of the colony. Now, we conduct the pre-defined varroa scenario once more, but short step-by-step instructions offer the possibility to explore colony and mite dynamics as well as foraging activities at different time steps of simulation run.
6. Thus, enter '127' into the input field 'X_days' and start the simulation run by clicking run X days.
Take a look at the age-cohorts of the colony, the actual number of healthy and DWV infected phoretic mites invading the colony and the current proportions of healthy and diseased foragers.
In total 26 Varroa mites actually infest the honeybee colony. However, 25 Varroa mites invade unsealed cells of the last larval stage on that day (output monitor: mites in cells).
Moreover, one phoretic mites is healthy (phoretic mites healthy) and no mites are infected by  The following windows appear: The windows of the larvae and pupae cohort display the ID (who), age and ploidy (1: drone, 2: worker) of the actual immature cohort. In addition, the ID of the mite organizer which invaded this actual brood cohort is displayed (invaded by mite organizer id 8. Click the right mouse button on the first Varroa sign in the following way in order to inspect the first mite organizer: Choose miteorganiser1268 and click inspect miteorganiser1268.
The following window appears: The window of the mite organizer shows how many brood cells of one worker or drone brood cohort are infested with zero, one, two or more mites whereas the maximum number of mites invading one brood cell depends on the chosen model of mite reproduction (e.g. 'Martin': four mites). However, the ID (who) of the mite organizer (1268)  year is printed underneath the panel of the age-cohorts.
On that day in total 526 Varroa mites invade the honeybee colony. 465 mites have entered brood cells just prior to capping and reproduce within these sealed drone and worker cells.
The remaining mites live attached to adult bees ( fig. 2.20). Moreover, the proportion of infected mites increases as illustrated by 'Generic plot 7', and no foragers were infected by DWV-transferring mites as adults (rate healthy foragers; 'Generic plot 5'). 11. In the input field 'X days' enter '90' and click run X days.
The tick counter should display 582.
The panel of age-cohorts shows that worker and drone brood cells are infested by mites.
Moreover, adult drones and workers have been infected as pupae as well as adult bees. Thus, the number of mites invading the colony increases to 4004 on that day. The majority of foragers are healthy (rate healthy foragers: 0.855), but some foragers have been infected as pupae (diseased) or adults (DWV-carriers) as illustrated by 'Generic plot 5'. Inspect one mite organizer and some age-cohorts of larvae, pupae and in-hive workers by clicking the right mouse button on any age-cohort or mite organizer.
Then select the marked mite organizer or age-cohort and click inspect … cohort/ mite organizer … (ID).
Now, the rate of healthy mites is very low (0.01336) and 142 mites actually fall from the combs (mite fall). Furthermore, the majority of in-hive workers have been infected as adults as illustrated by figure 2.22. For instance, the in-hive bee cohort 9528 consists of 92 in-hive workers including 69 workers infected as pupae, 7 workers attacked by mites as adults and only 16 healthy workers ( fig. 2.23 left). Moreover, the window of the miteorganiser 9557  13. Now enter '150' into 'X_days' and click run X days.
The simulation run attains day 1007 and the majority of active foragers are infected by the deformed wing virus transmitted by infected mites. Thus, foragers of the colony are either carriers of the virus or already diseased ('Generic plot 5'). Furthermore, the queens's egglaying rate, the amount of worker and drone brood and the total number of adult workers and drones decreases over the simulation period ('Generic plot 1, 2 & 6'), whereas the mite infestation distinctly increases ('Generic plot 4'). Also, take a look at the current number of healthy and infected phoretic mites, the number of mites infesting brood cells of the colony and the mite fall ( fig. 2.24).
15. Subsequently, perform a second Varroa simulation experiment by reducing the distances of the foraging sources to the hive. For this purpose, enter the following distances of patches to the hive: '1000' and '300' metres into the input boxes DISTANCE_R and DISTANCE_G.
Next, click Setup in order to initialize these new parameter values. Now, click Run and wait until the colony died. Inspect the outputs of the model run.
What has happened?
Note: You can also change the virus type of infected phoretic mites, the number of initial healthy and virus-infected mites as well as the mite reproduction model in order to examine the effect of different Varroa infestation conditions on colony survival.
Furthermore, you may perform simulation experiments combining Varroa infestation and different foraging conditions. investigations of e.g. factors affecting performance and foraging behaviour of honeybee colonies briefly and quickly. As or real experiments, the 'design of an experiment' have to be specified. Thus, we suggest that you give a complete and self-explanatory title and description of your simulation experiments performed with BEEHAVE by answering the following questions and briefly document your assumptions and chosen parameters: -What is the purpose of the simulation experiment?
-Specify the settings of the simulation experiment: Which model versions, parameter values, initial values of state variables were used?
Which assumptions were applied?
Define the ranges and values over those specific parameters were varied! Which references (literature, field data, and expert opinion) for used parameters were turned into account?
-What currencies (output variables) did you measured?
-What is the time horizon for each simulation?
-How many repetitions were run for each simulation?
-What processes in the model are stochastic?
Probably, you like to create a short overview of your simulation experiments and may use the following suggested box head (table 3.2).