Modeling Colony Collapse Disorder in honeybees as a contagion

Honeybee pollination accounts annually for over $14 billion in United States agriculture alone. Within the past decade there has been a mysterious mass die-off of honeybees, an estimated 10 million beehives and sometimes as much as 90% of an apiary. There is still no consensus on what causes this phenomenon, called Colony Collapse Disorder, or CCD. Several mathematical models have studied CCD by only focusing on infection dynamics. We created a model to account for both healthy hive dynamics and hive extinction due to CCD, modeling CCD via a transmissible infection brought to the hive by foragers. The system of three ordinary differential equations accounts for multiple hive population behaviors including Allee effects and colony collapse. Numerical analysis leads to critical hive sizes for multiple scenarios and highlights the role of accelerated forager recruitment in emptying hives during colony collapse. 2010 Mathematics Subject Classification. Primary: 92D30, 92D40; Secondary: 37F99.


Introduction
The honeybee industry has great economic significance.Honeybees (Apis mellifera) play a dominant role in pollination, being one of the primary managed pollinators available for field and outdoor fruit crops [17].Honeybee pollination accounts annually for $14 billion in United States agriculture alone [16].One widely-cited study estimates that pollinators altogether provide over $200 billion to the global economy [10].Much of the food production for the world is dependent on honeybees; in the United States, for example, one third of a person's diet comes from insect-pollinated plants, and honeybees are responsible for 80% of that.With expanding food production there is a further increase in the need for honeybees [12].
Since 2006, a massive loss of honeybee colonies has been reported.Some apiaries have lost up to 90% of their colonies [16].In this time the average overwinter loss of honeybee colonies in the United States has exceeded 30% consistently [26].Similar hive loss has been a concern in Canada, throughout Europe, and in Japan.The cause of this colony loss is not yet known, and it has been termed Colony Collapse Disorder, or CCD.Surveys of pathogens associated with colony collapse events have identified many disease organisms present and several newly described bee pathogens have been linked with CCD, but at the time of writing no definite single agent has been identified as the cause of CCD [3,16].The emerging consensus is that CCD is not caused by any single factor but is the result of a complex combination of multiple factors, including certain agricultural pesticides, beekeeper-applied chemicals, poor nutrition, pathogens, and parasites [24].In CCD, the adult bees vacate the collapsing hives in droves, leaving behind the queen, her brood, and frames full of honey and pollen.None of the absconding bees are found dead near the hive.Strangely, the abandoned stores of honey, which would normally be 'robbed' by neighboring bees or other organisms, remain untouched [22].
While previous mathematical models have tried to understand this phenomenon, they typically focused only on disease dynamics.To the authors' knowledge, at the time of writing there was no model that accounted for healthy dynamics alongside infection dynamics.DeGrandi-Hoffman et al. produced the first time-based model of honeybee colony growth [19].Their goal was only to include healthy hive dynamics without introducing any sort of infection.Models created after this looked mostly at infections that are well known such as Varroa mites [5] and the parasite Nosema ceranae [11].Khoury et al. established a model to study different death rates of foragers and the impact these had on colony growth and development [13,14].They then linked their results to CCD.Eberl et al. also studied a model connecting Varroa mites to CCD [9].They found the importance of thresholds for hive worker bees to maintain and take care of the brood.Sumpter and Martin studied Varroa infestation's role on viral epidemics, finding that sufficiently large mite infestations may make hives vulnerable to collapse from viral epidemics [21].
This study investigates the scenario that CCD is provoked by a transmissible pathogen or contaminant introduced by foraging activity, in order to understand how colony collapse might occur as a perturbation of normal hive dynamics.The model presented in this study therefore begins by accounting for healthy hive dynamics.Since CCD only seems to affect worker bees, as they are the ones leaving, the model only takes into account hive bee and forager bee classes.Infection dynamics, postulated on pathogen/contaminant introduction by foragers, are then incorporated using a term consistent with both direct (bee-to-bee) and indirect (via contaminated plant vectors) transmission since there is no consensus as to the precise transmission route for an infection leading to CCD.The model represents a simple framework to explore these dynamics.Classical qualitative analysis identifies the criteria for each possible outcome, while numerical analysis and parameter estimates provide context for the operative range of the model.This analysis leads to some interesting biological results in terms of the parameters of the model, in particular egg laying and maturation and forager recruitment rates.Results highlight the influence of the transition rate from hive bee to forager bee, along with the queen's reproductive rate, in determining the fate of an affected colony.

