The Unique ergodic stationary distribution of two stochastic SEIVS epidemic models with higher order perturbation

1.   Introduction

Honeybees, Apis mellifera, are essential for maintaining biodiversity in undisturbed ecosystems, and for the production of economically valuable agricultural crops worldwide [1]. One third of the food that we consume relies on pollination by bees for its production. Pollination by honeybees is valued between $\$$15 and $\$$20 billion annually in the US [2]. If higher-order industrial sectors that rely on the production of honey bee pollinated crops are included, the value of honey bees is an additional US$\$$10.3 to $\$$21.1 billion [3]. The work of Dr. vanEngelsdorp and Dr. Meixner [4] suggests that 1) although the global number of bees has been increased in the past 50 years, the geographical differences are huge even within the same region. And 2) the number of colonies in Europe and North America are both decreased. For the past 10 years, beekeepers have lost between 30–45% of their colonies annually with the most recent losses in 2019–2020 exceeding 40% (https://beeinformed.org/wp-content/uploads/2020/06/BIP_2019_2020_Losses_Abstract.pdf). Although there are more and more lost colonies, the reason why the number has not actually been decreased is that the beekeepers rapidly replace the lost colonies [4]. And beekeepers may adopt cold storage to help the colony survive over the winter smoothly [5]. These treatment methods are accompanied by relatively large costs in the context of apiculture. Therefore, the loss of the colony is a global economic, agricultural, and ecological crisis.

There is no single factor responsible for colony losses. Colonies die at different times of year due to a combination of stressors such as pesticides, lack of food sources and malnutrition, parasites such as Varroa destructor, Nosema ceranae, and viruses/pathogens including bacterial diseases. Varroa mites are considered to be one of the primary causes of colony losses worldwide [6-8]. Varroa can cause colonies to perish because parasitism compromises the bee's immune system and reduces the longevity of adult bees parasitized during development. Varroa also transmits and activates several viruses during feeding including Acute Bee Paralysis Virus (ABPV), Deformed Wing Virus (DWV) and Israeli Acute Paralysis Virus (IAPV). These viruses differ in their symptoms and virulence, but all reduce brood and adult bee survival and suppress colony growth. Colonies that are heavily infested eventually collapse from the combination of parasitism and disease especially overwintering in the U. S., Canada and Europe [4,7,9]. To prevent colony losses, beekeepers currently are applying miticides 5–7 times each year. Even with these measures, overwintering colony losses can exceed 30% [10]. Colonies also can be lost from Nosema infections. Nosema is a microsporidium and because of reduced metabolic capacities, relies heavily on its host to furnish energy for growth and reproduction [11,12]. Nosema infects adult bees and damages the midgut causing reduced nutrient absorption. Infected adult bees show symptoms of nutritional stress such as reduced larval feeding, earlier on set of foraging and reduced longevity of adult bees [13,14].

In addition to stress from pathogens and poor nutrition, colonies also can be lost from the impacts of pesticides. Pesticides can be toxic to honey bees due to their effects on the nervous system. However, pesticides also can have sublethal effects that exert stress that eventually lead to colony losses. For example, ingredients in pesticides can suppress the immune system, thus opening the way for parasitic infections and viral diseases [15-18].

The first mathematical model that incorporated colony population dynamics driven by biological and environmental factors was the BEEPOP model developed by DeGrandi-Hoffman et al. (1989) [19]. Other honeybee models followed (see review paper [20]) and can be classified as colony 'models' addressing within-hive colony dynamics (e.g., see [19,21-24]), 'Varroa models' exploring interaction between honeybees and Varroa mites (e.g., see [25-29]), also were developed as were 'foraging models' simulating the dynamics of honeybees collecting nectar and pollen to sustain colony growth (e.g., see [30-34]).

There is an increasing recognition of the importance of developing mathematically tractable models to obtain useful insights on ecological processes and factors that lead to colony losses. In this paper, we provide a review of recent progress in mathematical modeling of honeybee colony population dynamics since the review by [20]. We summarize existing models, and discuss their potential to identify combinations of factors leading to colony losses. Our review will focus on 1) honeybee only models (Section 2) with age and/or food/nutrient components; 2) honeybee-parasites and/or virus interactions models (Section 3); 3) models with effects from pesticides (Section 4); and 4) computational models (Section 5). We end by suggesting future directions for mathematical model development and applications. For convenience, in the Appendix, we also provide the equations of all models and tables that provide features of mathematical models that we review.

2.   Honeybee population systems

Honeybees are social insects that live in colonies consisting of a single egg laying queen, zero to several thousand reproductive males, and tens of thousands of reproductively sterile female workers, and 10,000–30,000 eggs, larvae and pupae (brood) [35-37]. Figure 1 provides a typical life cycle of a worker honeybee. Colonies rear brood and expand their populations in the spring and summer. Brood rearing declines in the fall and colony populations decrease in preparation for overwintering [38]. The population of workers in a colony is divided into two major temporal/behavioral subcastes: hive bees that remain in the colony and foragers that collect resources such as nectar and pollen from flowering plants, water to cool the hive, and plant resins (propolis) that has antibiotic properties and serves as part of the colony-level immune system) [39]. In general, tasks that worker bees perform are based on their age. The youngest adult workers (0–3 days old) clean cells in the brood nest, and usually by day-three begin feeding larvae. When workers are 10–14 days old, they participate in nest construction, food processing, and guarding the entrance of the nest. In the third or fourth week of life, workers transition to foraging outside the nest.

Though the dynamics of a colony population are complex, models that incorporate simple rules can reveal biologically relevant mechanics that drive population growth. For example, Dennis and Kemp (2016) [23] proposed and studied a deterministic, one-state-variable model that formulates the recruitment and loss rates of colony members in terms of density dependent social components as follows:

where $A$ is the number of adult bees, $\lambda$ is the upper maximum as adult population grows, and $\theta$ is the number of bees at which new worker population is half of maximum. The mortality function $(\frac{\alpha\beta}{\beta+A}+\mu)$ sets a minimum of per-individual mortality rate ($\mu$) and a maximum of per-individual mortality rate ($\mu+\alpha$) for small adult population values; and $\beta$ represents social effects on bees mortality. The model (2.1) incorporates cooperative effects in both birth and mortality terms that lead to strong Allee effects where there exists a critical size in adult bee numbers, below which mortality is greater than recruitment, with ensuing colony loss. The model suggests that the strong Allee effect may be the ultimate cause of colony death.

Kang et al. (2016) [40] developed a simple model that includes the cooperative effects in the honeybee colony as follows:

where $H$ is the honeybee population, $r$ indicates daily egg-laying rate of queen, $\sqrt{K}$ indicates colony size where brood survival rate is half maximum, $d$ indicates mortality rate of bees. The model (2.2) used a Holling-Type III functional response ($\frac{H^2}{K+H^2}$) to model collaborative efforts of hive and foraging bees to guarantee the survival of newborns. This approach is motivated by the work of Kang et al. (2011) [41] and Real (1977) [42]. We would like to point out that Dr. Eberl in his 2010 work [43] introduced a sigmoidal Hill function ($\frac{H^n}{K^n+H^n}$) to model eclosion. Both work of Kang et al. (2016) [40] and Eberl et al. (2010) [43] have been applied in new research on honeybee population dynamics (e.g, see [22,29,44-48,50]). The nonlinear term $\frac{H^2}{K+H^2}$ in model (2.2) can generate a saddle-node bifurcation in the system such that it induces an Allee threshold above which the colony can survive and below which the colony collapses. The model (2.2) has been extended to include age-structure [22], and the interactions with mites [29,48,49] and viruses [40].

The models (2.1) and (2.2) focus on colony populations without age-structure [23,40]. Honeybees go through different stages, thus it seems to be more realistic to propose and study population models with age-structure. In the next subsection, we will review honeybee models with age-structure, and the related models with the food resources and/or nutrition factors.

2.1.   Models with age/stage structures

In this subsection, we will review two types of compartmental models with age structure: the models with delay terms that include varied brood stages and the compartment models without delay.

Motivated by biological stages of honeybee illustrated in Figure 1 and the honeybee-only model (2.2) from Kang et al. (2016) [40], Sweeney et al. (2019) [49] proposed the model (2.3) that includes brood larvae $B_l$, brood pupae $B_p$, and hive bees $H$ as below:

whose dynamics are similar to the model (2.4) by Chen et al. (2020) [22]:

where $B$ indicates brood stages, $H$ indicates adult stage, $r$ indicates daily egg-laying rate of queen, $\sqrt{K}$ indicates colony size at which brood survival rate is half maximum, $\alpha$ is the regulation effects of brood, and $d_b$ and $d_h$ indicate mortality rate of brood and adults, respectively. Both models assume that the new brood $\frac{r H^2}{K+H^2+\alpha B}$ or $\frac{r H^2}{K+H^2+\alpha B_l}$ is positively regulated by adults and negatively regulated by brood $B$ or $B_l$.

