A model of wild bee populations accounting for spatial heterogeneity and climate‐induced temporal variability of food resources at the landscape level

Abstract The viability of wild bee populations and the pollination services that they provide are driven by the availability of food resources during their activity period and within the surroundings of their nesting sites. Changes in climate and land use influence the availability of these resources and are major threats to declining bee populations. Because wild bees may be vulnerable to interactions between these threats, spatially explicit models of population dynamics that capture how bee populations jointly respond to land use at a landscape scale and weather are needed. Here, we developed a spatially and temporally explicit theoretical model of wild bee populations aiming for a middle ground between the existing mapping of visitation rates using foraging equations and more refined agent‐based modeling. The model is developed for Bombus sp. and captures within‐season colony dynamics. The model describes mechanistically foraging at the colony level and temporal population dynamics for an average colony at the landscape level. Stages in population dynamics are temperature‐dependent triggered by a theoretical generalized seasonal progression, which can be informed by growing degree days. The purpose of the LandscapePhenoBee model is to evaluate the impact of system changes and within‐season variability in resources on bee population sizes and crop visitation rates. In a simulation study, we used the model to evaluate the impact of the shortage of food resources in the landscape arising from extreme drought events in different types of landscapes (ranging from different proportions of semi‐natural habitats and early and late flowering crops) on bumblebee populations.


| INTRODUC TI ON
Bees are key pollinators in many natural ecosystems and provide pollination services for agricultural crops (Klein et al., 2007). However, domestic and wild bees are threatened by consequences of climate change and land-use change (IPBES, 2016;Soroye et al., 2020;Vanbergen & the Pollinators Initiative, 2013) and have been negatively affected in many parts of the world, at the local and regional scale (Biesmeijer et al., 2006;Potts et al., 2010). A good quantitative understanding of global change effects on pollinators is important to understand the consequences for pollination services and the need for conservation. Pollinator models are increasingly used in ecosystem service mapping and impact assessments to relate land-use patterns to the population size of pollinators such as bees (Becher et al., 2018;Gardner et al., 2021;Koh et al., 2016). The consideration of variability in weather and climate has received less attention, despite its importance in driving demographic rates, population sizes (Selwood et al., 2015), and distributions (Aguirre-Gutiérrez et al., 2017). Indeed, land use and climate change have interactive effects on pollinators (Oliver et al., 2015;Prestele et al., 2021) and thus should be considered jointly.
Climate change effects include gradual changes in annual means, altered seasonal variation, and increased frequency of extreme events, such as heatwaves and extended drought periods (IPCC, 2021). The phenologies of many ecological processes are moderated by temperature and thus sensitive to climate change.
Bumblebees (genus Bombus) contribute critically to the provision of crop pollination (Table S2 from Kleijn et al., 2015). The presence of floral (i.e., semi-natural habitats (SNH) and flowering crops) and nesting resources (i.e., natural and SNH, field edges, forests, meadows, and permanent grasslands) in the landscape during their life span is key for their colony development. This includes among others the ability to produce large foragers (Persson & Smith, 2011), to produce males and new queens at the end of the season (Rundlöf et al., 2014), and to survive during winter (Persson & Smith, 2013). Intensive agricultural management in large arable fields reduces the availability of nesting sites, with reduced crop diversity (Aizen et al., 2019) being associated with the dominance of individual flowering crops which may cause bottlenecks in terms of foraging resources for pollinators outside of the flowering period of these crops. Understanding how spatial and temporal variability of resources driven by land use and climate change interact at the landscape level and affect pollinator populations is crucial to help ensure that crop demands and pollinator supplies are well-matched (Settele et al., 2016). Bees are affected by climate change and there are adaptive limits of this pollinator group to track climate change (Kerr et al., 2015;Potts et al., 2010;Prestele et al., 2021;Sirois-Delisle & Kerr, 2018). While the availability of habitats rich in resources for bees, such as SNH, may be able to offset the deleterious effects of climate change on bee communities (Papanikolaou et al., 2017), the interactions between the effects of land use and climate are still poorly studied.
Spatially explicit models of pollinators produce bee visitation rates from proxies of bee abundance and floral resources at a landscape scale, which is used as a representation for the supply of pollination (e.g., InVEST pollination module using Lonsdorf et al., 2009).
Visitation rates can be derived from central place foraging theory (Lonsdorf et al., 2009;Olsson et al., 2015), assuming that fitness is entirely dependent on the distribution of floral resources around the nest as derived from a spatially explicit land use map. These models produce indices of bee visitation rates for fixed floral resources dividing the flying season into two (Häussler et al., 2017) or three periods (see Gardner et al. (2020)). This allows for a limited variability in resources or visitation rates within and between seasons, and therefore makes it difficult to study more fine-scale variability in floral resources driven by climate (e.g., the different start of flowering).
Interactions between pollinators and land use, and thereby changes in nesting and floral resources, require models combining foraging theory with population dynamics. For these models to be sensitive to variability in resources in the landscape, there is a need to model interactions at relevant spatial and temporal scales and include changes in growing conditions due to climate conditions (Johansson & Bolmgren, 2019).
To fill this gap, we developed a spatially and temporally explicit model for wild bees, referred to from now on as the LandscapePhenoBee model, that uses the foraging function of an existing pollination model developed by Häussler et al. (2017).
LandscapePhenoBee can be used to explain and project population dynamics based on changes in colony dynamics driven by landscape components and weather-induced variability in resources. The landscape components consider variation in landscape-scale land use covers both differences in composition (proportions of different habitat types) and configuration (e.g., field sizes). Weather-induced variability of resources entails that the phenological growth of plant resources is triggered by a generalized seasonal progression, but also that extreme weather events (e.g., droughts), can influence the growth of resources. The different growth development of different stages at the colony level is also triggered by a theoretical generalized seasonal progression.
The aim of this study is to use the LandscapePhenoBee model to explore the effect of the temporary, drought-induced shortage in food resources on the population viability of bumblebees, evaluated by (1) population size and (2) production of queens, and on (3) the pollination services provided by bumblebees in different types of landscapes, ranging from simple (landscapes with a low proportion of SNH) to complex (landscapes with a high proportion of SNH) agricultural landscapes, and including early and late flowering crops. We expect that (1) bumblebee populations to be more severely affected by drought in less complex landscapes since the distance between a nest and floral resources is on average larger in a simple compared to complex landscapes and (2) the presence of early flowering crops will have an effect on the colony dynamics, expecting higher production of bumblebee workers with landscapes with lots of early flowering crops, while landscapes with more late-flowering crops will have a positive effect on the production of queens. The influence of the population parameters is evaluated by sensitivity analysis.

