Social Buffering of Pesticides in Bumblebees: Agent-Based Modeling of the Effects of Colony Size and Neonicotinoid Exposure on Behavior Within Nests

Neonicotinoids are a globally prevalent class of pesticides that can negatively impact bees and the pollination services they provide. While there is evidence suggesting that colony size may play an important role in mitigating neonicotinoid exposure in bees, mechanisms underlying these effects are not well understood. Here, a recently developed, agent-based computational model is used to investigate how the effects of sub-lethal neonicotinoid exposure on intranest behavior of bumblebees (Bombus impatiens) are modulated by colony size. Simulations from the model, parameterized using empirical data on bumblebee workers exposed to imidacloprid (a common neonicotinoid pesticide), suggest that colony size has significant effects on behavior neonicotinoid-sensitivity within bumblebee nests. Specifically, differences are reduced between treated and untreated workers in larger colonies for several key aspects of behavior within nests. Our results suggest that changes in both number of workers and nest architecture may contribute to making larger colonies less sensitive to pesticide exposure.

Colony size is a factor in sensitivity to environmental stressors as it plays a role in the resilience and conservation of social insects (Straub et al., 2015;Fisher et al., 2019). Colony size may affect susceptibility to neonicotinoid pesticides. One study (Rundlöf et al., 2015) found no significant effects of neonicotinoid exposure on honeybees (large colony sizes, >10,000), intermediate effects on bumblebees (moderate colony size, < 500 workers), and very strong effects on solitary bees, under the same exposure conditions. Some negative effects of imidacloprid on honeybees were shown to be stronger in smaller colonies (Wu-Smart and Spivak, 2016).
The effects of colony size on worker behavior within nests in response to pesticide exposure are not well understood. Colony size may improve resilience to pesticides as a direct effect of larger population size by decreasing the chances of colony failure when a small number of workers are lost to disease or death or behavioral impairment (Straub et al., 2015). Alternatively, larger colonies could mitigate the effects of exposure on the behavior of individuals. Changes in interaction rate or nest architecture associated with larger colony size could modulate the effects of pesticide exposure on individuals and colony performance through the division of labor, behavioral development, or social interactions associated with large social insect colonies.
Our goal is to elucidate the mechanisms underlying the effects of neonicotinoid exposure on worker behavior within bumblebee nests and how these effects are modulated by colony size. We use our spatially-explicit agent-based model called BeeNestABM (Crall et al., 2018b;Ford Versypt et al., 2018) to test the hypothesis that the effects of exposure on worker behavior within bumblebee nests are reduced in larger colonies. Bumblebee colonies range from a single individual during colony founding by a solitary queen up to several hundred workers (Michener and Laberge, 1954). Inside colonies, interactions with nestmates and nest architecture can have important effects on worker behavior. Interactions with nestmates directly affect activity and task switching in colonies (Renner and Nieh, 2008;Crall et al., 2018b) and could modulate behavioral impacts of pesticides. The physical structure of the nest (e.g., food pots and developing brood) can modulate worker behavior by carrying information on larval development and hunger (Heinrich, 1974;Dornhaus and Chittka, 2005;den Boer and Duchateau, 2006). Because neonicotinoids directly reduce activity, which indirectly alters spatial and social dynamics (Crall et al., 2018b), increased interactions with nestmates or nest structures in larger colonies could mitigate the impacts of neonicotinoid exposure.

MODEL DESCRIPTION
BeeNestABM was parameterized using empirical data (Crall et al., 2018b) and simulates bumblebee activity and interactions within a nest. The code and documentation are available (Ford Versypt et al., 2018). Here, we briefly describe the model.
BeeNestABM simulates movements of bees within a nest on short time scales allowing for interactions with nestmates and nest structures. The state variables at each time step and for each bee are the x-and y-coordinates, heading angle, speed, and activity state. The simulation is initialized with a user-specified number of bees. Bees are initialized with random positions inside a nest arena (0.2 × 0.2 m for all simulations here). The model is iterated in steps of 0.5 s. Nest structures are placed in 1 cm apart in a square grid covering a centered, user-specified fraction of the arena.
The model consists of three processes that occur at each time step. First, distances between bees and the nest structures ( Figure 1A) are calculated. Next, the activity state for each bee (Figures 1C-G) is updated. Finally, the coordinates for each active bee are updated using the movement rules ( Figure 1B).