Both the three component model (2.3) and the two component model (2.4) exhibit strong Allee effects generated from collaborative effects of colony members, and colony survival when the egg-laying rate is above certain threshold. The dynamics comparison suggests that the additional age structure may not affect dynamics much. To further explore the effects of age structure, both age structure models are compared to the adult bee only model (2.2). It seems that the model with or without age structure has equilibrium dynamics, and additional age structure components in models may not contribute much in population dynamics. However, age structure may play an important role when additional mortality factors such as parasitism are included in the model. We will come back to this topic in Section 3.

Egg laying by the queen varies with season and colony conditions [29]. Chen et al. (2020) [22] included a variable egg laying to Messan's data-based model [29,50], and used the following periodic function ($r(t)$) (2.5) for queen's egg-laying rate:

where $\tau$ indicates the time in the juvenile period; $\psi$ is the time of the maximum laying rate; and $r_0$ is the baseline egg-laying rate. The study in Chen et al. (2020) [22] shows that seasonality may inhibit the survival of bee colonies, and at times prevent it from collapse.

Khoury et al. (2011) [21] developed a compartmental ordinary differential equation (ODE) model to determine the impact of forager death rates on colony growth and survival. The model includes age polyethism and division of labor where the adult worker population is divided into hive bees $H$ and foragers $F$:

This model assumes that the eclosion term $L*\frac{H+F}{H+F+\omega}$ (i.e., adult emergence from pupal state [51]) approaches the maximum egg laying rate $L$ as colony size gets large; and $\alpha -\sigma\left(\frac{F}{F+H} \right)$ is the recruitment function, $\alpha$ indicates the maximum rate that hive bees will become foragers when there are no foragers existent in the colony. Foragers can revert back to workers through social inhibition which is modelled by $\sigma\left(\frac{F}{F+H} \right)$ indicating that too many foragers can reduce the rate of hive bees becoming foragers. There is a critical threshold forager death rate beneath which colonies regulate a stable population size, and above which population decline and colony failure is predicted. The model includes recruitment of hive bees into the forager population when forager death rates are high. The recruitment of hive bees simulates precocious foraging that ultimately shortens the lifespan of adult bees, and if sustained for a long enough period, causes the colony to die. Their model provides more detailed observation of honeybee population dynamics, that the model shows that the adults switch between hives and foragers due to supply and demand, and provides a theoretical framework for experimental investigation of the colony failure. Therefore, much of the studies after this are based on the theories of Khoury et al. (2011) [21] and expanding their models to further study the reasons for the failure of honeybee colony caused by more factors [30,32,52,88,89,98,107,108,110].

Based on the model above (2.6), Booton et al. (2017) [52] proposed a hive-forager model with external stress to investigate the impact on the regulatory processes of recruitment to the forager class, social inhibition and the laying rate of the queen.

where $H$ is the in-hive worker population and $F$ is forager population. The model (2.7) predicts that constant density-dependent stress acting through an Allee effect on the hive can result in sudden catastrophic switches in dynamical behaviour and the eventual collapse of the hive. In addition, the model suggests that increased stress levels can be counteracted by a higher egg laying rate, lower levels of forager recruitment or lower levels of natural mortality of foragers. Increasing social inhibition cannot maintain the colony under high levels of stress.

2.2.   Honeybee population dynamics with food/nutrition effects

The interactions between food and population dynamics in honeybee colonies are very complex. Food collection is affected by the forager population size and food availability in the colony and the environment outside the hive. If food resources are limited in the colony, hive bees ($H$) will be recruited to become foragers ($F$) (precocious foragers). Precocious foraging shortens the lifespan of adult workers and affects the age structure of the colony and the amount of brood it can rear. These dynamics were captured in a model by Khoury et al. (2013) [30] when they extended their model (2.6) in Khoury et al. (2011) [21] by including food availability ($f$) and brood population ($B$) to explore how food availability interacts with the intrinsic demographic processes, and used time delay ($\tau$) in maturation term to present the days of adult bees emerging from pupation:

The study of the model (2.8) above suggests that both food availability and forager death rate have very strong influences on colony growth and development. It predicts that colony population eventually stabilizes at an equilibrium size determined by the forager death rate, and is no longer affected by increasing food availability, whereas food stores continued to increase. The dynamics of the model suggest that different forager death rates result in qualitatively and quantitatively different colony outcomes, which range from a stable population with an excess of food stores, to a stable population with limited food stores, to zero population with residual food stores. Their new model adds the influence of food to the original recruitment term, it is used to explore the dynamic of colony by interaction between the food storage and foragers mortality, and also provides a mathematical foundation for more complex performance of the honeybee colony [32,89].

Precocious foraging occurs under chronic stress from lack of food, disease, or loss of foragers. Precocious foraging reduces effective foraging and forager longevity thus leading to low food reserves and high mortality rates. In order to identify the most effective strategies to improve colony resilience, Perry et al. (2015) [32] extended the model in Khoury et al. (2013) [30] to the following one, the population of uncapped brood is denoted by $B$, hive bees by $H$, foragers by $F$, and the weight of stored food by $f$:

which assumes that the amount of trips taken and transitional survival ($T(a)$) are a function of the age at onset of foraging. The death rate of foragers is a function of the age at onset of foraging. The model (2.9) by Perry et al. (2015) [32] assumes that food storage depends on foraging trips which in turn depends on the age of onset of foraging. This assumption resulted in a breakdown in the division of labor and loss of the adult population, leaving only brood, food, and few adults in the hive. This work explains social processes that may drive rapid depopulation of a colony, and explores possible strategies to prevent colony failure.

In addition to the above-discussed food and the dynamics system of honeybee population, there are also regulatory mechanisms for task allocation via the "common stomach" (sometimes called "social crop"). The concept of a "common stomach", where all bees in a colony have similar amounts of food in their crops, has been incorporated in varied models for honeybees [53-56]. Honeybees communicate substance/information through a common stomach to store or retrieve pollen or nectar [53]. To study the economic benefits of pollen and nectar foraging strategies to both plants and honeybees under different weather conditions, Schmickl and Karsai (2016) [55] proposed and studied a detailed model that includes the nectar flows and nectar qualities offered by plants under distinct raining conditions. This model [55] simulates three strategies (adaptive, fixed and proactive) under different rains periods that affect protein stores and consumption rates in the colony. The adaptive strategy is similar to the principles that honeybees follow in that there is relatively a low-level of pollen foraging in good weather. During periods of rain when foraging activity stops, there is decreased pollen consumption and brood cannibalism. Vigorous pollen foraging occurs after the rainy period. Fixed strategy is that bees forage on pollen at a fixed recruitment rate and will not conserve protein when pollen is scarce. The proactive strategy is based on the fixed strategy, except that foraging rates enable an excess of pollen to be stored for periods of inclement weather. The model indicates that the adaptive strategy benefits bees and plants in that increased pollen foraging can compensate for lost pollen collection during the rain, while increasing pollination rates. Brood survival also was greater in colonies employing the adaptive strategy, but this depended upon the rain pattern with tropical patterns having less impact than continental. The model also indicates that the proactive strategy may cause colonies to optimize growth, but because it is at the expense of nectar collection the colony will be low on nectar stores causing starvation over the winter.

Schmickl and Karsai (2017) [56] extended the model of Schmickl and Karsai (2016) [55] to study the resilience and homeostatic regulation of honeybee colonies by including two "common stomachs" and the related foraging behavior. The first common stomach is the saturation of workers and brood with proteins that affect both the influx and the consumption of pollen. The second common stomach is the saturation of nectar stores, which affects the ratio of loaded to unloaded foraging bees returning to the hive, which in turn regulates the recruitment of nectar foraging bees. This model has a longer time scale than the previous study [55]. The study of Schmickl and Karsai (2017) [56] combined the simulation results with actual experimental data to explore the dynamics of adult forager populations, brood population, and protein and sugar stocks in the hive under varied environmental conditions.

The growth, development, productivity, and health of a honeybee colony is dependent upon fulfilling the nutritional demands of larvae and adult workers [57]. Vitellogenin (VG) is an egg yolk protein which is the primary source of amino acids [58], and it appears to be one of the most important regulators of immunity and longevity of honeybees [59-61]. The recent modeling work on the effect of honeybee nutrition on colony survival by Messan et al. (2018) [50] investigated how VG may affect colony dynamics (please see the detailed model in Appendix B model (B.1)). The non-linear differential equation system describes how population dynamics of colonies are influenced pollen stores and VG. Their study [50] focused on the effects of VG content per nurse bee and how it interacts with the essential demographic and allocation processes within the colony to influence colony growth. They showed that the size of both the brood and worker populations are directly dependent upon the increase of levels of VG titers in nurse bees; and the coexistence of both brood and worker populations is also dependent upon the quantity of food fed to the brood' (e.g., pollen collected and converted to VG and available foragers). In addition, the study in Messan et al. (2018) [50] explored how seasonal changes in pollen collection improves the prediction of long term colony dynamics.