| Theoretical model description
The LandscapePhenoBee model is designed to simulate wild bee species that are central place foragers, for example, a bumblebee species. We considered a fictive common early active bumblebee that stands for several species that are known to be important crop visiting species, including Bombus terrestris, Bombus lucorum, and Bombus lapidarius (Kleijn et al., 2015). The foraging module builds on an existing bee foraging model developed by Häussler et al. (2017) but expresses visitation rates with a higher time resolution (per week, instead of two time periods) and focuses on the colony and population dynamics at the landscape scale within the year (instead of as in Häussler et al. between years) (see Table S1 with a more extended comparison between the two models). The population dynamics are expressed for an average population in the landscape, which is informed by the total amount of resources gathered by nests in the landscape. The model produces the development of the population size and pollination services in the total landscape as outputs. Input and output model parameters are described in the Tables S1 and S2.
2.1.1 | Spatially explicit resources in the landscape Spatial variability in nesting and floral resources is represented by spatially explicit maps with a resolution of 10 × 10 m. The size of a landscape (2010 × 2010 m) was considered to be sufficiently large to account for bee foraging distance, which is normally about 1000 m or F I G U R E 1 Representation of the spatial-temporal explicit model developed in this study including the climate and landscape inputs and the outputs for the wild bee population dynamics and the ecosystem service. The simulated landscapes include four different types of patches: SNH, early flowering crop, late-flowering crop, and non-bee habitat (habitat that does not provide any resources for bees) less for bumblebees (Osborne et al., 2008). The landscape consists of three types of habitat patches: early-and late-flowering crops, and SNH. The rest of the grid consists of a matrix of non-suitable habitats for the bees (representing for example non-flowering arable land, or sealed urban areas), from this point forward referred to as non-bee habitat. The resolution, landscape size, and categorization of habitats can be altered.
Semi-natural habitats are the only habitats (among those represented) where nesting is considered possible. The nests are assigned randomly in the landscape, and the density of nests in the landscape is determined by the nest density parameter. The nest allocation is fixed for a given season. We specified that each SNH cell can have either no nests or one nest (but this can be changed if the model runs in a different resolution). The total number of nests in a landscape is proportional to the amount of SNH in the landscape (the more SNH in the landscape, the more nests).
We created artificial landscapes (e.g., Figure 1) with different proportions and configurations of four land-use classes (SNH, early and late flowering crop, and non-bee habitat) using the R package landscapeR (Masante, 2016;Thomas et al., 2020). The number and size of the patches of SNH were used to randomly create landscapes with different proportions of SNH between 5% and 25% which was used to describe spatial heterogeneity between simple and complex landscapes (Holzschuh et al., 2016;Scheper et al., 2014). The average size of crop patches was set to 500 m 2 , and the proportion of flowering crops was set such that these constitute together with SNH, 60% of the area in the landscape.
Within the proportion of crops, fields of early-and late-flowering crops (MFC) were randomly assigned to achieve a certain proportion of early versus late. Consequently, the presence of these two crops was negatively correlated (the more early-flowering crop, the less late-flowering crop).
Temporal variability in floral resources is considered by mapping floral resources with a weekly temporal resolution based on the floral phenology model (see Section 2.1.2). The floral phenology model provides the start, peak, and end of flowering in each habitat type based on simulated growing degree days (GDD) such that for each week t, cell i has the floral resource value F(t,i) ( Figure 1, and see Figure S1).

