Simulating the Effects of Pesticides on Honey Bee ( Apis mellifera L.) Colonies with BeePop+

: The US Environmental Protection Agency (USEPA) employs a tiered process for assessing risks of pesticides to bees. The model discussed in this paper focuses on honey bees ( Apis mellifera L.). If risks to honey bees are identiﬁed at the ﬁrst tier based on exposure and toxicity data for individual adult and larval honey bees, then effects are evaluated in higher-tier studies using honey bee colonies. Colony-level studies require large amounts of resources (to conduct and review) and can yield data complicated by the inherent variability of colonies, which are inﬂuenced by factors that cannot readily be controlled, including weather, pests, diseases, available forage, and bee management practices. To better interpret these data, the USEPA and the US Department of Agriculture (USDA) developed a simulation model, BeePop+, that assesses potential honey bee colony-level effects of pesticides. Here, we describe this model using the population model guidance, use, interpretation, and development for ecological risk assessment (Pop-GUIDE) framework, which is a conceptual framework for the development and evaluation of population models. Within the context of Pop-GUIDE, BeePop+ is considered a “realistic-precise” model and reﬂects the inherent variability of colony response to pesticide exposure by simulating many outcomes. This model meets the desired features needed for use in pesticide risk assessments as its required data inputs are typically available, it is applicable to different US locations, and the outputs are both relevant to USEPA’s protection goals for honey bees and are consistent with the outcomes of empirical studies. This model has also been evaluated using available empirical colony-level data; however, additional evaluation with other studies may still be done in the future prior to completing implementation.


Introduction
Bees (Hymenoptera: Apoidea: Anthophila), especially honey bees (Apis mellifera L.), are essential pollinators of farmed and natural landscapes, without which the productivity of around 80% of food crops would be seriously compromised [1]. The economic dependence of agriculture on pollination services is significant (USD 14.2-23.8 billion), but the higher-order economic dependence of industrial sectors fueled by crop production is also substantial (USD 10.3−34.0 billion) [2,3]. The value of crops produced by managed honey bee pollination influences multiple socioeconomic sectors, generating jobs and revenue for rural areas and for numerous industrial sectors through equipment and machinery manufacturers, agrochemical companies, food processing, and shipping to name just a few. Almond production alone requires more than two million colonies to pollinate the nearly 1 million acres (4000 km 2 ) of bearing almond trees [4]. The almond crop is worth USD 2.2 billion and adds an estimated USD 21.5 billion to the California economy [5,6] and accounts for 80% of the global almond market. Honey bee-pollinated crops support export markets that help balance trade deficits [7]. From the perspective of human nutrition, honey bee-pollinated crops are essential to human health [8,9]. Globally, there is concern over declines in honey bee health based on annual colony losses [10][11][12][13]. The losses are above the rate that beekeepers indicate is sustainable [14]. Growers that rely upon bees for pollination services are also concerned about subsequent reductions in yields. Declines in honey bee health have been attributed to four major factors: pesticides, pests (e.g., Varroa mites, Varroa destructor Anderson and Trueman), disease, nutrition, and genetics [15]. While no single factor has been identified as a sole cause of the losses, there is a recognition that interactions among the factors are reducing colony survival [10]. Both the United States Department of Agriculture (USDA) and the United States Environmental Protection Agency (USEPA) have been engaged in activities and research to better understand causes of colony losses and how to mitigate them and improve honey bee health [16]. The USEPA assesses risks of pesticides to honey bees using a tiered approach [17] ( Figure 1). The first tier relies upon laboratory toxicity data quantifying pesticide effects to individual adult and larvae from acute and chronic exposure. Tier I exposure estimates are based on default conservative assumptions for direct spray to forager bees and exposures via consumption of residues in pollen and nectar that are incorporated into the BeeREX model [17,18]. If risk is identified at Tier I and additional information is needed to define the risks to bees, empirical exposure data (i.e., measured residues in pollen and nectar) may be used to refine concentrations in nectar and pollen. At Tier II, toxicity studies involving honey bee colonies may be used to understand how effects on individual bees observed during laboratory-based studies manifest at the colony level. Tier II toxicity studies generally involve small (nucleus) or full-size colonies kept under "semi-field" conditions, where exposure to a pesticide is controlled (i.e., through confinement of colonies to a treated area using a tunnel or through feeding colonies known concentrations of pesticide). Tier III effects studies involve placement of full-sized colonies in treated fields where a pesticide is applied according to the label, and bees forage freely on the crop and other plants in adjacent areas. Tier III studies are intended to be representative of "real world" exposure conditions; however, they may be limited in their utility as (1) studies are influenced by test site factors (e.g., landscape, weather) and inter-colony variability, (2) the results may lack relevancy to other crops for which the pesticide is used (e.g., a study with cotton may not be representative for almonds, exposures from seed treatments may differ from foliar sprays), and (3) the resources required to appropriately replicate and review large-scale complex studies may be high. Another challenge with Tier II and III studies is the presence of a high degree of variability in measured endpoints, which often makes interpretation difficult (i.e., it may be difficult to isolate effects that are attributable to pesticide exposures). In addition, test conditions may result in stress to bee colonies that can impact measured endpoints (e.g., confinement stress in Tier II enclosure studies).
As a supplement to the current tiered risk assessment framework, the USEPA is investigating the use of a simulation model that can efficiently and reliably account for honey bee colony dynamics in response to pesticides [19]. The USEPA and the USDA have added pesticide exposure/effects capabilities to an existing model (BeePop/VarroaPop) to develop an updated simulation model (BeePop+) that assesses potential colony-level effects of pesticides. This model was designed to either replace the need for a higher tier empirical study or help design or interpret results from colony-level studies. In addition to pesticide effects, the model accounts for other colony stressors, including Varroa and weather. We use the population model guidance, use, interpretation, and development for ecological risk assessment [20] as a framework for documenting key components of BeePop+ development, implementation, and documentation. The Pop-GUIDE serves as a transparent, stepwise process that guides population model development while being consistent with risk assessment objectives and data availability constraints. We also use the Pop-GUIDE framework to describe potential applications of BeePop+. As a supplement to the current tiered risk assessment framework, the USEPA is investigating the use of a simulation model that can efficiently and reliably account for honey bee colony dynamics in response to pesticides [19]. The USEPA and the USDA have added pesticide exposure/effects capabilities to an existing model (BeePop/VarroaPop) to develop an updated simulation model (BeePop+) that assesses potential colony-level effects of pesticides. This model was designed to either replace the need for a higher tier empirical study or help design or interpret results from colony-level studies. In addition to pesticide effects, the model accounts for other colony stressors, including Varroa and weather. We use the population model guidance, use, interpretation, and development for ecological risk assessment [20] as a framework for documenting key components of BeePop+ development, implementation, and documentation. The Pop-GUIDE serves as a transparent, stepwise process that guides population model development while being consistent with risk assessment objectives and data availability constraints. We also use the Pop-GUIDE framework to describe potential applications of BeePop+.