The model proposed and studied in Messan et al. (2018) [50] is unique and different than most models in the literature that study the population dynamics of honeybees and nutrient stores. For example, Schmickl and Crailsheim (2007) [62] created a simulation model to study the population and resource dynamics of a honeybee colony and included the effect of division of labor in the hive. Their modeling approach is too complicated to be mathematical trackable. The model presented in Khoury et al. (2013) [30] has similar assumptions in the model of Messan et al. (2018) [50] with one of the main differences being that the assumption of the transition from nurse bee to forager. The model of Messan et al. (2018) [50] assumes that the transition from nurse bee to foragers increases by a transition rate that depends on the levels of VG in the nurse bee and not only by the absence of stored food and reduced social inhibition. In BEEHAVE by Becher et al. (2014) [33], they assume that the level of pollen and nectar stores in the colony affects the age at which workers initiate foraging activities while they do not take into account mechanisms changing at the molecular [60,63-65] and physiological level [60,66] as the model of Messan et al. (2018) [50].

The study by Messan et al. (2018) [50] presented basic but important assumptions that can help us understand and have greater insights in the complexity of honeybee population dynamics given their nutritional status. Their model can be extended by limiting the influx of pollen in the late summer to determine the effects on adult lifespan and colony survival. The effects that brood cannibalism or breaks in brood rearing at different times of year can have on colony survival also can be studied through simulations. Other mechanism to be considered is the transition back from forager to nurse bee under certain environmental conditions [67], and the interplay of diseases or infections, such as Nosema ceranae which is known to alter VG levels, and therefore alter normal age polyethism causing colony imbalance [68,69].

Social insect groups have evolved a collective immune defense against parasites. These "social immune systems" are built on the cooperation of individual members in the colony to fight the spread of disease. Collective defense can prevent the entry and spread of parasites. In the study of the social immune system, it is necessary to carry out a complex integration of information about the parasites and the internal state of the colony [70]. Laomettachit et al. (2021) [71] proposed a system (please see the detailed model in Appendix B model (B.2)) that describes the spread of infection within the colony when honeybee members engage in the trophallactic activity to distribute nectar to study the social immunity in honeybee society. Their model considered the social segregation of worker bees and a hygienic response. The segregation limits the forgers that exposed to pathogens interacting with bees in the colony, and hygienic response is the healthy nurses exterminate infected bees to control the infections spreading. Their study suggested that the segregation is the first line of defense in pathogens into the colony, and the hygienic response is the second line of defense. Their study also found the egg-laying rate is a critical factor in maintain healthy colonies against infection. In summary, they suggested the egg-laying rate, that is decreasing or pause combining the infections in the early spring, can damage the social immunity, and even cause colony collapse.

3.   Honeybee and parasites/virus population systems

Honeybee colonies provide a favorable environment for parasites and pathogens because of the warm temperature and the high concentration of hosts, which may contribute to honeybees disappearing in large numbers across the globe [72]. Studies showed that most of the losses could be generally attributed to Varroa mites who have been implicated as the main culprit in dying colonies [73-75]. In the following subsections, we will provide review of mathematical modeling that includes Varroa alone and with virus transmission. We also discuss a model of the effects of Nosema on colony dynamics.

3.1.   Honeybee-mite systems

Varroa mites are a parasite of both adult bees and the developing brood and can shorten the life span of their hosts [79]. Adult female mites undergo two phases in their life cycle, the phoretic and reproductive phases. During the phoretic phase, female Varroa feed on adult bees and move from bee to bee as they pass one another in the colony. Mite reproduction can occur only if honeybee brood is available [24]. The reproductive phase begins when a mature female leaves her adult host, enters a brood cell containing a larva shortly before it is capped, and sequesters herself in the bottom of the cell. Egg-laying commences about 60 hours after a cell is capped, and both mother and offspring feed on the host's fat body [80]. Mature offspring mate within the cell, but only mature females survive outside the cell. In the average temperate climate, mite populations can increase 12-fold in colonies which have brood half of the year and 800-fold in colonies which have brood year-round. The period of mite population increase continues through the spring and summer and peaks in the fall when brood rearing is nearly done. This makes the mites very difficult to control, especially in warmer climates where colonies maintain brood year-round [28].

Sweeney et al. (2019) in his undergraduate honor thesis [49] developed a realistic five-dimensional honeybee-mite interaction model with age structure in both honeybee and mite (see Figure 2, and the detailed math model in Appendix B model (2.3)). This model includes three delay terms: one delay term $\frac{e^{-d_l\tau_l}rH(t-\tau_l)^2}{K+ H(t-\tau_l)^2}$ describes the transition from larvae to pupae, and the second one $\frac{e^{-\int_{t-\tau_p}^t \left[d_p+\frac{\alpha_p M_p(s)}{a+B_p(s)}\right] ds}rH(t-\tau_p-\tau_l)^2}{K+ H(t-\tau_p-\tau_l)^2}$ describes the transition from pupae to adult; and the third delay term $c\alpha_p \frac{B_p(t-\tau_p)}{a+B_p(t-\tau_p)}M_p(t-\tau_p)$ describes the newborn mites coming from the parasitism of honeybee brood $B$. The modeling approach of those delay terms is adapted from Chen et al. (2020) [22].

The advantages of the honeybee-mite interaction model proposed and studied by Sweeney et al. (2019) [49] are that it includes 1) three stages in honeybees: larvae, pupae, and adult bees; 2) two stages in mites: phoretic stage and reproductive stage; 3) mites at the phoretic stage are parasitizing adult bees $H$ which contributes to the additional mortality rate; and 4) mites at the reproductive stage are parasitizing pupae $B_p$ to reproduce. The simple analysis and the simulations suggest that parasitism combined with the stage structures can destabilize the system causing fluctuating dynamics.

The disadvantage of the proposed model by Sweeney et al. (2019) [49] is the difficulty of rigorous mathematical analysis. To explore the mechanisms that generate instability of the honeybee-mite interaction dynamics, Messan et al. (2021) [29] simplified the model of Sweeney et al. (2019) [49] to the following three-dimensional honeybee-mite interaction model:

where the brood $B$ consists of both larvae $B_l$ and pupa $B_p$ in model (B.4), and no stage structure for the mites. When there is no mite in the system, i.e., $M = 0$, the model (3.1) reduces to the two-stage honeybee model (2.4) studied by Chen et al. (2020) [22]. Both the work by Chen et al. (2020) [22] and the model (B.4) without mites by Sweeney et al. (2019) [49] suggest that time delay does not affect dynamics, while the study of models (3.1) and (B.4) shows that the large time delay in the presence of parasite mites (i.e., $M > 0$) could destabilize population dynamics and drive the colony to collapse. In addition, both models [22,29] use real data to show the effects of seasonality Eq (2.5), and their study suggested seasonality is a key factor for colony's collapse or survival.

Both age-structured honeybee models without mites in models (3.1) and (B.4) suggest that the delay or age-structure terms may not play significant roles in population dynamics. Kang et al. (2016) [40] proposed a simple honeybee only model (2.2) that has been extend to the honeybee-mite interaction model without age-structure.

where $S_h$ indicates the population of honeybees (includes all brood and adults) and $S_m$ indicates the population of mites. Their model can have bistability and go through the subcritical Hopf-bifurcation that leads to the extinction of colony. The study shows that parasitism could destabilize population dynamics of honeybees.

Messan et al. (2017) [48] extends the model (3.2) to the following simple two-patch honeybee-mite interaction model (B.5), the detailed math model in Appendix B (B.5), to explore the mite migration effects on the population dynamics of honeybees and mites. This study suggests that migration between two patches can generate extremely complicated dynamics including multiple attractors, and plays the key role in determining whether honeybee colony survives or collapses. For example, the study shows that the median migration rate can potentially enable mites and honeybees to coexist, while a high migration rate can rescue one colony in two patches where both were going to collapse.

3.2.   Honeybee-virus systems

Viruses are one of the main threats to the health and survival of honeybees. More than 20 viruses have been identified to infect honey bees worldwide [33]. The viruses most commonly detected are: Acute Bee Paralysis Virus (ABPV), Black Queen Cell Virus (BQCV), Deformed Wing Virus (DWV), Israel Acute Paralysis Virus (IAPV), Kashmir Bee Virus (KBV), Sacbrood Virus (SBV) and Chronic Bee Paralysis Virus (CBPV) [81]. For a complete discussion of symptoms associated with each virus and the mode of transmission see [82]. Though several different viruses can infect worker bees throughout the year [83,84], by the fall, DWV often is the predominant virus detected [85]. This is because of differences in the quantitative dynamics among the viruses transmitted by Varroa. Those dynamics can generate a shifting succession of virus infections that ultimately leaves DWV as the predominant infection [86]. High levels of DWV in the fall are a harbinger of colony loss over the winter [87].