| Floral phenological model
The influence of climate on population dynamics and temporal variability in the availability of floral resources for different land resources during the year is considered by linking the LandscapePhenoBee model to a theoretical generalized seasonal progression from 0 to 100, mimicking a sigmoid function of cumulative GDD for northern hemisphere context (see Figure S4). If daily observed temperatures are available, the model allows for a simple calculation of GDD and uses it as model input. In this manuscript, we present a seasonal progression with an arbitrary GDD.
We modeled the temporal dynamics in floral resources for each land use separately as a function of the generalized seasonal progression following a sigmoid function product of a cumulative standard normal distribution and defined for day y at time t as: where Φ is the probability function for the normal distribution ( Figure S4). For each land use type h, the start and end of floral resources are given by a theoretical gdd start,h and gdd end,h The floral resources in cell i at time t is where z h (t) is a standardized value between 0 and 1 of the day of the year corresponding to time t, derived by the following expression The maximum floral resources in habitat h, f max,h , is a theoretical parameter that has assigned a value between 0 and 1 to capture how the floral resources in different habitats relate to each other, chosen to have the following relations: early MFC > late MFC > SNH. The reason was to simulate the development of resources that these different habitats provide along the season, a peak of resources early in the season by early MFC, a lower peak of resources later in the season by late MFC, and lower but constant floral resources provided by SNH (see also Figure 3a).

| Bee foraging
The number of foraging bees from the nest in cell i at time t is X(t,i).
The foraging bees are initially overwintering queens, that is, X(t,i) is 1, and later workers are produced by the queens. The foraging function is an exponential kernel reweighted by the floral resources (Häussler et al., 2017). The rate at which each cell j is visited by foraging bees from cell i during week t is: where F(t,j) is the floral value of cell j, d(i,j) is the Euclidean distance between cells i and j, γ is the mean dispersal distance when foraging; U i is the set of cells reachable from cell i. The denominator in Equation 4 weighs the attractiveness of cell i compared to the total attractiveness of the cells in the landscape and by foraging distances. In this way, a cell further away from cell j compared to cell i but with higher floral resources compared to i can be receiving more visits. Thus, the resources collected correspond to the distance-weighted resource values from cells in which bees are nesting. The spatial layer is treated as a toroid, that is, the edges are connected to each other such that a bee moving beyond the boundary of the spatial layer appears at the opposing edge.
Under the assumption that there is no depletion of floral resources in the landscape, the resources collected by foragers emerging from the nest in cell i at time t is where U i are the cells within reach from cell i.
The average amount of resources gathered per nest in the landscape at time t is

| Bee colony and population dynamics
The model covers the active period of bees in the season for a given year, that is, from the emergence of queens in the spring until the production of daughter queens at the end of the summer.
Bee population dynamics are modeled as two main stages within a season, corresponding to who is foraging: either overwintering queens (1) or workers (2) (see Figure 2). In turn, these two stages are subdivided into stages depending on what is being reproduced: Stage A1: Overwintered queens are foraging, and workers are not yet being produced.
Stage A2: Overwintered queens are foraging, and workers are produced.
Stage B1: Workers are foraging and being produced.
Stage B2: Workers are foraging, and daughter queens are being produced.
Stage B3: No new individuals are produced, and the mortality of workers and daughter queens increases.
Queens emerge when theoretical GDD reaches φ (Table S2, Figure S4) and start producing workers after 2 weeks (from A1 to A2). The number of foragers during stage A1 and A2 are one per nest, that is, X(t,i) = 1 if there is a nest in cell i, and 0 otherwise.
In stage A2, the average number of workers per nest at time t, W(t), is given by: where δ t < 1 is the survival of the workers from time t−1 and w(t) is the average number of workers produced per nest at time t. Growth depends on the resources gathered during the two previous time periods (t−1 and t−2) according to a plateau function: where α t is the maximum number of workers that can be produced per nest at time t and β t is a parameter determining for which amount of resources half of the potential number of workers are being produced at time t. Survival rate and the two growth parameters are constant values δ, α, and β, respectively, during stages A2-B2.
Workers take over the foraging, that is, there is a transition from stages A2 to B1, when a theoretical GDD has reached μ (where μ > φ)  Figure S4). During this period the number of foragers The colonies begin production of daughter queens, that is, the transition from B1 to B2, when theoretical GDD reaches ψ (where ψ > μ) (Table S2, Figure S4)