Materials and Methods
The Pop-GUIDE process consists of five linked phases that are followed in a stepwise manner (see Figure 1 of [20]). The BeePop+ model described here is based on previously developed modelling platforms, and most, but not all, elements of these phases already existed (i.e., BeePop+ represents additions to the model VARROAPOP [21], which was based on BeePop [22]). More detail on model history and documentation is provided in the Supplementary Materials. Therefore, the Pop-GUIDE framework was used to structure the presentation of the novel elements of this continuing effort. This included documenting the implementation of the pesticide component in the context of the honey bee colony conceptual model through a process focused on identifying model objectives, compiling relevant parameterization data, and considering alternative algorithms to represent key processes. This serves as a basis for the current model version and will also help guide further development efforts of the colony model.
The first phase (i.e., identifying model objectives) considers regulatory statutes as a context for evaluating trade-offs associated with the required degree of realism and precision from the model output. We leverage the Pop-GUIDE decision tree [20], with some modifications, to determine these requirements. Desired features of a colony simulation model for regulatory purposes include: (1) accounting for major factors that determine Tiered approach used by USEPA to assess pesticide risks to honey bees. C = contact exposure; O = oral exposure.

Materials and Methods
The Pop-GUIDE process consists of five linked phases that are followed in a stepwise manner (see Figure 1 of [20]). The BeePop+ model described here is based on previously developed modelling platforms, and most, but not all, elements of these phases already existed (i.e., BeePop+ represents additions to the model VARROAPOP [21], which was based on BeePop [22]). More detail on model history and documentation is provided in the Supplementary Materials. Therefore, the Pop-GUIDE framework was used to structure the presentation of the novel elements of this continuing effort. This included documenting the implementation of the pesticide component in the context of the honey bee colony conceptual model through a process focused on identifying model objectives, compiling relevant parameterization data, and considering alternative algorithms to represent key processes. This serves as a basis for the current model version and will also help guide further development efforts of the colony model.
The first phase (i.e., identifying model objectives) considers regulatory statutes as a context for evaluating trade-offs associated with the required degree of realism and precision from the model output. We leverage the Pop-GUIDE decision tree [20], with some modifications, to determine these requirements. Desired features of a colony simulation model for regulatory purposes include: (1) accounting for major factors that determine colony dynamics; (2) quantitatively linking various measurement endpoints to assessment endpoints; (3) readily parameterizable, using existing biological and pesticide-specific data; (4) accounting for geographic variation in pesticide use pattern, crops and climate, and (5) being scientifically defensible, transparent, well-documented, and publicly available. Surveys of the available literature indicated that several colony models were available that met some, but not all, requirements [19,23]. Although some candidate models included many of the desired features, none implement the primary feature of quantifying pesticide exposure and effects on the colony. Ultimately, an existing USDA model originally implemented as BEEPOP [22], and later modified as VARROAPOP [21], was selected for further development. The VARROAPOP model met most of the desired features identified by USEPA and has been subsequently modified to quantify pesticide exposure and effects at the colony level [24,25]. The modified version of VARROAPOP, referred to as BeePop+, is the model discussed in this paper.
The second phase consists of compiling pertinent biological, chemical, and environmental data. The last decade has seen a steady increase in the publication of relevant data on exposure and effects of pesticides on honey bees. We organized available types of data into tables as suggested by Raimondo et al. [20], structuring information based on organism, population, spatial, and habitat characteristics. As a regulatory agency, the USEPA also has authority to request data for pesticide registration decisions, so we also focus on those data from unpublished studies submitted in support of pesticide registration.
The third phase includes a series of steps that elaborate on central modeling concepts in addition to exposure and effects, including: life history representation; organism-level processes that relate to growth and reproduction; population influences such as spatial factors; density-dependence and behavior; and other factors that range from diet to interspecies interactions. This phase is central to identifying the core algorithms to construct the population model and associated exposure and effects models.
The fourth phase brings these elements together into a unified conceptual model that provides a graphical and textual summary of the population compartments, with important functions and other linkages between compartments clearly identified. The core of this conceptual model is the life history representation of the honey bee colony. The level of detail of the conceptual model is heavily influenced by the levels of realism and precision identified in the first phase, in combination with the available data for parameterization from the second phase.
Finally, the fifth phase includes model implementation and evaluation. This step includes model parameterization and calibration with subsequent performance evaluation based on output behavior and statistical fit to observations. We summarize two studies where successively improved versions of BeePop+ have been implemented via a broad sensitivity analysis [24] and fit of model results to field data from an empirical study where honey bee colonies were exposed to the insecticide clothianidin through diet [25]. These studies allowed for simulations of colony demography and mortality from ingestion of pesticide-contaminated food under a range of conditions. These performance evaluations qualitatively evaluate the realism of the colony population dynamics (with and without a pesticide stressor) and quantitatively fit observational data that includes colony population size, structure, and responses over time, including responses to different levels of exposure. Specification of objectives are necessary to inform the needed levels of realism, complexity and acceptable levels of uncertainty required by the model. An important early step in this process is determining the needed complexity for the model. The USEPA uses a tiered weight-of-evidence approach to identify data collection needs and evaluate ecological risks from pesticides (more detail in Section 1). For BeePop+, model objectives are largely informed by the relevant regulatory statute, i.e., the Federal Insecticide Fungicide Rodenticide Act (FIFRA) and USEPA's pesticide risk assessment method for bees [17] and by its location within the tiered bee risk assessment framework. Therefore, the complexity of the model is heavily influenced by the USEPA objective to utilize the model as part of a weight-of-evidence analysis to assess pesticide risks to honey bees. BeePop+ was designed to help translate individual-level exposure and effects information into potential colony level effects using data that are typically available. It is possible that this model could be sufficient to address risk assessment questions without the need for a higher tier empirical study. Alternatively, this model may be useful in developing specific hypotheses that can be integrated into higher tier empirical study designs.