Kang et al. (2016) [40] extended their honeybee only model (Eq (2.2)) to the honeybee-virus model (3.2) shown below:

which assumes that infections only spread horizontally with additional mortality $\mu_h$ in infected honeybees, new adults bees are all healthy, and the virus-infected bees have reduced performance that is denoted by $\rho$. The dynamics of model (3.3) seem to be relatively simple: i.e., no limit cycle. The study of Kang et al. (2016) [40] shows that 1) honeybees can persist with proper initial conditions when the virus transmission rate among honeybees $\beta_h$ is not large, 2) both honeybees and the virus can coexist if $\beta_h$ is in the medium range, and 3) the large virus transmission rate $\beta_h$ among honeybees can drive the colony to extinction.

Kribs and Mitchell (2014) [31] investigated the scenario that CCD is provoked by a transmissible pathogen or contaminant introduced by foraging activities to understand how colony collapse might occur as a perturbation of normal hive dynamics. In the model, $H$ indicates hive bees, $F$ indicates forager bees, and $I$ indicates infected bees.

where virus infection brought into the colony is by the foragers and spreads in the hive via both direct bee-to-bee transmission and indirect (via contaminated plant vectors) transmission. The disease dynamics assumed in this paper are very general, and is not particularly discussed in relation to any vector or virus. But this model still helps us to understand the dynamics of healthy colonies and colonies collapsing because of CCD, especially the impact of accelerated forager recruitment on the collapse of the colony, and the balance between recruitment rate and laying egg rate.

Betti et al. (2016) [88] proposed a PDE model that takes into account age structure and seasonal effects on the dynamics of a colony infected with a pathogen.

where $a$ is age and $t$ is time, for example, $H_S(a, t)$ means susceptible hive bees of age $a$ at time $t$ and $\mu(a)$ is natural death rate. The model assumes that the hive bees and foragers spread the infection at a rate ($\beta$) between or within classes. Figure 1 shows the processes between hive bees and foragers are social inhibition and transition. The term $u(a)$ represents similarly in the recruitment of Khoury et al. (2011) [21]: $u(a) = \alpha(\frac{a}{a+k})^2(1-\frac{\sigma}{N}\int{(F_s+F_I)da})H_V(a-a_R)$ with $N$ being the total number of bees, $\alpha$ being the maximum rate of recruitment, $k$ being that recruitment will be half the maximum rate of recruitment and $\frac{1}{\sigma}$ being the maximum proportion of bees that are foraging at any given time. The term $H_V(a-a_R)$ is the Heaviside function to reflect that recruitment cannot begin before age $a_R$. When determining the initial value, the food factor is introduced, $H_S(0, t) = LS$, where $L$ indicates the daily egg-laying rate of the queen, and $S$ is the survivability function that includes food to determine how many eggs can survive.

The PDE model (3.5) shows that the age-dependent recruitment rate affects the dynamic of disease within the colony. The natural age-structure has a certain control effect on disease transmission. In other words, if there is an interference with the age-structure, it may also increase the extent of infection. The modeling work by Chen et al. (2020) [22] and Messan et al. (2021) [29] show that time delay alone may not have huge impacts on the dynamic of bee-only model, while delay combined with mites may lead to colony collapsing. The model by Betti et al. (2016) [88] showed that the maturation time from hive bees to foragers affects both colony dynamics and infection dynamics if a pathogen is present in the colony. The age dependent model in Betti et al. (2016) [88] captures the dynamics of 'spring dwindling' which occurs due to the fact that worker bees emerged from the previous fall are dying faster than the emergence of new workers in the spring. In addition, the model of Betti et al. (2016) [88] also includes seasonality by adding a transition model after end-of-winter to expose the process of colony aging from winter months to active young bees. The work of Betti et al. (2016) [88] indicates that colonies are particularly vulnerable to collapse during spring dwindling especially if disease or other stress factors are present; and the declining population of bee colony in the early spring is the natural consequence due to seasonal vulnerability.

3.3.   Honeybee-mite-virus systems

The host-parasite relationship between honeybees and Varroa mites has been complicated by the mite's close association with a wide range of honeybee viral pathogens. While Varroa parasitism can weaken honey bee immunity and deplete nutrient reserves in their host, the most serious threat is transmission of viruses [82,90,91]. Varroa transmit numerous viruses during feeding on developing and adult bees. Mathematical models are powerful tools that could help us obtain insights on how mites and/or viruses contribute to the collapse of honeybee colonies.

Several models have captured the mechanics of virus infections in colonies in the presence of mites [40,92]. Kang et al. (2016) [40] proposed and studied a honeybee-mite-virus model (B.6) shown in Appendix B model (B.6)) that incorporates parasitic interactions between honeybees and Varroa mites with terms describing virus transmission between bees and mites at different lifestages. The honeybee-mite-virus model (B.6) incorporated colony's internal organization (i.e., division of labor) that can generate Allee effects, and five virus transmission terms between honeybees and mites at different stages of Varroa mites. The routes of transmission include from: honeybees to honeybees, adult honeybees to phoretic mites, honeybee brood to reproductive mites, reproductive mites to honeybee brood, and honeybees to phoretic mites.

The model (B.6) captures the synergistic effects of parasitism and virus infections on honeybee population dynamics and persistence of virus in the colony. The components included in the honeybee-mite-virus model (B.6) allow us to explore the dynamics of the: honeybee only model (2.1), honeybee-mite model (3.2), honeybee-virus model (3.3), and honeybee population with both mites and viruses. The model revealed that colony survival is strongly affected by initial conditions (i.e., colony size) due to Allee effects experienced by the honeybee population. High adult to brood ratios enable coexistence of Varroa and virus in the population. However, low adult to brood ratios have destabilizing effects on the system which generate fluctuating dynamics that lead to a catastrophic event where both honeybees and mites become extinct. This model outcome is similar to colony losses over the winter that are commonly seen under conditions of high Varroa and virus (particularly DWV) titers in the fall [87]. Additionally, persistent viral infection in late winter and early spring can lead to high mortality of worker bees, and an inability to rear sufficient amounts of brood (Allee effects) to enable the colony to build in the spring.

The model (B.6) by Kang et al. (2016) [40] applies to disease/virus in the general fashion. Dénes and A. Ibrahim (2019) [93] studied a honeybee-mites-virus system with susceptible bees, bees infected with mites and bees infected with viruses (please see the detailed math model in Appendix B model (B.7)). Their results suggest that in order to eradicate disease/mites, the related transmission rates needs to be decreased. In the following subsections, we will provide reviews on literature work that have been applied to specific viruses such as ABPV and DWV.

3.3.1.   Honeybee-mite-ABPV/DWV/AFB systems

Acute Bee Paralysis Virus (ABPV) and Deformed Wing Virus (DWV) are transmitted by Varroa[82,94]. ABPV is deadly in the sense that brood infected with it do not metamorphose to the adult stage. DWV causes death of the pupae and reduced longevity of adult bees [90,91]. And American Foulbrood (AFB) is a fatal bacterial disease of honey bee brood, and AFB is considered a very destructive and lethal brood disease [95].

Sumpter et al. (2004) [92] proposed a model that used SIR-like framework to study the effect of a honeybee colony dynamics with Varroa mites and virus. They studied ABPV and DWV separately. The model assumed the total number of mites is constant, therefore the death and birth rates of mites are not important in their model. Since viruses can exist in different stages of the life cycle of bees, seasonal effects are included and reflected in the values of different parameters for each season. A critical mite load that a colony can tolerate was estimated, and it was determined that as the number of mites exceeds the critical mite load, the disease free equilibrium becomes unstable, which means an epidemic occurs. Recommendations for reducing colony losses generated from the model were that mite and virus levels need to be controlled in the summer, and delaying control until fall increases the likelihood of colony loss.

The process of bee eclosion depends on how well the brood is cared for during development. For that purpose, a sufficient number of healthy adult worker bees are required. The relationship between healthy adult workers and eclosion was modeled by Eberl et al. (2010) [43] by extending the ABPV model by Sumpter and Martin (2004) [92]. They formulated brood maintenance as a Sigmoidal Hill function which introduced Allee effects and led to the collapse equilibrium being locally stable which was not present in Sumpter and Martin (2004) [92]. The model predicted colony death at points where the uninfected population drops below a certain threshold that is dependent on the brood maintenance terms. We would like to point out that Eberl et al. (2010) [43] also includes the seasonality through simulating different parameters values in four seasons. Their simulations suggested that an infected colony may collapse after several years in the spring. This result is in line with the work by Betti et al. (2016) [88].