| Drought
In this simulation study, we introduced a drought event as a reduction in floral resources that was defined by a fixed starting point early in the season and had a duration between 1 and 4 weeks (see Figure 3, Figures S2 and S3). The reduction of floral resources was simulated in a way that penalizes the growth of resources: if in normal conditions the growth is positive, the growth in drought conditions will be close to zero, while if the growth is negative, under drought conditions this will translate in a 50% larger reduction of growth. Drought and no-drought conditions were evaluated as a case-control setup, which allowed us to study the effect of drought by comparing the model results with and without drought while the rest of the model design, including the generated landscape and distribution of floral resources in space, was kept constant within case-control pairs. The start of the drought was defined by an arbitrary GDD above 5°C, which translates into a corresponding day of the year, following the other events triggered by arbitrary GDD (see Figure S4).

| Simulation design
To study the effect of landscape heterogeneity on the impact of   included the landscape identification number as a random factor (to be able to evaluate the case-control setup with the drought). We assessed the significance of the main effects using likelihood-ratio tests comparing models with and without the effect (Table S5).

| Sensitivity analysis
The The parameters fmax,h were not used in the sensitivity analysis because these characterize the input to the population model and are therefore not parameters of the population model as such, which was the target of the sensitivity analysis in the study. The nesting density parameter was included in the sensitivity parameter because it is used to derive the initial population size and can therefore be seen as a parameter of the population model.

| Regression analysis
All three output quantities were significantly affected by the pro-

portions of SNH and MFC and if there is a drought or not (p-values
for all LR statistics were <.01, see Table S8

| Sensitivity analysis
The sensitivity analysis showed that the parameter survival rate (δ) had the highest influence on the three output quantities ( Figure S5) as well as estimated effects from %SNH and % Early flowering crop and Drought ( Figure S6). The parameter, temperature-dependent timing for the start of queen production (ψ), and the two growth rate parameters (α and β) also had a high influence on the model outcome.
The nesting density and proportion of workers foraging (pw) had a minor influence. The parameters, arbitrary GDD, for workers to start foraging (μ), the proportion of workers produced compared to new queens (ε), and arbitrary GDD for queen emergence (φ), had the least influence on the estimates of the effects of SNH and early flowering crop and drought on all three output quantities. See Tables S4-S7 and Figures S5-S9 for the results of the sensitivity analysis).

| DISCUSS ION
The

| Sensitivity analysis
From the sensitivity analysis, we found that population dynamic parameters related to survival, growth, and the time for the start of queen production were the most influential. It is not surprising that