Results and Discussion
BeePop+ is designed to account for major factors that determine colony dynamics in a wide variety of locations across the US. USEPA's ecological risk assessments (ERAs) for pesticide registration are generally conducted at a national scale, just as honey bee colonies are widely used for pollination services and honey production throughout the US. It is well known that beekeepers move colonies across the country to meet the pollination needs of various crops (e.g., almonds in February, apples in May). As weather is known to impact honey bee colony dynamics [22], and bee-attractive crops are often limited in their geography, this model is designed to account for differences in use patterns for a pesticide active ingredient across crops, climate, and landscape.
To be useful to risk assessors, model parameterization needs to be derived from data that are readily available to USEPA (e.g., acute and chronic toxicity data for individual bees) and quantitatively links effects across various measurement endpoints to assessment endpoints of regulatory interest and ultimately to protection goals. Model outputs include indicators of colony size such as number of adults, number of individuals in different brood stages (i.e., eggs, larvae, pupae), colony survival, and production of colony products (e.g., honey).
Our use of the Pop-GUIDE decision tree (see Figure 2 of [20]) to determine trade-offs among generality, realism, and precision led us to describe this approach as 'realisticprecise'. Although this model does not simulate an endangered or threatened species (as discussed in the decision tree [20]), the model does represent a single species (i.e., the Western honey bee) that aligns with the protection goals of the USEPA method for assessing risks of pesticides to bees (discussed further below). Additionally, due to its position in the tiered assessment process, it is necessary to account for variability that may be attributed to different locations relevant to pesticide use sites (e.g., weather, application timing); therefore, a greater level of realism and precision is needed. This model may be broadly applied to the US and was not designed for a specific type of location; however, to simulate colony responses it requires some locally specific data, such as crop phenology and weather and pesticide use patterns.
One limitation of the model is that it is assumed that bees forage completely on pesticide treated areas. Honey bees forage over wide distances and collect pollen and nectar from a variety of plants. Although honey bees do not necessarily forage 100% on a single crop, one generality of this model is that they do so. Model realism could increase much more if landscape characteristics (e.g., other habitats within foraging range of colonies) were considered; however, this was seen to be overly complex when considering the utility of this model in a Tier II assessment. The proportion of pollen and nectar obtained from different plants varies widely based on time and landscape. A simplifying assumption of the model is that bees consume 100% of their pollen and nectar from a simulated crop. This assumption can be evaluated by calculating the amount of forage needed to result in a colony level impact and comparing that to available data documenting the proportion of pollen collected from single crops. Although one could argue that this model may be more general because of the simplified assumption related to landscape and foraging, it is still considered realistic because it can account for location and species-specific features.

How Model Outputs Link to Measurement Endpoints and Protection Goals
The USEPA would use the model to determine whether a proposed use of a pesticide impacts identified protection goals that include provision of pollination services, production of colony products, and preservation of bee diversity (Table 1) [17]. These protection goals are linked to assessment endpoints that describe how the goal will be considered and measurement endpoints that quantify effects levels in risk assessments [26]. BeePop+ can estimate multiple measurement endpoints considered to be of regulatory relevance, including widely recognized measures of colony health (e.g., abundance of eggs, larvae, pupae, adults) and hive quality characteristics (i.e., honey, pollen stores). Additionally, the model can represent the variability in colony responses to pesticides via Monte Carlo simulation of single or multiple colonies with varying biological parameters (e.g., queen egg-laying rate, worker lifespan) [24,25]. Therefore, the model allows for the assessment of pesticide impacts to measurement endpoints across different honey bee colonies while accounting for stochastic environmental and biological conditions. Since this model is focused on honey bee colonies, it does not directly address the protection goal related to preservation of bee diversity.