Healthy hive dynamics
The simple hive model includes only two classes: hive bees (H) and forager bees (F ), assumed uninfected.The bee life cycle sketched in the model begins with hive bees feeding and cleaning the brood so that they reach pupation and emerge as adults: eggs are laid by the queen at a constant rate L, and can therefore later mature into young hive bees at a maximum rate L, mediated by the efforts of the hive worker class.As the hive bees feed the brood and clean the hive, the emergence rate is multiplied by a saturation function H H+Ω , where Ω is the number of hive bees needed for the emergence rate to reach 1/2 L. This fraction gives the proportion of eggs which survive to eclosion, given as an increasing function of the number of hive workers available to feed and protect the brood.Thus the eclosion rate is given by L H H+Ω .(Khoury et al. [14,15] use a similar saturation form to model eclosion but consider it instead as a function of the entire hive population, not just the hive bees.We consider here the hive bees' role to be the critical one in fostering emergence, since food shortages prompt hive bees to forage, as discussed next.) As their brains and wings develop, hive bees become forager bees, able to bring food to the hive; this occurs naturally after about 3 weeks (represented in the model by a rate γ), but the process can be accelerated if the forager class is depleted and the hive needs more food.The maximum additional maturation rate is given by α, multiplied by a factor Φ F +Φ which measures the hive's need for more foragers, where Φ is the number of foragers at which the additional maturation rate is α/2.(In comparison, Khoury et al. [14,15] consider instead a maximum forager recruitment rate reduced for large forager classes via social inhibition.)Finally, forager bees die from natural mortality at a rate µ 1 .(All the given parameters are taken as constant over the span of a season, the timescale on which colony collapse is observed.)Together, these assumptions lead to the system

Infection dynamics
Since there is no consensus yet on the means by which individual bees become infected, we construct a model for a hypothetical infection consistent with either direct bee-to-bee transmission or vector transmission via contaminated flowers.We first consider that [healthy] forager bees F become infected via contact with contaminated plants, and that plants become contaminated via contact with infected bees.This assumption agrees with the observed dynamics of CCD as foragers are the ones most affected by the phenomenon.We further assume that forager bees interact very little with hive bees while foraging, with only the exchange of collected materials happening.If we further assume that the rate of bee-plant contacts is saturated in plants but not in bees (i.e., flowers are not continuously occupied by bees during daylight hours, and if the flower density increases, bees cannot visit more flowers per day than they are already visiting), then the populations of infected bees I and contaminated plant vectors V evolve according to the equations where β is the rate of bee infection, µ 2 is the per capita mortality rate of infected bees (µ 2 ≥ µ 1 ), c is the rate at which plants clear contamination (e.g., due to rain and wind), all in units of 1/time, and β is the rate of flower contamination (in units of plants per bee per day).Plant density P is assumed constant.
To simplify the model, we use a timescale argument to eliminate the plant vectors, claiming that plant contamination occurs on a faster timescale than the bee epidemic in a hive, and so we allow dV /dt → 0 and observe the resulting equilibrium value, This allows us to eliminate V , replacing V /P in the system by I/(I + K).In this way we obtain a simpler model which can also be used under the alternative assumption that bees become infected by direct contact with each other (in interactions such as communicating about feeding site locations, rather than indirectly via feeding sites acting as vectors), at a rate which saturates as the number of infected bees increases, with K as the half-saturation constant.Thus, finally, the model is where Ω, Φ, K are saturation constants measured in bees and µ 2 ≥ µ 1 .

Equilibria
The equilibria of system (1)-( 3) can be identified and analyzed via standard methods, using the Jacobian matrix to determine conditions for local asymptotic stability.Details are presented in the appendix, and results summarized here.The system has up to five biologically meaningful equilibria: hive extinction, up to two disease-free equilibria, and up to two endemic equilibria.
The second of each pair, if it exists, is unstable, so this corresponds to three outcomes for the hive: colony collapse, a healthy hive, or survival in an endemic state.Equilibrium values and the conditions for [biologically meaningful] existence and for local asymptotic stability (LAS) are given in Table 1.For convenience, the state variables are rescaled: Existence and stability criteria are likewise given in terms of the following notation: demographic parameters and infection-related parameters (a, b, c, d have units of bees/time; r and k are dimensionless.)Of these, a is the only one not automatically positive, but by inspection of (1), if a ≤ 0 then dH/dt < 0 and the hive dies, so   This same analysis (using a next-generation operator method [7,25]) yields the infection's basic reproductive number, R 0 = βF * /µ 2 K, where F * is the equilibrium value from the disease-free equilibrium (DFE1 in Table 1).A consequence of this expression is that R 0 is only defined when this equilibrium exists (is nonnegative).
The juxtaposition of the criteria given in Table 1 leads to some interesting dynamics, including regions of bistability and even tristability in parameter space.Existence and stability of diseasefree equilibria depend on the demographical quantities a (a measure of the maximum eclosion rate relative to the baseline maturation rate, signaling the health of the hive worker class), b (a measure of the baseline forager death rate) and c (a measure of the accelerated maturation rate when the forager class is depleted).Figure 1 summarizes these criteria in terms of the a-b parameter plane, while Figure 2 illustrates them in terms of a-b-c parameter space.Basic hive persistence occurs for a > c, at the boundary of which a transcritical bifurcation occurs (shifting stability from hive extinction to hive persistence); however, a saddle-node bifurcation to the left of a = c creates a region where an Allee effect operates, and the hive only survives above a threshold critical size.
The endemic equilibria are likewise introduced via a line of transcritical bifurcations, with an adjacent curve of saddle-node bifurcations, determined as functions not only of a, b and c, but of the epidemiological measures d (baseline infection rate), r (a reciprocal of R 0 , r = f * 1 /R 0 with f *   2): XE stable in shaded region, DFE1 stable in vertically dashed region, EE1 stable in horizontally dashed region (CCD labels the region of colony collapse due to CCD) three new regions in the positive quadrant of (a, b) space (cf. Figure 3): one in which the hive survives in an endemic state, one in which an Allee effect separates the endemic state (EE1) from hive extinction, and a third in which the hive goes extinct because of the infection.This last region (marked ** in Figure 3) is classic colony collapse.The placement of the additional saddle-node bifurcation curve creating the second endemic equilibrium-and, in particular, at what point it meets the line c = (a − br)(r + 1)-depends on the value of c relative to d, r and k.As shown in the appendix, this point has coordinates endemic Allee extends into colony collapse region and thus is in the positive quadrant iff c > dr(r + 1) 2 /(r + k).In general, there are four cases, determined by where b 0 falls with regard to 0, c/(r + 1) 2 , and c/(r + 1), described in Table 2. Figure 4 illustrates the last and most complex of these (with the endemic saddle-node bifurcation curve crossing the horizontal axis between the disease-free saddle-node bifurcation curve and the endemic transcritical bifurcation line, and continuing up to meet the line in the region where a > c).However, the constraint b < dr/k prevents some of the resulting behaviors from occurring, since the two saddle-node bifurcation curves intersect at b = dr/k (Cases 3 and 4, see Figure 4).For b 0 > 0, in the interval 0 < b < min(b 0 , dr/k) the endemic saddle-node bifurcation extends the region of the endemic Allee effect into the CCD region (where R 0 > 1).For such parameter values, the accelerated forager recruitment is strong enough to sustain an endemic population when the colony would collapse under normal forager recruitment.More complex behavior is possible for b > dr/k, although this is of no biological interest since it requires (µ 1 > µ 2 ) that infected bees remain in the hive more, not less, time than healthy foragers.This behavior extends the endemic Allee effect further, into the natural extinction region (where R 0 is undefined), only possible because infection prolongs forager lifetimes when µ 1 > µ 2 .The same phenomenon is at work in the other effects (as b 0 increases to c r+1 and beyond): first a region of tristability, in which extinction, healthy hive survival, and endemic survival are all possible outcomes, and then a region exhibiting classic "backward bifurcation" effects, in which a large enough initial prevalence enables the infection to persist but hive survival is guaranteed regardless.
The boundaries of regions in Figures 3 and 4 correspond to bifurcations; for an additional perspective, the bifurcation diagrams in Figure 5 illustrate the five different possible combinations of these bifurcations as a increases, for fixed values of b, c, d, r, k (corresponding to horizontal "slices" of Figures 3 and 4), under the constraint b < dr/k.
In summary, the model exhibits several different behaviors depending on demographic (a, b, c) and infection-related (d, r, k) parameters, including a healthy hive, an endemic state, and hive extinction.Of most interest is the region in parameter space representing colony collapse due to the infection.The following section develops some parameter estimates in order to study specific scenarios which might cause this collapse.

Parameter estimation and quantitative analysis
Some model parameters were taken or estimated from the literature.Bodenheimer estimated an egg-laying rate of 1000-2000 eggs/day [2], so we take the average of the two extremes (1500) for L. A honeybee's brain [27] and wings typically mature, enabling it to forage, in about 3 weeks (1/γ = 21 days), although the period can be shortened to as little as 1 week (7 days) [1].This maximum rate γ + α then leads to a value of α = 1 7 − 1 21 = 2 21 days −1 .Finally, a typical forager survives about 3 weeks (1/µ 1 = 21 days) before its wings wear out, stranding it [1].
Ranges for the remaining parameters were estimated more heuristically.For instance, the death or departure rate for infected forager bees µ 2 should certainly exceed the mortality rate µ 1 of uninfected foragers, but in fact infected foragers are likely to develop symptoms of infection much sooner after contamination than that.For numerical purposes we estimate that an infected bee discovers/manifests its illness within a day of contagion, and therefore suppose µ 2 ≈ 1/day.Likewise in estimating the rate β of potentially infectious contacts made by foragers, whether with contaminated flowers or other infected bees, the primary exposure is in the foraging process, where bees come into contact with up to thousands of flowers per day (see estimation for K below) and potentially thousands of other foragers from the hive.At a bare minimum, therefore, we expect β ≥ 1/day.
For the half-saturation constants we consider typical hive dynamics to come up with a range or upper bound for possible values.For Ω, we consider that the hive must have a certain number of hive bees in order for eclosion to occur.Mature hive populations vary from 15,000-25,000 in the spring, at the start of foraging season, to 50,000-60,000 at summer's end.If we assume that a hive of 50,000 bees has an eclosion (brood emergence rate) of at least 90%, this makes H H+Ω ≥ 0.9 with H = 25, 000 bees if we assume the hive population is split evenly between hive and forager bees (assuming the drone population is comparatively very small), leading to an upper bound of Ω ≤ 2778 bees.Assuming 90% eclosion is reached for smaller hives than 50,000 bees leads to an even lower estimate for Ω.For Φ we similarly assume that a mature hive of 50,000 bees uses no more than 10% additional forager recruitment, so that under a 50/50 hive worker/forager split, F = 25, 000 bees makes Φ Φ+F ≤ 1  10 .This again leads to an upper bound of Φ ≤ 2778 bees; if additional forager recruitment tapers off faster, the bound for Φ is even lower.For simplicity we will round both bounds to 3000 bees.
For K we apply data from flowers to the plant-based definition K = P c/ β given in developing the model, with plant density P , flower contamination rate β, and plant contamination clearance rate c.California almonds are estimated to bloom at a density of about 2 million blooms per acre [23], and two hives are commonly used per acre, making each hive responsible for pollinating 1 million flowers (so let P = 10 6 plants).Foragers are commonly cited to make around thirty trips a day and can visit up to 100 flowers in a single trip, making β ≈ 3000 plants/bee/day.We assume the plants take about a week to clear any contamination, making c =1 7 /day.Using these data we get an estimate for K of about 50 bees.
Parameter estimates and ranges are summarized in Table 3.
These estimates produce values of a ≥ 1357, b ≤ 143, and c ≤ 286, all unaffected by the gross uncertainty in infection-related parameters.Although the half-saturation constants Ω and Φ also have considerable uncertainty, reducing them from their upper bound of 3000 bees/day preserves the ordering b < c < a as long as Φ ≤ 2Ω.Also, given that K << Φ and µ 2 ≤ β, r = µ 2 K/βΦ << 1, so that not only is b < c but also b < c/(r + 1)2 .The significance of these inequalities can be seen by returning to Figure 3. a > c guarantees hive survival in some form; b < c/(r + 1) 2 ensures that the infection will not die out (R 0 > 1 when the hive survives), and for lower values of a corresponds to the region where colony collapse occurs (marked ** in the figure).The two inequalities together place the parameter set in the region where a sole endemic equilibrium is the global attractor.
The parameter estimates are therefore consistent with CCD save that the egg production rate is about 3-4 times as high as in the CCD region.Since widespread CCD has been observed, one possible inference is a missing feedback mechanism in the model, perhaps the dependence of the queen on the workers' food production, which is obviously severely curtailed when a large proportion of foragers are infected or gone.
As one illustration of how the model accounts for colony collapse, we consider a scenario in which the egg-laying/maximum eclosion rate L has been reduced to 630 bees/day.Figure 6 illustrates the dynamics that result when a hive operating at the disease-free equilibrium (stable in the (H, F )subspace where I = 0) sees one of its foragers infected at time t = 50 days  3, except Ω = Φ = 1500 bees).Hive workers in dark gray, healthy foragers in light gray, infected foragers in black.
quickly exhaust first the forager class and then the hive worker class, until it is no longer possible to sustain the hive and the colony collapses.In a second, related scenario, L remains at its original value and one infected forager is introduced at time 0, leading to an endemic state, until at time 50 the rate L is reduced to 210 bees/day, shifting the dynamics to a region in which the infection causes colony collapse (Figure 7; note healthy and infected forager curves are nearly equal).Both scenarios highlight the vulnerability to CCD of hives operating under an Allee effect.