Activity State Updates
Activity-switching probabilities depend on location on or off the nest structures ( Figure 1D), proximity to nestmates ( Figure 1E), and treatment group (Figures 1C-G). Workers are assigned probabilities calculated previously (Crall et al., 2018b Supplement) from the means of empirical data of automatically tracked (Crall et al., 2015) bumblebees. These probabilities represent the probabilities of switching between active (moving) and inactive (stationary) between time steps, calculated for different combinations of social interaction state and location (Figure 1). Probabilities for treated bees are set equivalent to controls for parameters that did not differ significantly between control and treated workers (Figures 1F,G). Treated bees have activity-switching probabilities that make bees less likely to initiate movement either spontaneously or after physical contact with a nestmate and more likely to become immobilized, with stronger effects when treated bees are located off the nest structures (Figures 1F,G).
In the dataset (Crall et al., 2018b), workers were exposed orally to an acute dose of either 0 ng (control) or 1.0 ng imidacloprid at 87 ppb (treated) in a sucrose solution and allowed to behave freely in a nest arena. 1.0 ng was approximately equal to the cumulative imidacloprid consumed per worker in a single day of chronic feeding on nectar with environmentally realistic imidacloprid concentrations (6 ppb). While that study found qualitatively similar effects on behavior in workers exposed to imidacloprid either acutely or chronically, the behavioral effects during chronic exposure varied with time of day. To limit model complexity, we did not incorporate circadian effects  Figure  S4 and data from Figure S5 in Crall et al. (2018b).
but focused on acute exposure data collected at the same time of day.

Movement Rules
Previously inactive bees are set at a speed sampled from the empirical speed distribution. Bees previously moving are updated to the sampled speed if it is within specified limits of the current speed; otherwise, there is a 10% probability of switching to the sampled value from the current speed. The heading angle for each bee is the weighted angular mean of a random walk angle fluctuating within 90 • of the current angle and the angle toward the environmental stimuli determined by the net displacement between the bee and nest structures ( Figure 1B). The environmental stimuli represent attractions to the nest structures with weights estimated using empirical data (µ = 0.067 and 0.052 for control and treated bees, respectively). After location updates, movement is truncated for each bee that would move outside the perimeter.