Since the transmission of ABPV depends on the Varroa mites, it is important to track the mite population. Therefore, Ratti et al. (2012) [44] (see also Devillers (2014) [97]) adopted the framework of Eberl et al. (2010) [43] and introduced logistic growth for mites (please see the detailed model in Appendix B model (B.8)). As a result, mite-induced death rates have also been considered. Ratti et al. (2012) [44] assumed that the carrying capacity ($\alpha (H_S+H_I)$) is dependent on the total bee population which is a proxy for the brood size. For a model with logistic growth of mites, infestations may lead to a die-off of the colony or to bee-mite coexistence. The bee colony that otherwise coexists with mites might start dying when the virus is introduced. This die-off may take several years. Whether or not the bee population can fight off a virus epidemic depends on model parameters describing disease transmission and how fast infected bees die, as well as on the size of the mite population in the stable bee-mite coexistence equilibrium. They studied the model in constant-coefficients analytically and then investigated in computer simulations whether or not these findings carry over to the more general transient case in which the parameters are taken as seasonal averages. The model did not take into account seasonal or yearly fluctuations in conditions.

Research has shown that seasonality has the potential to change the dynamics of infectious diseases since it creates complex fluctuations in the host and vector populations. Most infectious diseases change with temperature, rainfall or other seasonal effects most of which are periodic and predictable. The effects of seasonality on honeybee and thus mite populations have been reported in the literature. Ratti et al. (2015) [45] adopted the framework of Ratti et al. (2012) [44] to study the seasonal effects on the mite and virus-infested colony by considering the coefficients to be periodic functions of time with a period of 1 year. They used Floquet theory to analyze the stability of virus-free periodic solutions. For controlling vector-borne disease, it is required to control the vector. Therefore, Ratti et al. (2015) [45] studied the effect of varroacide application by introducing additional death rates ($\delta_i$) in the bees and mite equations. In addition to the existence of periodic endemic solution in the absence of virus, it was found that the success of varroacide treatment depends on the cumulative efficacy over the year, rather than a specific time course. They also found that in order to keep disease under control, it is not necessary to eradicate mites completely.

Forager bees from diseased colonies experience reduced homing and flight capabilities which makes it important to track forager populations in models as in Khoury et al. (2011) [21]. Ratti et al. (2017) [47] constructed a model combining their ABPV model (from Ratti et al. (2015) [45]) with a division of labor component from Khoury et al. (2011) [21] to study forager loss. In their model (B.9) (please see the detailed model in Appendix B model (B.9)), they assumed that the total bee population was divided into hive bees and foragers, all newly emerged bees are uninfected, and only horizontal virus transmission occurs vectored by Varroa. Ratti et al. (2017) [47] categorized the population into uninfected hive bees ($H$), uninfected forager bees ($F$), infected hive bees ($H_I$), virus-free mites ($M_S$) and virus-carrying mites ($M_I$) (see Eq (B.9)). Ratti et al. (2017) [47] The model predicted that mite and bees can co-exist without treatment in virus-free conditions. If virus exists and a varroacide is applied causing a decline in virus-carrying mites that exceeds the mite birth rate, the disease-free equilibrium can be stable. Simulation indicates that infected colonies eventually collapse in the absence of Varroa treatment.

By following the work of Henry et al. (2012) [98], Ratti et al. (2017) [47] incorporated the increased forager mortality due to exposure to pesticides as a result of which the bees fail to return to the hive. The idea was to study the dynamics of the colony under multiple stressors (mites, virus, and pesticides). Though they investigated the effect of increased forager loss particularly due to pesticides, the framework can be easily adapted to forager loss due to other factors. They also suggested that a honeybee colony that experiences periods of increased forager loss early in the season can withstand these losses and recover to (almost) normal strength, if the losses remain below a defined threshold. If the losses exceed this threshold the colony fails rapidly. The level of external forager loss that can be tolerated is only slightly higher for healthy colonies compared to colonies with Varroa and ABPV.

Britton et al. (2021) [99] proposed the first model to study the dynamics of DWV that classified bees as covertly infected ($H_X$), chronically overtly infected ($H_Y$), and acutely overtly infected ($H_Z$), and mites as covert infected ($M_U$) and overt infected ($M_V$) based on the infection levels. The model equations for spring, summer, and fall are given by

The mite model for winter is given by

Moreover, they modeled the spontaneous transition from one infected class to the other, and included logistic growth for mites for all seasons except winter [44,45,47]. Since there is no brood in winter, the mite population dynamics is much simpler and hence the mite equations are reduced to having only negative linear terms (see Eq (3.7)). Mites are no longer considered to be a mechanical vector; the virus actually replicates in the body of the mites. For the bee dynamics, they factored in the vertical as well as horizontal transmission. When possible, they calculated the $R_0$ using next generation matrix. The model did not incorporate seasonal variations analytically or computationally. The focus of this paper was on gaining the insight into the disease dynamics by computing the reproduction number using weighted directed graphs with the next generation method. They found that in the absence of mites, the condition for overt infection leads to unrealistic results. They also found that the overt infection can be maintained either by the transmission of overt virus between adult bees (for unrealistic parameters) or by vertical transmission in mites coupled with their interactions with acutely overtly infected bees. This demonstrates the importance of considering transmission involving mites and chronically overtly infected bees.

Datta et al. (2013) [100] studied the spreading of American foulbrood (AFB) by using an SIR model and estimated parameter values by using Markov chain Monte Carlo (MCMC) scheme based on data of AFB epidemic in the island of Jersey bee colony during the summer of 2010. They found that both distance-based and owner-based transmissions contribute to the spread of AFB. And the simulations of the stochastic SIR model, they showed that the earlier inspection, the better control of the epidemic, and the higher possibility of AFB extinction.

Jatulan et al. (2015) [101] build a compartmental ordinary differential equations system (please see the detailed model in Appendix B model (B.10)) to study American Foulbrood (AFB) infection within a honeybee colony. Their model is different from others in the sense that it includes how inside environment such as contaminated cells and outside environment infect bees. Their study suggested that AFB infection eventually leads to colony collapse; and infection thresholds were predicted based on the stability of the equilibrium states.

3.4.   Honeybee-Nosema systems

Nosema ceranae is an intracellular parasite of honeybees. However, its effects on colony mortality remains poorly understood. However, Nosema will cause precocious foraging and learning and memory deficits [105]. Infected bees defecate and deposit spores in the hive on the wax comb and hive bees can ingest them while cleaning comb cells. During colony confinement in the winter, bees will not defecate and can have millions of spores their digestive tract. There is also evidence that Nosema can be spread via oral-oral transmission, through trophallaxis [106]. Hive bees mouth parts get infected with Nosema while cleaning the hive and it gets transmitted to the uninfected hive bees during trophallaxis.

Betti et al. (2014) [89] presented a generic model with a working example of Nosema to study the dynamics of the spread of disease within a bee colony and the effects on colony demographics. The model for active season (spring, summer and fall) is following:

where $H_S$ and $F_S$ indicate the susceptible hive bees and foragers respectively, $H_I$ and $F_I$ are the infected hive bees and foragers respectively, and $f$ is food. Parameter $L$ is the egg-laying rate, $S$ is the proportion of eggs that survive to adult bees, $R$ is the proportion of recruitment from maturing hive bees to foraging duties. And all $\beta_{i, j}$ is the contact between stage $i$ bees and stage $j$ bees. Parameter $m$ indicates natural death rate, $c$ os the food collecting rate, $\gamma_A$ and $\gamma_L$ are the food consumption by adults and larvae, and $N = F_I+H_I+F_S+H_S$. The age structure part of this model is based on Khoury et al. (2011 & 2013) [21,30]. The study considered a direct transmission from bee to bee (e.g., by feeding). Adult bees will contract Nosema either from eating food contaminated by infected bees, or while removing fecal matter from the colony. The model also considers seasonality; winter is one season and the remaining (spring, summer and fall) are combined as another season, that is the model we showed above. The parameters are taken to be constant in the two seasons. And in the winter season, they proposed a corresponding model which is no newborn term ($L = 0$) because Queen stop laying egg in the winter, no recruitment term ($R = 0$) because all bees stay at the colony, and no new food ($c = 0$) because no foragers go outside to collect food. The model indicated that transmission rate and death rate due to disease are the key factors that determine colony survival. In comparison, the transmission rate has a greater impact on survival than the mortality. The model also predicted that a time interval of approximately 20 days before the onset of winter is when the colony is most affected by the onset of infection. From this point of view, the transition between winter and non-winter seems to be an important aspect of seasonality and colony survival. Therefore, the transition model between the winter and non-winter is also worth learning, that Betti et al. (2016) [88] proposed 21-day transition equations we discussed before.

Petric et al. (2016) [107] specifically developed a model for Nosema ceranae (please see the detailed model in Appendix B model (B.11)) based on the assumption that infected bees deposit spores to an environmental reservoir in the hive created by defecation and hive bees ingest those spores while cleaning comb cells. The model incorporate two stressors: 1) infection of a honeybee colony with Nosema ceranae following Khoury et al. (2011) [21] and 2) increased forager loss. Seasonality was incorporated in the same way as in Ratti et al. (2015 & 2017) [45,47]. Due to lack of data at that time, they performed sensitivity analysis of the model instead of the complete parameterization. The model indicates that the colony can survive if each stressor alone is applied but the combination of infection and increased forager loss causes the colony to fail.