Discussion
In order to understand how Colony Collapse Disorder affects normal honeybee population dynamics, we constructed a simple model tracking hive bees and forager bees, with nonlinear transition rates reflecting the influence of the worker classes on these processes, and introducing CCD via a hypothetical contagion spread to foragers.Saturation in the eclosion, additional forager recruitment, and infection rates creates a rich bifurcation structure which yields complex behaviors, including Allee effects and so-called backward bifurcations.As many as three possible outcomes can occur: healthy hive survival, endemic persistence, and colony collapse, with one specific region (marked ** in Figure 3) representing CCD: extinction due to infection.This extends previous results such as those of [9,14] to include CCD through an explicit infection.
The primary condition for ensuring hive survival independently of initial hive size (a > c) is that the egg-laying rate exceed the maximum forager recruitment rate at the egg-laying half-saturation threshold, a measure of the robustness of hive dynamics independent of the effects of CCD.If this criterion fails to hold, then extinction is always possible; if, however, bees live long enough on average to become foragers under accelerated recruitment (b < c), then hive survival may still be possible under an Allee effect (through a backward bifurcation in the healthy hive dynamics), where a large enough worker population can sustain the hive.
Infection dynamics add complexity to the hive dynamics.On a base level, when R 0 > 1 any hive survival becomes an endemic state.However, CCD also affects the hive's ability to survive, numerous factors including the proportions of hive workers in the hive and of eclosion saturation could explain this.
draining it to the point of collapse in circumstances when a large enough healthy hive would normally survive (long-lived bees but a borderline egg production rate).This is classic colony collapse.If the additional hive-to-forager recruitment (as measured by c) is great enough, the range of parameter values which cause colony collapse is reduced but not eliminated (via the endemic saddle-node bifurcation).This peculiar outcome derives from CCD's depletion of the forager class spurring greater recruitment into the forager class (hive survival normally depends on the worker population).
Although the rescaled model exhibits multistability which makes all 7 possible combinations of the three outcomes (healthy hive, endemic state, and extinction) possible for some parameter values near the boundaries of the simple survival/extinction and healthy/infected hive thresholdsmost strikingly superimposing the Allee effect and the infection-related backward bifurcation in a region of tristability, where initial hive size and outbreak magnitude determine whether the hive survives as well as whether the infection persists-the constraint on relative departure rates of infected foragers allows only an Allee effect and classic R 0 -based infection persistence, apart from the colony collapse outcome.The model's primary limitations are the focus on hive bees' role in fostering eclosion and the foragers' role in spreading infection through encounters at feeding sites outside the hive and in collaborative efforts to exploit them.
Our rough numerical estimates place a typical commercial hive in a scenario where the infection persists but should not drive the colony to collapse; however, a weakened egg laying/maximum eclosion rate, as may happen if a depleted forager class fails to bring in sufficient food, would bring the hive into CCD territory.Future work may entail incorporation of this dependence of eclosion on the forager population; Khoury et al. recently studied possible mechanisms for modeling food in this context [15].In the meantime, the present model allows us not only to identify through R 0 what criteria enable CCD-promoting infections to invade a hive, but the critical thresholds for a hive to survive such an invasion: namely, the balance between egg production and maximal forager recruitment (a demographic tug-of-war on the hive worker class).The depletion of the hive via infection accelerating forager loss, which in turn accelerates hive worker recruitment to replace lost foragers, is consistent with the conclusions of [14] but goes further to illustrate explicitly how contagion can lead to rapid colony collapse.Although CCD is modelled here via a generic contagion, and no control methods are explicitly incorporated, the association of low eclosion and even lower forager mortality rates, relative to emergency hive-to-forager recruitment (b << a << c), with colony collapse (cf.Figures 3 and 4) highlights the vulnerability of hives which rely on long-lived foragers to maintain the hive worker-forager balance when eclosion is slow, since a CCD-related infection outbreak can upset that balance by draining the forager class below the Allee threshold.In such a case, only boosting eclosion can avert the hive's vulnerability to CCD.
Future work may also consider contagion of CCD-related infection agents (whether parasites, viruses, or contaminants) among multiple hives with shared food sources, as is common in commercial uses of honeybees as pollinators.
• Taking (ii) and (iv) in ( 2 Local stability analysis involves the system's Jacobian matrix, For the XE this simplifies to since µ 1 , µ 2 > 0, the XE is thus locally asymptotically stable (LAS) iff L < (α+γ)Ω; that is, a small hive dies out when the maximum birth rate of new bees is less than the maximum hive-to-forager recruitment rate at the hatching half-saturation constant.

A Disease-free equilibria
We have equilibrium conditions (i), (iv) and (v).Substituting (i) I * = 0 in (iv) and solving for H * yields which we in turn substitute into (v): Multiplying through by (F * + Φ) and moving everything to one side allows us to rewrite this condition as a quadratic in F * /Φ, f (F * /Φ) = 0, where  When c < a (unique DFE), the criterion (F * /Φ)+1 > c/a to make H * > in ( 6) is automatically satisfied.When instead b < a < c < (a+b) 2 /4b, the criterion can be shown to be satisfied as follows: shows the criterion holds for the smaller F * , which implies that it must hold for the larger F * as well.
The step where square roots are taken in the above progression is justified To determine stability of DFEs, we consider the Jacobian The last entry gives one eigenvalue, which is negative iff βF * /µ 2 K < 1.This criterion addresses disease dynamics, and a quick calculation using a next-generation operator method ( [7] or [25]) yields R 0 = βF * /µ 2 K, so this is the condition R 0 < 1 (a DFE can only be LAS if the disease dynamics are weak).The remaining 2 × 2 submatrix involving the hive population dynamics decouples from the disease dynamics and can be analyzed using the 2-D Routh-Hurwitz criteria (negative trace and positive determinant).The trace of the submatrix can be seen to be negative if one uses (iv) to observe that the upper left entry is − LH * (H * +Ω) 2 < 0. Since the entry in row 2, column 1 can similarly be rewritten L H * +Ω using (iv), the determinant of the submatrix is now Thus the stability condition that the determinant be positive simplifies to In the case where the DFE is unique, we have c < a, so that (a + b) 2 − 4bc > (a − b) 2 ; thus the unique positive F * has satisfying the condition.
In the case where two DFEs exist, we have c < (a + b) 2 /4b, so that the larger F * has again satisfying the condition.However, the smaller F * can be shown not to obey the condition, applying the triangle inequality Thus the smaller DFE is unstable when it exists, whereas the larger DFE is LAS (when it exists) iff R 0 < 1.

B Endemic equilibria
We have equilibrium conditions (ii), (iv) and (vi).Rescaling h * = H * /Ω, f * = F * /Φ and i * = I * /K and defining d = βΦ, r = µ 2 K/βΦ and k = K/Φ, we rewrite (ii) as f * = r(i * + 1), from which Note that if a ≤ br then B ′ , C ′ > 0 and g has no real roots.We therefore assume otherwise in the analysis that follows.We also observe that each positive root of g corresponds to a biologically meaningful EE (i.e., all components positive) by noting that, for i * > 0, from (ii) f * = r(i * +1) > 0, and then from (vi (0 < C ′ < (B ′ ) 2 /4A ′ ), and no real roots otherwise.We also note that the above interval always exists (i.e., the lower bound is not above the upper bound) as long as a > br, by algebra similar to that used in analyzing the DFE(s).3 and 4 there is a surprising region of XE/DFE/EE tristability, and in case 4 there is even a region of DFE/EE bistability.
However, the constraint b < dr/k precludes the behaviors unique to Cases 3 and 4, since some algebra shows that b = dr/k precisely where θ 0 = θ 3 , which for Cases 3 and 4 occurs below c (r+1) 2 .Thus for Cases 3 and 4, the behavior of the original (unrescaled) system is limited to the region 0 < b < dr/k which features colony collapse but never allows DFE1 and EE1 to be stable simultaneously, nor does it allow EE1/EE2 to exist when DFE1/DFE2 do not exist as well-both of those outcomes would require µ 2 < µ 1 , infected forager bees remaining longer in the hive than their uninfected counterparts.Thus the only five bifurcation diagrams possible (in a) are those depicted in Figure 5.

Figure 1 :
Figure 1: Conditions for existence of DFEs in the a-b plane for fixed c: where no DFEs exist, the hive perishes; where one DFE exists, the hive persists; and where two exist, an Allee effect operates (the second DFE being unstable)

Figure 2 :
Figure 2: Conditions for existence of DFEs in a-b-c space: no DFEs above the curved surface, one below the plane, and two between the plane and the curved surface

Figure 5 :
Figure 5: Bifurcation diagrams showing F * , I * as functions of a. Solid curves denote stable equilibria; dashed curves illustrate unstable ones.The extinction equilibrium is in black, disease-free equilibria are medium-dark gray (DFE1) and medium-light gray (DFE2), and endemic equilibria are light gray (EE1) and dark gray (EE2).The two upper graphs illustrate colony collapse due to CCD.

2 −
convenient to introduce the following notation to simplify: a = L − γΩ, b = µ 1 Φ, c = αΩ.Then B = (b − a)/b, C = (c − a)/b, and the condition that the expression for H * in (6) be positive is (F * /Φ) + 1 > c/a.Note that b, c > 0, and that a > 0 in all interesting cases, since if a ≤ 0 the hive will die out (as dH/dt < 0).The roots of f are (−B ± √ B 4C)/2.If C < 0 (c < a), the discriminant dominates |B|, and f has one positive and one negative root.If C = 0 (c = a), f has one zero root and one root of sign opposite B. If 0 < C < B 2 /4 (a < c < (a + b) 2 /4b), |B| dominates the discriminant, and f has two roots of sign opposite B. (Note that a ≤ (a + b) 2 /4b ⇔ 4ab ≤ a 2 + 2ab + b 2 ⇔ 0 ≤ a 2 − 2ab + b 2 = (a − b) 2 , so this interval always exists.)If C = B 2 /4 (c = (a + b) 2 /4b), f has one root of sign opposite B. If C > B 2 /4 (c > (a + b) 2 /4b), f has no real roots.The condition B < 0 that when f has a root, f has a positive root, is a > b.Thus f has one positive root when c < a and two positive roots when b < a < c < (a + b) 2 /4b.(Note in both cases a positive root requires a > 0.) The boundary curve (a + b) 2 = 4bc crosses the line a = b in the a-b plane at the origin and the point (c, c), so for any fixed value of c the existence of DFEs in the first quadrant is given as in Figure 1: a single DFE to the right of the vertical line a = c and a pair of DFEs below the boundary curve between the origin and the line a = c.Note that in the latter case b < c.

Figure 2
gives a three-dimensional sketch of these regions in a-b-c parameter space.

Figure 8 :
Figure 8: Existence of two endemic equilibria in the a-b plane for fixed c, r, d, k: (a) top left, case 1, where the region 2EE (θ 3 < a < θ 1 for 0 ≤ b < b 0 ) is empty; (b) top right, case 2, where the region 2EE fits entirely inside the region 2DFE (θ 0 < a < θ 1 for 0 ≤ b < c (r+1) 2 ); (c) bottom left, case 3, where 2EE extends beyond 2DFE but remains within a < c; (d) bottom right, case 4, where 2EE extends beyond the line a = c.All terms are as defined in the main text.

Table 1 :
Equilibria of system (1)-(3) and criteria for their existence and stability (LAS).Stability criteria implicitly assume existence criteria.All criteria are proven in the appendix, except stability of endemic equilibria which was verified numerically.

Table 2 :
Cases for where the two endemic equilibrium bifurcation curves meet, determining the maximum a and b values for which EE2 can exist.c, d, r, k, b 0 are as defined in the main text.

Table 3 :
Parameter definitions, estimates, and units 1.The infection dynamics