. Temporal and Spatial Considerations
Several major spatial and temporal considerations exist that are important for pesticide risk assessments for bees. Honey bees have pronounced variability in seasonal activity patterns that can influence exposure and susceptibility to pesticides in the environment. They are additionally influenced by weather and temperature, with colonies in different geographic locations having different activity patterns. Stressor dynamics also vary, and the timing of pesticide applications is generally dependent upon interactions between crop production cycles, environmental conditions, associated crop pests, and disease pressure. Finally, an important management factor is that the seasonal transport of bees across the US in response to agricultural pollination needs or to support honey production influences the timing and types of pesticides to which a colony is exposed.
In temperate climates, honey bees remain inside the colony for the duration of cold winter periods (overwinter) and will only forage outside of the colony once certain ambient temperature and photoperiod thresholds are met in the spring [22]. The resources that support a colony during the overwintering phase are accumulated during the previous year. Within a single colony, the population structure, population size, honey stores, and pollen stores change over the course of the year with resource availability and weather conditions. In general, in early spring, the number of individuals within a colony increases exponentially. The size of the colony generally levels off in the late summer and then decreases in the fall in preparation for winter and in response to declining food resources. The timing of pesticide exposure within these seasonal activity patterns may significantly alter impacts on the colony. For example, exposure that causes mortality to foragers during periods of high foraging activity may be particularly harmful, as these periods are critical for provisioning a colony to survive overwintering. Alternately, exposure to a pesticide that causes high larval mortality may be most harmful in the spring, when large numbers of larvae are being reared to rebuild the adult bee population that declines in the winter and when foraging resumes in the spring.
Honey bee colonies respond to weather conditions, with wind velocity, temperature, and rainfall thresholds generally accepted as restricting the ability of honey bees to collect nectar and pollen [22,27,28]. These activity limits can influence colony resources due to a lack of foraging activity and may also result in a lower egg-laying rate of the queen. These influences are most clearly observed in early spring, when colonies in warmer southern locations typically have greater numbers of bees compared to colonies in colder locations.
During agricultural growing periods, the interaction of temporal, spatial, and habitat factors influence the foraging behavior and success of active honey bee colonies, along with their subsequent exposure potential to pesticides. Pollen and nectar resource availability depends upon the plant species assemblage within the foraging range, the blooming phenology of different plants in the assemblage, and the proximity of the colony to those plants over the course of the season. Temporal "forage gaps" in the availability of blooming plants and long flight distances to suitable foraging resources can result in negative impacts on population dynamics [29]. In agricultural landscapes where honey bees are used to pollinate crops, both pesticide application timing and the plant phenology affect bee colony dynamics by influencing pesticide exposures. If spray applications occur when bees are foraging, bees are likely to receive the greatest possible exposure to a pesticide. Pesticides applied prior to blooming may not result in direct contact with foraging bees, but still may generate significant bee exposures if they are persistent on plant surfaces or are systemic and translocated to pollen and nectar [30,31].
Furthermore, the presence of honey bee colonies on agricultural landscapes is often dependent on the transport of colonies for pollination services, with colonies often used for the pollination of multiple crop types within a single season. As different growing systems may employ different types of pesticides, and honey bee colonies are moved to pollinate different crops, exposures throughout a season can vary [32]. Thus, pesticide exposure conditions for honey bees exist at the intersections of seasonal colony activity and the type and timing of pesticide applications.
In addition, honey bees are exposed to multiple stressors besides pesticides. Disease and pest pressures are influenced by increased intraregional colony transport for pollination services and intercontinental trade of bee and bee products, with important spatial and temporal fluctuations in stressor levels. For example, infectious spores of the gut parasite Nosema sp. decrease in number from spring to summer but peak in the late fall and winter when the worker bee population is lower and confined inside the colony. Alternatively, the parasitic Varroa mite reproduces on developing broods and steadily increases throughout the spring and summer, peaking in the late summer [33]. These parasites and the viruses they vector can act as co-stressors during exposure to pesticides and may interact to influence the susceptibility of a colony to pesticides [34,35].