Muhammad and Eberl (2020) [108] extended the model of Petric et al. (2016) [107] by including direct transmission (e.g., by oral-oral route) following Betti et al. (2014) [89] (please see the detailed model in Appendix B model (B.12)). The model indicated that including direct transmission of a disease in a time-periodic model of a discrete structured population can induce what appears to be a chaotic transition from an endemic state to colony failure in a small range of the transmission rate. They also concluded that the observed chaotic transition is heavily influenced by the Holling-Type III functional form of the eclosion terms.

Not only Nosema can spread within a colony, but also between colonies, therefore there is a competitive relationship between different species of Nosema. Natsopoulou et al. (2015) [109] presented an ODE model that focused on Nosema by investigating the within-host competition between two different Nosema species (Nosema apis and Nosema ceranae) as below:

where $A$ indicates Nosema apis and $C$ indicates Nosema ceranae. They found that the competition is asymmetric, and the order of infection determines the outcome of the interspecific competition. It was also found that the prevalence of Nosema ceranae was higher causing the growth of Nosema apis to be inhibited.

Since Nosema has adverse affects on the food consumption and foraging ability, Comper and Eberl (2020) [110] constructed a model (B.13) (please see the detailed model in Appendix B) to study the effect of food in a colony infested with Nosema. They categorized food as harvested food i.e., honey and externally supplied food i.e., sugar syrup. They considered sugar syrup to be given to bees before winter as a way to make up for the honey that was harvested. They used the framework for division of labor from Khoury et al. (2011) [21] and indirect Nosema transmission from Petric et al. (2016) [107]. Seasonality was taken into consideration and was adapted in the same manner as in Petric et al. (2016) [107], Ratti et al. (2015) [45], and Ratti et al. (2017) [47]. Comper and Eberl (2020) [110] found the same three possible outcomes as in Petric et al. (2016) [107]. The honey yield was found to be directly proportional to the Nosema strength in the colony. One of the main parameters was found to be the supplementary sugar syrup. They noticed that if the supply of sugar syrup falls below a certain level, the colony does not survive. Honey loss cannot be made up by using supplementary sugar syrup. Honey loss due to disease can be repaired by applying a Nosema treatment if the level of infection in the colony is moderate. However, treatment does not have any effect if the infection levels are high.

4.   Pesticide effects on honeybee population

Pesticides are an environmental stressor causing population reduction. Honeybees are exposed to pesticides while foraging pesticide-treated crops. Foragers bring back the pesticides to the hive which stays in the wax comb and in the pollen and nectar stores. Thus, honeybees do not get exposed only during the flowering period but much longer [111]. The effects of pesticide exposure include significant decrease in survival of foragers, reduced flight performance, compromised immunity, difficulty in protein digestion, and bee feeding [18].

Pesticides also have adverse effects on the memory, learning ability and navigational skills of the foragers. Henry et al. (2012) [98] study the foragers homing failure due to pesticides exposure. In addition to real experiments, they also combined model simulations by Khoury et al. (2011) [21] with their experimental data to predict the dynamics of colonies. They assumed the forager mortality rate includes daily natural mortality and additional death rate because of post-exposure homing failure. Henry et al. (2012) [98] found regardless of the queen's egg-laying rate, the population exposed to the treated nectar will decline during pollination, and will hardly recover afterwards. This paper is an excellent application of mathematical models and helps to give parameter values closer to reality for the mathematical model.

As mentioned in Section 3.3.1, Ratti et al. (2017) [47] following the work of Henry et al. (2012) [98] studied the effect of pesticides via loss of foragers in a colony infested with mites and ABPV. The objective was to investigate how the colony dynamics would change if a colony, that is already under stress (due to mites and ABPV), experiences increased forager loss (e.g., due to pesticides). They studied the first 3 months of a beekeeping season (i.e., the spring season) encompassing the blooming period of crops that are treated with pesticides and visited by bees. Increased forager loss was formulated in the model by using a linear death term (see model (B.9)). The study concluded that a honeybee colony that experiences periods of increased forager loss early in the season can withstand these losses and recover to (almost) normal strength, if the losses remain below a threshold. If the losses exceed this threshold the colony fails.

Magal et al. (2019) [112] developed a honeybee difference equation model that tracks the foragers leaving and entering the hive each day to study the effect of pesticides leading to homing failure. The number of uncontaminated ($U$) and contaminated ($C$) forager bees each day is given by

where $F(U, C)$ and $G(U, C)$ represents the daily change in uncontaminated ($U$) and contaminated ($C$) foragers respectively and are given by

The function $V(U, C)$ indicates daily total uncontaminated and contaminated foragers returning home and is given by

where $p$ is the fraction of contaminated forager bees that succeeded in returning to the hive at the end the day, $\beta$ is the rate of new foragers per day, $\mu$ is the natural morality rate of forager bees and $\alpha$ is the rate of contamination. One of the key factors of colony survive is the critical forager population viability threshold. If the population is below the threshold, the colony dies off rapidly, otherwise the population stabilizes. This dynamic of their model is incorporated as an Allee effect, that a sufficient number of young hive bees is necessary to offset the forager loss. Magal et al. (2019) [112] discussed three equilibria (one trivial and two interior) by varying the parameter $p$. By varying this parameter, the model captured sublethal effects of pesticides leading to colony loss due to the reduction of homing ability by foragers.

Magal et al. (2020) [113] extended their previous model in Magal et al. (2019) [112] (see model (4.2)) by including spatial components to study colony collapse due to pesticides. The model consists of differential and difference equations for the spatial distribution of uncontaminated and contaminated foraging bees. In addition to considering whether contaminated foragers can successfully return to the colony, they also considered a new condition of whether foragers will return to the previously visited foraging site. They assumed that foragers revisit the locations based on their memory and not on the waggle dance. Furthermore, they assumed that a proportion of foragers will return to the hive and start over the next day without the memory from the previous day. Included was a proportion of foragers that use their memory from the previous day. Their model focused on the density of uncontaminated ($u (t, x, y)$) and contaminated forager bees ($c (t, x, y)$) at time $t$ and location $(x, y) \in \mathbb{R}^2$:

where $\mu(x, y)$ indicates the mortality rate of the bees, and $\alpha(x, y)$ indicate the rate of contamination of the forager bees by pesticides. The parameter $\varepsilon$ represents the diffusion rate of uncontaminated foragers. They found that a high fraction of contaminated foragers that fail to return home brings the population below the thresh old and leads to colony collapse. Also, if the fraction of all foragers that return to previous foraging sites is high and spatial variation is present in the contaminated environment, foragers who visit contaminated sites multiple times have a higher probability of becoming contaminated, and colony collapse ensues.

Rumkee et al. (2015) [114] used BEEHAVE model to study the effects of pesticides on different life stages of honeybees. They explored the dynamics of bees through five stresses: increasing the mortality of larvae, hive bees and foragers, and reducing the egg laying rate. They set control mortality by multiplied and daily background mortalities by different percentages. Included was a proposed LIS50 to indicate the level of stress leading to 50% colony mortality. Small LIS50 means this pressure exerts high. Their results suggested that the mortality of adults has a greater impact on the colony than the mortality of larvae or the egg-laying rate. Different stresses also have corresponding levels of influence on the colony dynamics under the seasonality, such as forager mortality in winter and in-hive mortality in summer.

5.   Computational models

The earliest computational model, DeGrandi-Hoffman et al. (1989) BEEPOP [19], included fundamental factors for colony growth including egg laying, brood development, worker aging from house bees to foragers, and number of foraging days before the worker was removed from the colony. The BEEPOP model incorporates weather conditions (temperature, photoperiod, wind speed and rainfall) into colony population size predictions. The model tracks numbers of workers and drones in all life stages, by treating the number of eggs a queen lays each day as a cohort. The cohort is aged from egg to adult forager and then removed from the colony population after a defined number of days of foraging. The daily number of eggs the queen lays is determined by temperature, daylength and colony size. More eggs are laid in the spring and summer than fall and winter in comparably sized colonies because temperatures are warmer and daylength is longer. All other conditions being equal, more eggs are laid by queens in large colonies compared with small colonies. The number of bees in all cohorts are summed daily and reported as total numbers of eggs, larvae, pupae and adult bees.

Factors affecting the survival of bees in the various lifestages were not included in the original BEEPOP model, with the exception of forager longevity. The number of days that an adult worker could forage was initialized at the start of a simulation. If foragers had a short lifespan (e.g., 4–6 days), this would simulate stress from factors such as poor nutrition, disease, or distant nectar and pollen sources. Though the BEEPOP model is simple in its structure, it did generate realistic predictions on colony population dynamics. For example, the incorporation of weather conditions into predictions of daily egg laying rates enabled the model to capture colony population dynamics in different latitudes and climates. Sensitivity analysis revealed that the queen's egg laying potential and worker longevity exert significant effects on colony population growth, and this has been supported through field studies by others [115].