| SNH and nesting opportunities
SNH provide both nesting habitats and a continuous provision of food resources throughout the whole season. As per the design of our study, the number of nests increases with the amount of SNH.
Therefore, the SNH determines the initial colony conditions for bumblebee population growth in the landscape. After that initial stage, flower resources will influence population growth. Estimating the location and the number of bumblebee nests in a landscape remains a challenge in pollination ecology, since nests cannot be easily surveyed (Knight et al., 2005). We did not simulate a specific bumblebee species, but a fictive common bumblebee that stands for  has been observed in assessing drought and short-term temperature increases in taxa of flower-visiting insects, since Oliver et al. (2013) showed that quantity and low degree of fragmentation of habitats reduce sensitivity to drought in butterflies. Oliver et al. (2013Oliver et al. ( , 2015 also assessed recovery from drought but included many other factors involved, such as inter-patch dynamics between generations, which are not directly comparable to our results. Papanikolaou et al. (2017) showed, using long-term monitoring data from Germany, that a higher amount of SNH can mitigate the negative effects of shortterm increases in temperature on wild bee species richness. It is unclear whether the patterns found were due to the disproportionate importance of food resources in SNH under drought conditions, or whether these habitats provide cooler microhabitats. We did not consider differences in microclimates between habitats, which may be relevant in explaining interactive effects on bees between landuse and climate extremes (as discussed in Papanikolaou et al., 2017).

| Effects of SNH and impacts of drought
Additionally, drought periods were implemented as a reduction in all floral resources in the landscape, and further model developments could consider drought resistance between different land-use types, allowing for example to explore different drought-resistant crop varieties.

| Effects of MFC and impacts of drought
Our simulation results show that the impact of short periods of drought early in the season is not severe on population size but has negative impacts on both queen production and pollination potential. This can be explained as short drought periods have an intermediate decrease in food resources ( Figures S2 and S3) and bees are able to compensate for the lack of food resources in space and time (maintaining the population size, i.e., number of workers). However, although population size can still reach a high number under drought conditions, the colonies may have missed the peak of early flowering or the peak is reduced, translating to a loss of resources available to produce queens.
The production of queens was higher in landscapes where there was a high enough proportion of late-flowering crops. Where the ratio of early-flowering crops to late-flowering crops was too high, the production of queens decreased (Figure 4e).
In our model simulations, the more SNH, the more nests in the landscape, and the higher chance that a nest is close to floral resources in crops, which can explain the positive effect on pollination potential results. Pollination services from bees have been previously measured as visits per unit area (Lonsdorf et al., 2009), visits per flower (Rader et al., 2012), improved crop yield  or assumed as a direct proxy of species richness (Perennes et al., 2021).

| The use of temperature sums in the model
GDD is the accumulation of temperature above a certain base temperature for each calendar day, making it a good indicator to account for both spatial and temporal variation in temperature, including the lower developmental threshold in which plant growth development and flowering are possible. Since bees are sensitive to temperature change (Martinet et al., 2020;Pawlikowski et al., 2020), GDD has a great potential to predict insect phenology, that is, the development of insects (Cayton et al., 2015).
In this study, we have used a theoretical generalized seasonal progression, to represent a theoretical GDD value to mark the start and end of flowering for the flowering period, spring bee emergence, and worker foraging. A choice of the model was to also use an arbitrary GDD to trigger the initiation of the production of daughter queens. While the emergence of spring queens is triggered by temperature (Alford, 1969;Goodwin, 1995), the production of new queens is a complex combination of factors involving resources in the landscape, temperature, and health of the colony (Goulson, 2003). There are several approaches to model the switch point including time (Crone & Williams, 2016), and assessing the daily ratio of larvae-worker below a certain threshold (e.g., 3 used in BumblebeeHAVE (Becher et al., 2018)). By adding GDD in our model as a switching point, we could add variability to the time of queen production, by shifting the day of production, instead of a fixed day in time as in Crone and Williams (2016). We acknowledge that using temperature as a unique switching condition to produce new queens is not perfect, and some studies do not find temperature as the main reason for switching point (Holland & Bourke, 2015;Vogt, 1986).
However, there is a lack of experimental studies that control for different temperatures, or field data and modeling to construct a better understanding of how bumblebees or other social insect pollinators respond to changes in temperatures (but see Zaragoza-Trello et al., 2021). Additional carefully controlled experimental studies, for example with variable temperature regimes, in combination with data from the field and modeling, should help construct a fuller understanding of how major social insect pollinators are likely to respond to climate change.
An additional aspect worth mentioning is that we did not consider in our model the production of males, and the complex colony dynamics of queen-worker conflict derived from the switching point of producing reproductives in the colony (Goulson, 2003). Given the nature of the model to produce a representation of the bumblebee colony cycle pattern, the number of daughter queens was slightly higher in landscapes that had more late-flowering crops, as shown in empirical studies (Rundlöf et al., 2014). Results from the simulation indicate that landscapes with more late-flowering crops would increase the number of workers, and therefore bring more resources to the colony and produce more queens.

| CON CLUS IONS
Climate and land-use changes are two drivers of bee population decline that should be considered in combination. The LandscapePhenoBee is a new mechanistic pollination model that can account for landscape heterogeneity and temporal variability in resources to bees yet keeping the model relatively simple with a few parameters. By introducing climate-induced temporal variability of food resources, this model contributes to the methodology of studying and predicting the impact of important drivers and extreme events on wild bees. As an example, the theoretical model shows the ability to qualitatively reproduce patterns observed in the very limited number of studies combining landscape-scale availability of resources and drought (Oliver et al., 2013(Oliver et al., , 2015Papanikolaou et al., 2017).

ACK N OWLED G M ENTS
We thank the handling editor, subject editor, and two anonymous reviewers for their constructive and helpful comments that greatly improved the manuscript. The authors are affiliated with the strategic research environment BECC (Biodiversity and Ecosystem Services in a Changing Climate).

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The code and simulated data that support the findings of this study are openly available in the GitHub repository link: https://