RESULTS AND DISCUSSION
We used BeeNestABM to examine the effects of colony size and exposure treatment intensities (the proportion of workers within the colony exposed to imidacloprid) on worker behavior. We ran simulations where two aspects of colony size (numbers of nest structures and of workers) were varied simultaneously from 4 to 441, spanning the ecologically relevant colony sizes for bumblebees, to maintain a 1:1 worker:nest structure ratio across a range of treatment intensities (Figure 2). To examine the relative importance of nest architecture and number of workers, we ran simulations varying these parameters independently (Figures S1,  S2). For each combination of colony size and exposure intensity, we estimated mean values across all simulations and individuals for four behavioral metrics (Figures 2A-D and  Figures S1, S2). Exposure intensity and colony size significantly affected worker behavior within nests. Larger colonies had higher activity levels (Figure 2A), mean distance to the nest center (Figure 2B), and proportions of time spent on the nest structure ( Figure 2C), with interaction rates (Figure 2D) decreasing at intermediate colony sizes. Higher treatment intensity was associated with substantial declines in activity and interaction rate (Figures 2A,D), a weak increase in D NC (Figure 2B), and a weak decrease in P nest (Figure 2C). Simulations isolating the effects of variation in nest size and number of workers separately suggest that the effects of colony size are driven primarily by changes in the size and FIGURE 2 | Effects of colony size on emergent behavior in bumblebee nests. Mean behavioral metrics (left column) and relative differences between treated and control workers (right column) as a function of colony size and treatment intensity (% of workers treated). Data are shown for four metrics: (A,E) activity, (B,F) distance to the nest center (D NC ), (C,G) proportion of time spent on the nest structure (P nest ), and (D,H) interaction rate (R int ). Values are the mean across all workers and replicate simulations at a given parameter set. The number of simulation iterations at each parameter set was varied with colony size (# of iterations = 10000/N, where N = # of workers). Phase spaces smoothed using local averaging over adjacent values (smoothing window: ±10% for treatment intensity, (sqrt(N-2)) 2 :(sqrt(N+2)) 2 for N nest /N workers ). For all simulations, the number of nest structure elements (N nest ) and the number of workers (N workers ) were varied simultaneously to maintain a 1:1 nest structure: worker ratio. Differences between control and treated workers (right column) were calculated as a percentage of mean absolute value, abs([trtctrl]/mean(trt, ctrl)])*100, where trt denotes the treated workers and ctrl denotes the control group.
number of nest structures, especially for D NC and (trivially) P nest (Figure S1).
We found evidence that the relative difference between control and treated workers varied significantly with colony size (Figures 2E-H). Relative differences between control and treated workers in activity (Figure 2E), D NC (Figure 2F), and P nest (Figure 2G) were reduced in larger colonies, while relative differences in interaction rate were greatest at intermediate colony sizes (Figure 2H). The reduced differences between control and treated workers in larger colonies are likely driven primarily by the weaker effects of imidacloprid on activity parameters when located on the nest structures (Figures 1F,G and Figure S2). We found qualitatively similar effects for absolute differences between control and treated workers for most metrics (Figure S3), with the exception of proportion of time spent on the nest structure, which varied most with colony size (Figure 2 and Figure S3).

CONCLUSIONS
Our results support the hypothesis that larger colonies buffer the effects of neonicotinoid exposure on behavior within bumblebee nests: differences between control and treated workers in several key emergent behavioral parameters were reduced in larger colonies (Figures 2A-D). The dominant effects of nest structure size (Figures S1, S2) suggest that nest architecture, rather than altered rates of interaction with nestmates, is the primary driver of reduced effects of imidacloprid exposure in large colonies, although higher worker numbers did have qualitatively similar, but quantitatively weaker, effects for several metrics ( Figure S2). The strong effects of nest architecture found here result directly from the empirical observation that neonicotinoids have particularly strong effects on activity when workers are located on the nest structure (Crall et al., 2018b). The specific mechanisms driving this patterns are not clear but suggest that some aspects of the nest structure (potentially including chemical signals derived from the brood, den Boer and Duchateau, 2006 or from nestmates, Richardson and Gorochowski, 2015 laid on the nest structure) play a direct role in regulating worker behavior. Being on the nest structure stimulates movement (Figure 1), potentially mitigating the negative effects of neonicotinoids on activity.
The mitigating effects of colony size likely translate into significant impacts on colony performance. Location on the nest structure is a strong proxy for direct brood care and maintenance behavior within the nest (Crall et al., 2018a). Our results lead to the testable predictions that (a) similar levels of exposure will lead to stronger effects on nest behavior in smaller colonies and that (b) this would lead to stronger reductions in per-worker nursing rates and brood growth in smaller colonies. For bumblebees, increased sensitivity at small colony sizes is especially relevant because exposure during colony initiation by a solitary queen or during early developmental stages can substantially impair colony performance (Baron et al., 2017;Leza et al., 2018).
Our model improves understanding of how colony size and pesticide exposure affect behavior within bumblebee nests. However, the scope of the model is limited to a simplified set of behaviors. In the future, the model could be expanded to more complex behaviors of bumblebees (Cameron, 1989;Crall et al., 2018a). Combining this model with colony dynamic processes at larger spatial and longer time scales, including foraging, growth, and landscape dynamics (Becher et al., 2018) could help assess how the impacts of neonicotinoids on behavior translate to long-term colony performance.

AUTHOR CONTRIBUTIONS
JC, BD, and AF contributed equally to this work. JC, BD, and AF designed the study and wrote the computational code. JC shared relevant experimental data from a previous study and conducted the in silico experiments. JC, BdB, BD, and AF analyzed results of the in silico experiments and wrote the manuscript.

ACKNOWLEDGMENTS
The authors acknowledge support from the Statistical and Applied Mathematical Sciences Institute (SAMSI) and the organizers of the workshop titled Data-Driven Modeling of Collective Behavior and Emergent Phenomena in Biology that enabled the collaboration. JC was supported in part by funding from the Rockefeller Foundation and the Winslow Foundation. BdB was supported by the NSF (IOS-1557913), the Alfred P. Sloan Foundation, The Klingenstein-Simons Fellowship, and the Smith Family Foundation.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00051/full#supplementary-material Figure S1 | Independent effects of nest structure size and worker number on emergent behavior in bumblebee nests. Mean behavioral metrics when the number of nest structural elements (N nest, left column) and number of workers (N workers , right column) are independently varied. Data are shown for four metrics: (A,E) activity, (B,F) distance to the nest center (D NC ), (C,G) proportion of time spent on the nest structure (P nest ), and (D,H) interaction rate (R int ). Values are the mean across all workers and replicate simulations at a given parameter set. The number of simulation iterations at each parameter set was varied with colony size (# of iterations = 10000/N, where N = # of workers). Phase spaces smoothed using local averaging over adjacent values (smoothing window: ±10% for treatment intensity, (sqrt(N-2)) 2 : (sqrt(N+2)) 2 for N nest /N workers ). For each metric, identical color scales are used for the left and right columns. Figure S2 | Independent effects of nest structure size and worker number on relative effect size of imidacloprid. Relative difference between treated and untreated workers when the number of nest structures (N nest, left column) and number of workers (N workers , right column) are independently varied. Data are shown for four metrics: (A,E) activity, (B,F) distance to the nest center (D NC ), (C,G) proportion of time spent on the nest structure (P nest ), and (D,H) interaction rate (R int ). Differences calculated as in Figure 2. The number of simulation iterations at each parameter set was varied with colony size (# of iterations = 10000/N, where N = # of workers). Phase spaces smoothed using local averaging over adjacent values (smoothing window: ±10% for treatment intensity, (sqrt(N-2)) 2 : (sqrt(N+2)) 2 for N nest /N workers ). For each metric, identical color scales are used for the left and right columns. Figure S3 | Absolute differences between treated and control workers across treatment intensity and colony size. Simulations are the same as in Figure 2. Relative differences for (A) activity, (B) distance to the nest center (D NC ), (C) proportion of time spent on the nest structure (P nest ), and (D) interaction rate (R int ) are shown as the absolute difference (rather than % difference in the right column of Figure S2).
Frontiers in Ecology and Evolution | www.frontiersin.org