The effects of Varroa mite parasitism were included in a later version of the BEEPOP model [28]. VARROAPOP considered the growth of the Varroa population in a colony and its impact on colony dynamics due to reduction in worker longevity. The model assumes that adult workers parasitized during development have reduced longevity as adults [94,116,117]. The reduction in longevity is a function of the number of foundress mites (i.e., reproductive female mites) that enter a worker cell. Workers that are parasitized by multiple foundresses during development have shorter lifespans than those that are singly infested. During times of colony growth in the spring and summer, there often are more brood cells to parasitize than mites to infest them, and so most cells are infested by a single foundress. However in the fall, brood rearing declines. Multiple infestations of brood cells can occur because the mite population has increased throughout the spring and summer, and there are fewer brood cells in the fall to infest. The multiple infestations of brood cells in the fall and the resulting reduction in worker longevity causes the colony age structure to be skewed to short-lived adult bees. If the population of short-lived adults is sufficiently large, the model predicts that the colony will perish in the spring when flight weather resumes. The death of colonies over the winter and in the spring due to mite infestations in the fall is a common occurrence [87] that is captured in VARROAPOP simulations.

The VARROAPOP model has been useful in the development of control strategies for Varroa. The model can simulate the application of miticides to reduce the population of phoretic mites. Specifics of a miticide application such as the date when the application occurs, period of effectiveness, and the percentage of phoretic mites that will die from exposure (this value if reduced over the period of effectiveness for the miticide) can be entered into a simulation. Conducting simulations where miticides were applied at different times of year under weather conditions occurring in different geographic locations generated recommendations on the optimum time and number of miticide applications required to control Varroa populations and limit their impact on colony survival. VARROAPOP predicted that a single well-timed miticide treatment in mid-August to mid-September could adequately control Varroa and avoid mite populations that could cause a colony's demise. The model also was used to develop a commercial miticide by determining the minimum period when the active ingredient needed to be effective and the percentage of phoretic mites that had to be removed from the population.

VARROAPOP was constructed to include the migration of mites into colonies on foragers. However, the impact that mite migration could have on actual Varroa population growth in colonies was not fully realized until comparisons were made between Varroa population growth occurring from reproduction alone, and actual mite population growth especially in the fall [28,48,118]. Actual mite populations beginning in September were significantly higher than those predicted from reproduction even though mite numbers were small when the colonies were established and accurately predicted by VARROAPOP throughout the summer. The disparity between actual mite numbers and those predicted by the model occurred when the frequency of foragers with mites on their bodies collected at colony entrances increased.

Most recently, VARROAPOP was modified to track nectar and pollen collection and use by the colony based on larval and adult populations.Other modifications enable users to assign mortality rates for any life stage for a defined period. The age of first foraging for worker bees also can be changed for defined intervals during a simulation. Multiple intervals can be described for a single simulation so that the effects of increased mortality due to various stress factors can evaluated. The feature enables simulations to include effects of nutritional stress (altered brood mortality and age of first foraging of adult workers), swarming [119] and pesticide exposure (increased mortality at any lifestage and reduced forager lifespan) at [120]. Conducting simulations using Monte Carlo and sensitivity analysis techniques to determine the impact of foliar, seed or soil pesticide applications indicated that queen egg laying potential and forager lifespan determined the survival of colonies exposed to pesticides. The model also indicated that contact toxicity to foragers that resulted in reduced longevity, pesticides application method and rate, and nectar contamination (increase mortality of all lifestages) were critical parameters for colony growth and survival. The model predicted that foliar applications had the greatest effect on colony dynamics as they reduced forager lifespan and that of colony members through nectar contamination. However, simulations also revealed that the relative impact of the parameters fluctuated throughout the simulation period in relation to the status of other inputs. The time-conditional results from the simulations were expected as colony size, brood production and age structure change throughout a colony's yearly cycle.

Building on the VARROAPOP model, Becher et al. (2014) [33] constructed a model that includes effects of Varroa mites, viruses and related impaired foraging behavior, changes in landscape structure, and pesticide exposure on colony growth and survival (the BEEHAVE model). Colony structure and brood production components of BEEHAVE are similar to BEEPOP and VARROAPOP in that they are generated from difference equations that create cohorts of individuals that are aged from egg to adult. Varroa populations are initialized as phoretic mites (i.e., mites on adult bees) and reproduce in brood cells. Mites can infect developing pupae with DWV and ABPV, or be virus-free and not infect their pupal host. The foraging model is an agent-based system simulating foraging on actual and hypothetical landscapes. Cohorts of workers that reach the age of first foraging join the foraging population. The bees' foraging decisions are dependent on the energetic efficiency of a food source (for nectar collection) or by the total trip duration (for pollen collection). Energy gains and losses of each foraging trip are calculated. Forager mortality occurs after predetermined distances or maximum age is achieved.

BEEHAVE was used to simulate the effects that combinations of stress factors such as Varroa, low forage availability, and increased forager mortality, as a proxy for pesticide exposure, might have on colony survival at different times of year. Using data from large-scale colony feeding studies for model calibration and validation, Abi-Akar et al. (2020) [121] identified factors that are critical to overwintering success for a range of landscape compositions, weather conditions, and supplemental feeding schedules. The model predicted that factors affecting over wintering survival were weather conditions, and the effects that supplemental feeding in the fall have on colony size and food stores.

The BEEHAVE model also was used to simulate the sublethal effects of pesticides, and determine when those effects would manifest in colony dynamics, and if they could be mitigated [123], when those effects would manifest in colony dynamics, and if the effects could be mitigated. The sublethal effects included impaired forager orientation (simulated by doubling the distance to forage during period of exposure), impaired floral handling time (doubling handling time during a foraging trip), effects on brood care (halving the brood to in-hive bee ratio) and increased forager mortality. The greatest colony-level impact occurred due to sublethal effects of brood care. However, timing of the event was critical and most pronounced in the spring when colonies were building. Handling time effects were most pronounced in late summer and had long-lasting colony level impacts. The effects could be related to the importance of collecting sufficient pollen and nectar to meet seasonal requirements. In the spring, pollen collection drives brood production and colony size. In late summer, nectar collection is needed to insure sufficient stores for overwintering. Simulations indicated that increasing the availability of non-contaminated forage could improve colony recovery following pesticide exposure by increasing brood production (spring) or nectar stores (late summer). BEEHAVE also was used to evaluate the effects of the neonicotinoid, clothianidin, on colony growth and survival. Simulations with mortality due to pollen contamination indicated adverse effects only at high levels of exposure not seen in field measurements. No colony effects were predicted when clothianidin concentrations similar to field measurements were used.

Furthering investigations into the relationship between forage availability and colony growth, Becher et al. (2016) [34] developed a software tool to examine how scout bees explore the floral landscape. BEESCOUT uses image files imported into the program to serve as a "forage map". The location of floral (i.e., food) patches relative to the colony are calculated and used to generate probabilities of scouts finding the food patches. Different scouting strategies can be used for food detection, and probabilities of scouts finding food sources can compared. Detection probabilities of forage patches are used as input for the BEEHAVE model, to explore realistic scenarios of colony growth and death in response to different stressors. The nectar and pollen resources than can be collected are incorporated into the BEEHAVE model by linking information on patch size, location food availability, and detection probability (BEESCOUT) with daily amounts of nectar and pollen provided by the patch, sugar concentration, and handling times to collect full nectar and pollen loads as initialized in BEEHAVE. The total amount on nectar and pollen provided on each day during the flowering period is then calculated based on the patch size. Predictions from the simulations where the food source was close to the colonies fit empirical data with reasonable accuracy. However, no empirical data were available for food sources further away. Still, predictions from the model show that landscape structure and connectivity of food sources can have a strong impact on colony growth. The model can be a valuable tool for designing landscape configurations that incorporate searching behavior of bees, detection probabilities, food availability, and foraging decisions that increase patch visitation rates and enhance colony development and survival.

The BEEHAVE model was designed to simulate honey bee colony dynamics, but Becher et al. (2018) [124] modified the model to simulate the growth, behavior and survival of six bumblebee species native to the UK. The model uses an agent-based approach to generate predictions based in a mapped landscape where nectar and pollen resources can be approximated based on the model of Becher et al. (2016) [34] and can include the effects of multiple stressors on colony growth and survival. The model captures colony dynamics using individual behaviors that are driven by spatio-temporal resource availability in various landscape compositions. The model incorporates energy budgets and depletion of resources in mapped landscapes so that interspecific and intraspecific competition can emerge. Simulations on the effects of pesticide exposure indicate the impact on colony survival depends on the population size at the time of exposure.