Assumptions and Sources of Uncertainty
A model should have a level of complexity and uncertainty that matches the available data and the needs of the assessment. The BeePop+ model is relatively complex according to the Pop-GUIDEcriteria in [20] since it accounts for dynamics within a honey bee colony and responses to weather. Uncertainty exists in the output of the BeePop+ model due to natural variability in bee and colony characteristics and uncertainty in the value of model parameters. Key factors of colony dynamics, such as the queen egg-laying rate, forager lifespan and efficiency, and pesticide sensitivity, are known to vary significantly across A. mellifera subspecies and even among colonies from the same genetic line. The best fitting model parameter values for these and other colony characteristics may also be poorly defined due to a scarcity of prior evidence, inter-study variability, and other sources of random measurement error. We can measure and account for uncertainty by identifying key model parameters via global sensitivity analysis [24], allowing these sensitive parameters to vary probabilistically, and inferring their distribution from empirical data using a Bayesian approach [25].
Because honey bee colonies are complex systems, simplifying assumptions must be made in the model structure, adding another source of uncertainty. For example, foragers gather resources from a complex landscape that is affected by plant cultivation, population dynamics and phenology, human activity, predation, mortality rate of foragers, and daily weather patterns, resulting in changes to the amount and quality of pollen and nectar gathered. A model may simplify the foraging process by assuming an average daily pollen and nectar load on days where weather is suitable for flight. Alternately, a simplified landscape might be considered where resource patches are defined, and foraging success is calculated based on factors like patch size, quality, and distance from the colony. In either case, it is important to consider how these simplifying assumptions contribute to uncertainty by examining model performance across a range of landscape conditions and complexities.
Monte Carlo simulation can be used to address variability in pesticide response, as it allows for variation of key parameters (e.g., queen's egg-laying rate, adult forager life span) and resulting influence of a pesticide across exposed colonies. When the model is run with and without pesticide exposure, inherent variability across colonies can be accounted for. This allows for the evaluation of conditions that may make a colony susceptible to pesticide exposure and the proportion (or probability) of colonies that may be impacted. BeePop+ also can account for the impact of Varroa mite infestation and pesticide exposure on simulated colonies. Therefore, when considering multiple factors and accounting for natural variability, the model can account for conditions that may influence a colony's response to pesticide exposure, resulting in a more robust evaluation. This further supports the decision to select a realistic-precise design for the BeePop+ model.
The model can also account for different locations of crops by simulating different weather related to specific regions and bloom time windows. The model can also run with default exposure estimates or empirical concentration data for pollen and nectar. As previously discussed, there may be influences of landscape outside of a treated field on the magnitude of exposure for simulated colonies. At this time, the variability in floral resources outside of the treated field and their influence on the magnitude of exposure (i.e., influence on the proportion of pollen and nectar collected from the pesticide treated crop) is not accounted for quantitatively in the model. The model also does not account for nutritional stress on colonies, instead assuming that sufficient pollen and nectar are available to meet a simulated colony's needs. This was determined to be outside of the desired complexity of the model and availability of data on non-crop floral resources across time and space. As mentioned in Section 3.1, the model assumes that 100% of pollen and/or nectar is collected from the treated crop. Alternatively, this conservative assumption can be evaluated by inversely determining the proportion of pollen and/or nectar from the treated site that is needed to determine a given exposure or potential risk.

Phase 2: Data Compilation
There are many factors that may influence the likelihood and magnitude of the effects of pesticide exposures on honey bee colonies. The factors can be divided into groups that include honey bee colony dynamics, exposure and toxicity of the assessed pesticide, and spatiotemporal factors. Data are compiled to represent these different factors in the colony simulation model. Some data are represented by internal model assumptions and others are parameters that represent a given pesticide or location. As suggested by Schmolke et al. [36], the compilation of available data necessary for a minimal conceptual model can be efficiently done with the construction of species (Table 2) and exposure effects ( Table 3) tables. Also considered below are additional factors that may influence responses of colonies to pesticides, including environmental and biological stressors, landscape composition, and beekeeping practices (Table 4). These tables include information about specified characteristics, their associated ecological concepts, and how they are specified (e.g., as a parameter) or accounted for (e.g., as a simplifying assumption) in the model.

Phase 3: Decision Steps
The decision steps leverage the previously compiled data and determine how to represent core conceptual model components such as life history structure, organism-level processes, population factors, and external impacts while accounting for the effects of chemical exposures. These steps balance available data, alternative algorithm options, and the importance of each component to the overall model objectives. These objectives are based on designing the model to have the ability to translate individual-level toxicity and exposure information to predict colony-level effects and outcomes and to inform the design and interpretation of semi-field and full-field colony-level empirical studies submitted for pesticide registration.
The first decision step (i.e., representation of the life cycle structure) is depicted in Figure 2. The model is resolved daily, with tracking of population sizes of eggs, larvae, and different adult classes. The second decision step (i.e., organism-level processes for reproduction (egg-laying), maturation, and mortality) is addressed through the original implementation of BEEPOP [22]. Therefore, the core population model has stressor-independent inputs for life history stage growth, promotion, resource collection, and mortality (Table 2). However, additional model inputs were needed for the revised model to have utility for pesticide evaluations. Additional classes of input parameters for risk assessment include pesticide physicochemical properties, toxicity information, method-specific application rates, and crop phenology information (Table 3). Table 2. Data representing life history information for honey bee colonies.

Stage-level energetics [17]
Different stages (adult, larvae, drones, queens) have specific energetic needs that are met by the consumption of nectar/honey and pollen.
Stage-specific energetic needs are specified directly as model parameters represented by consumption rates.

Active season colony-level energetics [17]
Sufficient energy and pollen availabilities are required to raise brood Nectar and pollen consumption rates for brood stages are included as parameters.
Overwintering colony-level energetics [37][38][39] Sufficient energy reserves and adult members (due to thermoregulation needs) are required to overwinter Overwinter dynamics are simulated by significantly reduced foraging and egg-laying behavior, as dictated by daylength and temperature.
Colony growth * [21,22] Colony growth is determined by the rate of egg-laying by the queen, the mortality of each stage of bees in the colony, and the density of the colony.
Egg-laying rate is defined as queen strength parameter and is influenced by a function incorporating photoperiod, temperature, and colony size. Natural stage-based mortality rates are included directly as parameters. Mortality also occurs due to pesticide exposure and is modeled via a pesticide exposure module. Queen replacement can be simulated to represent a common beekeeping practice. Queen longevity. Colony size increase due to egg-laying rates, inverse for mortality.
Colony population structure and stage-specific rates [22,29,[39][40][41] Stage transition rates and proportion of drones influenced by photoperiod, local climate, colony size, and resource availability. Pesticide tolerance varies by stage and across colonies.
Stage-specific populations (egg, larvae, pupae, adults (drones, workers, foragers) controlled by egg-laying rate and stage-specific development rates and division of duties. Parameter specifying proportion of eggs that become workers and drones influenced by time of year (photoperiod). Parameter specifying proportion of eligible foragers actively foraging. Parameter specifying amount of sperm influences laying of fertilized eggs (lack of sperm results in laying drones).
Colony resource collection (by foragers) [22,29,[42][43][44] Empirical data to determine thresholds for foraging activity are limited. Forager activity influenced by temperature, wind velocity, and rainfall. Forage range varies from several hundred meters to 5500 m from colony. Foraging range and time influenced by landscape composition.
Foraging activity limited to threshold as determined by weather parameters (e.g., rainfall, temperature * A honey bee colony is considered a super organism. Colony growth is used to describe the change in number of individuals instead of the term reproduction. For honey bees, the term reproduction would be more accurate when describing mating between drones and queens from different colonies. Table 3. Parameters associated with pesticide exposure and effects.

Characteristic Ecological Concept Specification/Parameter (Simplifying Assumptions)
Exposure pathway [45] Consumption of nectar and pollen contaminated with pesticide(s). Contact exposure for foraging bees.
In colony, bees (larvae and adults) have age and stage-specific food consumption rates. Their pesticide dose is a function of food intake rate and pesticide concentration in pollen and nectar. Does not account for in-colony transfer and transformation of pesticide residues. Pollen and nectar (honey) stores are assumed to be well mixed, but they are actually stored in individual cells with variable pesticide concentrations. Does not account for exposure via plant guttation water.
Timing of exposures [46] Contact exposure occurs during application. Dietary exposure happens while crops are blooming.
Migratory bee colonies (for pollination services) may have repeated exposures. Exposure may occur past the growing season when stored nectar and pollen are contaminated with persistent pesticide.

Exposure pattern within habitat
Approach assumes that colonies are only feeding on treated crops.
Available information indicates that bees forage on a variety of plants, not just crops. Model can be adjusted to calculate percent of foraging on treated crop needed to impact a colony.
Exposure profile by stage [18] Stage-specific food consumption rates account for different exposure levels among life stages.
BeeREX includes information on stage and age-specific food consumption rates. Food consumption rates can be varied in Monte Carlo simulation.

Representation of toxic effects
Standard laboratory-based toxicity information for adult and larval bees used to calculate magnitude of mortality resulting from specific doses. Accounts for exposures and effects of a pesticide active ingredient.
Magnitude of mortality in a cohort corresponds to magnitude of exposure and dose-response curve from standard LD 50 study. Does not account for sublethal effects; however, user can adjust adult longevity or decrease foraging activity if toxicity data are available.
Does not incorporate queen mortality resulting from pesticide exposure, user can adjust queen longevity/replacement. Queen exposure assumed to be lower than workers because workers consume pollen and nectar and queen consumes jelly. Does not account for effects of multiple pesticide active ingredient exposures. Table 4. Other factors.
Biological stressors associated with Varroa mites, including viral infection, are accounted for directly by Varroa infestation in pupal cells and indirectly by manipulating adult worker longevity rates. Other stressors are not considered.
Beekeeper management practices [50] The model simulates impacts of miticide treatments on reducing Varroa impacts. Supplemental feeding can be considered.
Interactions between in-colony medications/miticides and pesticides.
Environmental stressors [22] Egg laying rates influenced by temperature and photoperiod Foraging influenced by temperature, wind speed, and rainfall Weather parameters chosen to represent geographic location of interest. Some days may be partial foraging days.
Landscape composition [29] Assemblage of natural and cultivated plants determine nectar and pollen resource availability, which influences foraging success and distance of foraging. Availability of pollen and nectar varies over time based on phenology of plants.
Landscape composition and corresponding resource availability is temporo-spatially determined. Model assumes that bees forage 100% on treated crop; however, user can evaluate impact of this assumption by calculating % foraging on crop needed to impact colonies. Effect of landscape on forager lifespan and mortality is not considered.
The third decision step addresses population and spatial factors. Core colony demographics such as density dependence, foraging, population dynamics, and habitat characteristics were inherited from the original model. However, a spatial objective for regulatory use was needed to account for geographic variation in use rates, crops, climate, and landscapes. Therefore, the ability to input weather data used in another USEPA model [44] was added to account for influences of temperature and rain on foraging. Routines that adjust colony behavior for this range of conditions were refined. The USEPA has hundreds of weather files available for locations across the US and the format is implementable for other locations. Currently though, a foraging area for a colony does not have spatial variability in the types and amount of pollen available to it during a growing season. As discussed above, it was assumed that the foraging bees collect 100% of their pollen and nectar from the simulated crop.
based on designing the model to have the ability to translate individual-level toxicity and exposure information to predict colony-level effects and outcomes and to inform the design and interpretation of semi-field and full-field colony-level empirical studies submitted for pesticide registration.
The first decision step (i.e., representation of the life cycle structure) is depicted in Figure 2. The model is resolved daily, with tracking of population sizes of eggs, larvae, and different adult classes. The second decision step (i.e., organism-level processes for reproduction (egg-laying), maturation, and mortality) is addressed through the original implementation of BEEPOP [22]. Therefore, the core population model has stressor-independent inputs for life history stage growth, promotion, resource collection, and mortality (Table 2). However, additional model inputs were needed for the revised model to have utility for pesticide evaluations. Additional classes of input parameters for risk assessment include pesticide physicochemical properties, toxicity information, method-specific application rates, and crop phenology information (Table 3).

Figure 2.
Conceptual model depicting BeePop+, including colony dynamics, Varroa mite infestation, and mortality due to pesticide exposure. Adapted from [24]. The flow chart at the top of the diagram represents the process included in VARROAPOP.
The third decision step addresses population and spatial factors. Core colony demographics such as density dependence, foraging, population dynamics, and habitat characteristics were inherited from the original model. However, a spatial objective for regulatory use was needed to account for geographic variation in use rates, crops, climate, and landscapes. Therefore, the ability to input weather data used in another USEPA Figure 2. Conceptual model depicting BeePop+, including colony dynamics, Varroa mite infestation, and mortality due to pesticide exposure. Adapted from [24]. The flow chart at the top of the diagram represents the process included in VARROAPOP.
The fourth decision step is for external factors that can be considered in the model. Pollen and nectar exposure considerations fall under this category. Exposure concentrations were added to the model and estimated using the approaches outlined in the BeeREX model [18]. Alternatively, empirical pollen and nectar residue data for pesticides can be used to estimate exposures for colonies. Exposures in jelly are assumed to be orders of magnitude below those of pollen and nectar [19].
The fifth decision step relates to exposure routes, effects, and temporal resolution. Exposure and effects algorithms were added, resolved at a daily time step, to account for the accumulation and distribution of pesticide throughout the hive and mortality effects within each of the life history stages. The Supplementary Materials includes detailed descriptions of these additions to the original VARROAPOP model [21] that allows for consideration of pesticide exposure and associated effects.

Phase 4: Conceptual Model
As depicted in Figure 2, the conceptual model of BeePop+ is centered around the simulation of colony dynamics by accounting for the different life history stages while evaluating the effects of stressors. Stressors accounted for in the model include Varroa mite infestation, weather conditions, and pesticide exposure. Impacts of beekeeper practices on the colony, such as Varroa mite control treatment, re-queening, and supplemental feeding can also be simulated. Details underlying the simulation of honey bee colony and Varroa mite dynamics are available in [22] and [21], respectively. Full details on the approach to estimate pesticide exposure and mortality of larvae and worker bees are provided in the Supplementary Materials.
The pesticide component of BeePop+ was developed to be consistent with Tier II of the risk assessment framework (Figure 1). The model was designed to utilize available information from Tier I (i.e., individual-level toxicity and exposure information). As with a colony-level study, exposure is meant to be controlled (e.g., 100% exposure to a treated crop) so that the assessor can assess whether individual-level effects may rise to a colony-level effect. The model was also designed to evaluate whether a proposed use on a pesticide product label (including application rate information for specific crops) may impact honey bee colonies. Details on the model equations and parameters are provided in the Supplementary Materials.

Phase 5: Model Implementation
BeePop+ is a collaborative product by the USDA and USEPA that emerged from the previous simulation models BEEPOP [22] and VARROAPOP [21]. BEEPOP simulates the number of individuals within a single colony, with consideration of different life stages (eggs, larvae, pupae, and adults), types of adults (e.g., in-colony, foragers), and condition of the queen. Major influences on the number of individuals within a colony include queen strength, longevity of adult workers, and weather. Varroa mites affect colony dynamics by reducing the longevity of adult workers parasitized during development [51,52]. To capture this, VARROAPOP was developed as an extension of the BEEPOP colony population dynamics model. VARROAPOP couples mite population dynamics with colony growth based on queen egg-laying rates and worker longevity [21,23]. An overview of the model routine is shown in Figure 2. The flow chart at the top of the schematic diagram ( Figure 2) represents the VARROAPOP process. Further detail on core colony model implementation of BEEPOP and VARROAPOP can be found in the original publications and the Supplementary Materials.
For BeePop+, pesticide exposure and effects components are integrated into the existing VARROAPOP model. This component is depicted in Figure 2, described briefly in Section 3.3, and presented in detail in the Supplementary Materials. The resulting BeePop+ model can simulate impacts of multiple stressors (e.g., Varroa mites, weather, pesticides) on honey bee colonies. The model is designed to simulate impacts of pesticide exposures resulting from applications to agricultural crops. Intended model users include risk assessors and scientists designing and conducting honey bee colony studies with pesticides. This model meets most of the desired features identified by USEPA. At this time, the implementation and evaluation phases are still on going and this model is not yet ready for regulatory use. Two major areas of development are needed to complete these phases: additional evaluations of the model with empirical data and development of the final model tool and associated user's guidance.
The algorithms of BeePop+ have been implemented in the C++ programming language and are available as source code [53] which can be compiled as a shared library. Additionally, a VarroaPy wrapper [54] is available for BeePop+. The wrapper allows BeePop+ to be run programmatically via Python, enabling high-throughput analyses such as Monte Carlo sampling. Model inputs are specified via Python dictionaries and outputs are returned as Pandas DataFrame or JSON data formats. Future areas of development include creation of a graphical user interface (GUI) and guidance to facilitate use of the tool by risk assessors. Additional evaluation and peer review may also be needed.
The BeePop+ model has been formally evaluated twice. First, a sensitivity analysis was conducted on the BeePop+ model [24] to identify the key parameters affecting colony trajectory across several pesticide exposure scenarios. This analysis found that queen strength (indicative of egg-laying rate) is the primary determinant of colony size in both exposed and unexposed scenarios. While initial colonies with queens of low strength rarely have viable colony sizes, interactions between queens of moderate and high strength with high average forager lifespans allowed for stronger colonies that were more resilient to pesticide exposure. In addition, the relative importance of exposure-and toxicity-related parameters depended on the pesticide application method: foliar applications, adult contact LD50, and application rate were relatively important parameters, and soil applications of systemic pesticides and the pollen load returned by foragers was more important.
Second, Bayesian inference was used to estimate key parameter distributions and compare BEEPOP+ outputs to the results of an empirical colony feeding study testing an insecticide (i.e., clothianidin) [25]. Results showed good agreement between predicted and observed effects on colony size, with empirical bee counts falling within the model's 95% prediction intervals for 21 out of 24 datapoints (Figure 3). The model was found to underestimate the negative effects at a low exposure level (36 µg/kg) and overestimate colony recovery from negative effects at the highest exposure level (140 µg/kg). However, these deviations may have been due to sublethal effects and interactions with other stressors which are not currently considered in BeePop+. From this study, an adult oral median lethal dose (LD50) of 0.0189 µg clothianidin/bee was estimated (Figure 4), which falls close to the range of values from laboratory toxicology studies on individual bees (0.0026-0.0157 µg/bee) ( [55,56]).
One avenue of future evaluations could be to use the Bayesian approach with other colony feeding studies for clothianidin to (1) evaluate the robustness of the model to other weather and climate scenarios and (2) discern the impact of natural stressors on colony response to the same chemical. In addition, it would be informative to evaluate the model's predictability with other pesticides, particularly those with different modes of action or life stages of susceptibility. For example, as clothianidin's primary effect is on adult worker survival, it would be informative to evaluate the model with a pesticide targeting larval survival. Additional colony level toxicity studies for other pesticides (including different modes of action than clothianidin) are available for these types of evaluations. Another area of future evaluation could also include investigating the interactive effects of pesticides and other stressors. One avenue of future evaluations could be to use the Bayesian approach with other colony feeding studies for clothianidin to (1) evaluate the robustness of the model to other weather and climate scenarios and (2) discern the impact of natural stressors on colony response to the same chemical. In addition, it would be informative to evaluate the model's predictability with other pesticides, particularly those with different modes of action or life stages of susceptibility. For example, as clothianidin's primary effect is on adult worker survival, it would be informative to evaluate the model with a pesticide targeting larval survival. Additional colony level toxicity studies for other pesticides (including different modes of action than clothianidin) are available for these types of evaluations. Another area of future evaluation could also include investigating the interactive effects of pesticides and other stressors. . Adult oral dose-response curve for clothianidin inferred from BeePop+ and the empirical feeding study data. Solid blue lines represent the median prediction and shaded regions denote the 68%, 95%, and 99% prediction intervals (PIs). Reproduced from [25]. . Adult oral dose-response curve for clothianidin inferred from BeePop+ and the empirical feeding study data. Solid blue lines represent the median prediction and shaded regions denote the 68%, 95%, and 99% prediction intervals (PIs). Reproduced from [25].

Conclusions
Colony simulation models can be useful tools for assessing risks of pesticides to honey bees. The current approach for evaluating risks of pesticides to honey bees relies heavily upon resource intensive empirical studies. The BeePop+ model was designed to help reduce the reliance on these empirical studies, providing a tool that can evaluate colonylevel impacts from effects on individual bees. Additionally, this model may be useful as a research tool to test stressor-related hypotheses and ultimately to help better design those empirical studies when needed. The original colony dynamics model upon which BeePop+ is based, BEEPOP [22], provides a core set of features that combine growth factors including egg laying, brood development with subsequent worker aging, and weather data (temperature, rainfall, sunlight). Multiple models have since been implemented (reviewed in [57,58]) that provide the ability to simulate stressors on colony-level outcomes. BeePop+ implements the original BEEPOP model and many features that are found in subsequent models, including in-colony age-structured dynamics, effects of parasites (e.g., Varroa mites), temporal foraging for pollen/nectar collection, and impacts of weather. In this study, we have documented our decision-making process of adding the capability to BeePop+ to simulate pesticide exposures from forager-collected pollen and nectar and how these exposures affect the age structure and survival of the colony. This model is still in the evaluation and implementation phase of development and is not yet ready for regulatory use. Future work could involve additional comparisons of model predictions and empirical colony-level toxicity studies involving other pesticides.
Supplementary Materials: Supplementary materials can be downloaded at: https://www.mdpi. com/article/10.3390/ecologies3030022/s1, The supplementary materials includes a detailed description of the additions to the VARROAPOP model that account for pesticide exposure and effects on individuals within a simulated honey bee colony. This includes equations, parameter definitions, and assumed values. The first document (BeePOP+_supporting_infoS1) containing supplemental information includes the details on the pesticide component of the model that estimates pesticides exposures and effects on bees within a simulated colony. The second document (BeePOP+_supporting_infoS2) describes the approach used to simulate colony failure due to low pollen and nectar stores during winter. References [59][60][61][62] are cited in Supplementary Materials File.