Computational models have been developed to simulate the interactions between Varroa and the viruses they transmit, and predict the effects on colony survival. The first Varroa-virus model was developed by Martin (1998) [26] and Martin (2001) [38] and simulated the dynamics of a colony infested with Varroa that transmitted either DWV or ABPV to the parasitized bees. The model indicated that DWV would remain at low levels through much of the spring and summer with little effect on colony growth. However, DWV infections increased in the fall as there were more mites carrying the virus. This led to high levels of DWV-infected bees in the winter cluster and eventual colony death due to the reduced lifespan of infected bees and the effects on the colony age structure. In contrast, ABPV causes rapid death of infected bees, The model predicted that ABPV was difficult to establish in a colony especially at low mite populations. Colony death occurred only when Varroa populations were large with high levels of mites carrying the virus. The model accurately predicted that DWV would become more widely established in bee populations than ABPV and be the predominant cause of death in Varroa infested colonies.

Wilkinson and Smith (2002) [27] presented a difference equation model (with a time step of one day) for adult worker bees, worker bee brood, drone bee brood, phoretic mites, and non-phoretic mites. Unlike Martin (1998) [26], a simple but detailed sensitivity analysis was performed in which input variables were ranked in terms of relative importance. They found that viable reproduction in worker cells, percentage of drone brood, mite invasion of drone brood and mite reproduction in drone brood are the most important factors determining Varroa population growth rate.

Ratti et al. (2016) [46] studied the effect of swarming on a colony infested with mites and ABPV using discrete-continuous framework based on Ratti et al. (2015) [45]. They studied swarming due to overcrowding and at fixed time intervals. This study was done mainly computationally. The model predicted that in a mite-virus infested colony, the proportion of mites leaving the parent colony with the swarm may determine the fate of the colony. A colony, that otherwise dies off due to virus, can survive if the percentage of mites leaving the parent colony is above a critical value.

In addition to applying mathematical and computational models for simulation, machine learning also is being used to predict colony survival. Calovi et al. (2021) [125] used random forest (RF) matching learning method to evaluate the importance of weather, topography, land use, and management factors on the apiary and colony-level overwintering mortality. Predictions on overwintering survival from the model were more accurate at the colony than the apiary level. In colonies where Varroa was controlled, temperature (expressed in growing degree days) and precipitation the previous year were the most important predictors of colony survival over the winter. Topographic factors and landscape quality were not important predictors of overwintering survival. Beekeepers that provided data for the model may have fed their colonies protein and carbohydrate supplements during periods when flowering plants were unavailable making it difficult to quantify landscape effects. Since the important variables for predicting overwintering success were weather-related, a simpler predictive RF model was created without landscape or topographic variables. Predictions from the second model were as accurate as the full model. Because the model requires only weather conditions to generate predictions, it could be used to predict overwintering survival under projected climate change scenarios, and serve as a basis for developing decision support tools to effectively manage colonies in a changing climate.

6.   Gaps and future directions

Addressing the global crisis of honeybee losses will require an integrated approach to evaluate the combination of factors that culminate in the collapse of a colony population. Mathematical models can be powerful tools to provide insights on how various stresses contribute to colony losses. These theoretical tools are needed to identify mechanisms that can be difficult or costly to quantify in field and will be an essential component in developing strategies to improve colony health.

The most common cause of colony losses is Varroa mites and the viruses they transmit [7,87]. An effective method of control is the selection of Varroa tolerant bee stock. Modelling interactions between honeybees and Varroa could point to key combinations of traits that could enhance the ability of selected lines to reduce the growth of Varroa populations and viral transmission among bees. Increases in Varroa populations and virus levels are connected, and it might be assumed that the virus, especially DWV, is the cause of colony loss [94,126] due to effects on adult longevity and the age structure of the colony. Modeling the connection between virus transmission rates and longevity of infected workers could provide targets for the selection of bee stocks with traits that reduce viral transmission or replication. Models that incorporate spatial scales can provide landscape level information on the spread of mite and viruses among colonies in an apiary. This information could be valuable in generating apiary-level recommendations for the arrangement of hives and application of miticide treatments.

Climate change is an important factor contributing to the ongoing decline of pollinators [130]. Models that describe the effects of environmental factors on plant growth and flowering characteristics (i.e., phenology, nutrients in pollen and nectar) under current and future weather scenarios could be instrumental in predicting where crops and pollinator plantings could be established, and the traits they would require to survive. Climate change also could affect the yearly cycle of honeybee colonies especially if fall and winter temperatures are conducive for flight. Models that incorporate future weather conditions generated by climate projection models could determine the effects of late season foraging activity on the age structure of colonies, and predict the impact on overwinter survival with and without additional stress factors such as pathogens or poor nutrition.

Pesticides (fungicides and insecticides) are widely considered to be a key factor in the reduction of insect pollinators. The effects of pesticides vary depending upon modes of action. Responses to pesticides that are neurotoxins differ from those that are metabolic inhibitors or that disrupt development. As more information becomes available on sublethal effects of pesticides, the impact on colony survival during and after exposure could be predicted by incorporating the information into colony population dynamics models. Analyses to determine periods when colonies are most vulnerable to loss from sublethal effects also could be defined so that beekeepers could avoid areas where pesticides are applied during those times. Models that include effects of pathogens or poor nutrition that reduce brood rearing and adult bee longevity also could capture the interactions with pesticide exposure that culminate in colony loss. Recommendations for rescuing colonies exposed to pesticides by improving nutrition [134] or controlling pests such as Varroa might also be generated within the framework of colony dynamics.

Colony losses are economic liabilities to beekeepers that jeopardize the solvency of their businesses. Economic models are needed to support decision making in implementing management protocols, estimating returns on investment at a colony level, setting prices for colony rental for pollination and assessing the cost of colony losses at different times of year. An economic model of a commercial beekeeping operation has been reported [5], but others are needed that include different management strategies and revenue streams that could increase the profitability of beekeeping operations. Economic models also could be a tool to continually update the value of honey bee pollination in crop production and the agricultural economy as acreages of bee pollinated crops increase.

One goal of our paper is to review what we have done in mathematical modeling of honeybees and suggest future studies that could occur. We found that with few exceptions [22,29,50,56,129], the validation of theoretical models with field data is lacking. Model validation with experimental data is needed to provide confidence in model predictions. The complexities of factors that determine the survivorship of colonies will require the construction of models, but they must be based on field data that include determinants of colony health (e.g., colony size, mite loads, presence of disease). Collaborations between mathematicians and biologist are needed to construct models that put field data into a broader context and make use of their predictive power. These collaborations could be instrumental in reducing colony losses because the models would provide a deeper understanding of the shifting vulnerabilities of colonies throughout their seasonal cycle, and serve as a basis for management tools to preserve colony health.

Acknowledgments

This research is partially supported by NSF-DMS (Award Number 1716802 & 2052820); NSF- IOS/DMS (Award Number 1558127) and The James S. McDonnell Foundation 21st Century Science Initiative in Studying Complex Systems Scholar Award (UHC Scholar Award 220020472).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.   Summary for papers

B.   Model details

Messan et al. (2018): [50]

where $B$ indicates brood, $N$ indicates nurses, $F$ indicates foragers, $P$ indicates pollen storage and $V$ indicates the rate of change of VG in nurse bees.

Laomettachit et al. (2021): [71] Core model

where $B_s$ represents the number of healthy brood, $B_i$ represents the number of infected brood, $N_s$ represents the number of healthy nurses and $N_i$ represents the number of infected nurses.

Extended model

where $R_s$, $R_i$, $F_s$ and $F_i$ represent the number of nectar-receivers, infected nectar-receivers, healthy foragers, and infected foragers, respectively. Nectar-receivers and foragers can also change their nectar-loaded state between unloaded (subscript 0) and loaded (subscript 1).

Sweeney et al. (2019): [49]

Messan et al. (2017): [48]

where $\rho_{ij}$ represents the migration rate of mites from Patch $j$ to Patch $i$.

Kang et al. (2016): [40]

Dénes and A. Ibrahim (2019): [93]

where $H_s$ indicates susceptible hive bees, $F_s$ indicates susceptible forager bees, $H_m$ indicates hive bees infested by non-infectious vectors and $H_i$ indicates hive bees infested by infectious vectors.

Ratti et al. (2012): [44]

where $M_I$ indicates the number of mites that carry the virus, $H_S$ and $H_I$ indicate the number of honeybees that are virus free and the number of honeybees that are infected with the virus, respectively.

Ratti et al. (2017): [47]

Jatulan et al. (2015) : [101]

where

The parameter $B_s$ indicates the number of healthy brood in the hive before $n$ hours old, $B_o$ indicates healthy brood in the hive older than $n$ hours old, $B_i$ indicates infected brood in the hive, $C$ indicates infected cell combs, $A_s$ indicates adult bees without AFB spores and $A_i$ indicates adult bees with AFB spores.

Petric et al. (2016): [107]

where $H, H_I, F,$ and $F_I$ represents the number of uninfected hive bees, infected hive bees, uninfected foragers, and infected foragers respectively.

Muhammad and Eberl (2020): [108]

Comper and Eberl (2020): [110]