Colony collapse disorder is an increasingly common problem facing honeybees. Many worker bees (small female bees that gather pollen and nectar for the hive) suddenly disappear from an otherwise apparently healthy hive, causing a sudden, drastic decline. Lots of ideas have been suggested by conservationists; some are a mite concerned about parasites, while others drone on about beekeeping practices. But the buzz is that no definitive solution has yet been found to this crisis, which affects honey producers, the bees themselves and the plants that depend on them for pollination.
The dynamics of honeybee populations is an example of what is called a wicked problem: an intractable problem arising in a complex interacting system, with many contributing factors and no clear solution. To address this problem, first you have to factor in that honeybees are not all the same: there are different types of bee in the colony that have different roles and could respond differently to crises. Next comes the inter-dependence of bee and flowering plant populations, since these plants provide food for bees and depend on bees for pollination. Then there are lots of environmental factors to consider, which can come with further complications. For example, including insecticides in a model can be difficult because they affect not only the herbivores that people are trying to eradicate, but also the pollinators. To make matters worse, the ecosystems honeybees live in can be drastically altered by extreme weather and human activities, so you have to consider that as well $\dots{}$ The list of factors seems endless, yet potential intervention measures must be carefully considered from all these angles to prevent unexpected side effects from undermining conservation efforts.
Having a realistic model of hive dynamics could help to identify the causes of colony collapse disorder. It would give us a way of predicting what would happen to the bees under different conditions, without worrying about the ethics and feasibility of experiments. We (Isabel May Harrison, Quinn Thomazin and Ali Yasin, supervised by Tom Montenegro-Johnson, Ellen Jolley and Daniel Booth) took this on for our summer project at the University of Warwick and threw some maths at the problem to see where it could lead us, trying to make some progress addressing this mathematical conundrum. We especially wanted to focus on the effects of insecticides, since that’s a clear example where people could change their behaviour to help support bee populations. Our approach was to first consider the relationships between individual factors like the number of bees or available plants and how they influence each other’s dynamics. We then put these together to form a system of differential equations, which we tried to solve analytically, and used to run numerical simulations, in the hope that relatively simple expressions could give insights into the system.
The hive model
The starting point of our project was to make an ordinary differential equation (ODE) system for the population dynamics of a hive in isolation. We divided the bees of the hive into two groups: the “hive bees” which work in the hive and “forager bees” which leave the hive to gather food. We named the populations of each group $B_1$ and $B_2$ respectively, and the number of total bees $Z = B_1 + B_2$.
In reality, the factors that contribute to this system are practically infinite and not fully understood-but the process of modelling only requires us to consider just a few of the most important factors. Specifically, we chose eclosion (the process of bees changing from pupae to adults), death for each type of bee, and the way bees change between each type during their life cycle. With an understanding of the kinds of effects each of these factors has on each population total, we can come up with terms in a system of ODEs.
The rate of bees being born increases with the total size of the hive, equal to $\beta_iZ^2/(\kappa_i^2 + Z^2)$, where $\kappa_i$ is the rate at which the hive being too big slows the bee birth rate, and $\beta_i$ is the maximum rate at which bees can be born in season $i$. Notice that, as $Z$ grows large, the whole term approaches $\beta_i$. This is a clever behaviour of queen bees that has been observed: they limit the rate at which they lay eggs based on how many larvae the colony can support. All newly adult bees are hive bees, and become forager bees at a rate dependent on the balance of the populations in the hive. In our model, there is a standard chance that a given hive bee will convert in the next day, $\sigma_{1,i}$, and when there are too many forager bees, this rate is reduced by a second term. The loss (or gain) of hive bees due to type transfer is exactly equal to the corresponding gain (or loss) of forager bees. Finally, for each type of bee there is a constant per-bee death rate ($\eta_i$ and $\phi_i$). Since hive bees are younger than forager bees and mostly stay in the hive, the hive bee death rate, $\eta_i$, is basically zero except in the winter.
Putting all this together, we get equations for the changes in the number of hive and forager bees over time: \begin{align*} \frac{\mathrm{d}B_1}{\mathrm{d}t} &= \underbrace{\beta_i\frac{Z^2}{\kappa_i^2 + Z^2}}_{\text{eclosion}} – \underbrace{\sigma_{1,i} B_1 + \sigma_{2,i} \frac{B_2^2}{Z}}_{\text{bee type conversion}} – \underbrace{\eta_i B_1}_{\substack{\text{hive bee}\\\text{death}}}\\ \frac{\mathrm{d}B_2}{\mathrm{d}t} &= \underbrace{\sigma_{1,i} B_1 – \sigma_{2,i}\frac{B_2^2}{Z}}_{\text{bee type conversion}} – \underbrace{\phi_i B_2}_{\substack{\text{forager bee}\\\text{death}}} \end{align*}
In the notation above, any function with a subscript $i$ is seasonally dependent, but does not vary day-by-day. Before we could make any pretty plots based on these equations, we needed actual season by season numbers. We managed to find some representative values by combing through zoological papers for obscure facts like how much pollen an average honeybee can carry, but finding parameter values can easily be a huge aspect of the challenge of modelling.
This system in and of itself is interesting, and has behaviour that we expect like the bee hive population spiking in summer, and having a lower stable value during spring and autumn. Plotting based on values found from real life data, we can see a stable yearly cycle of populations.
The number of hive (blue) and forager (orange) bees over time
The effect of insecticide
Next, we wanted to explore how farming management practices such as use of insecticides might be changed to help to protect the bees. Looking at what happens to the population after a setback (like being exposed to insecticides), and identifying values of parameters where the long-term behaviour of the system changes from recovery to collapse, gives an idea of how intervention measures can be designed to improve effectiveness. We assumed that insecticides are sprayed every $T$ days, and the concentration of insecticides in the environment decreases exponentially between sprays. Bees come into contact with insecticides when they eat or forage. This is measured by the toxicity of a bee (the amount of toxins it is carrying, not its objectionable social attitude), $X$, where \begin{equation*} X \propto 2^{-t’/\tau}. \end{equation*} Here $t’ = t \text{ mod } T$ is the time elapsed since the most recent spray, and $\tau$ is the half-life of the insecticide in the environment. The constant of proportionality is different for hive and forager bees, since forager bees are likely to be exposed to more insecticides. This leads to an additional death term in the model, of the form \begin{equation*} \frac{\mathrm{d}B}{\mathrm{d}t} \propto – \frac{XB}{X+L}. \end{equation*} The death rate increases with $X$, and $L$ is the value of the toxicity that would be lethal to half of the bees.
We ran simulations varying the number of days between insecticide sprays, $T$. Since the parameters vary with season, we chose to consider the seasons separately to identify when the bees are the most vulnerable. Using summer parameters, when $T$ was at most 41 days, the average population gradually declined over a long timescale, then suddenly collapsed. Reviving the (mathematical) bees yet again and running the simulation with $T=42$, the colony survived—what a difference a day makes! More surprisingly, with spring parameters the first spray killed the entire population, and no value of $T$ could save the bees.
The number of hive (blue) and forager (orange) bees when spraying pesticides every 41 days (colony collapses)…
…and when spraying pesticides every 42 days (bees survive!)
To explain this behaviour, we assumed that insecticides are used frequently enough that the toxicity is approximately constant. The hive model simplifies to \begin{align*} \frac{\mathrm{d}B_1}{\mathrm{d}t} &= \underbrace{\beta \frac{Z^2}{\kappa^2 + Z^2}}_{\text{eclosion}} \;-\; \underbrace{\sigma_1 B_1}_{{\text{recruitment}}} \;+\; \underbrace{\sigma_2 \frac{B_2^2}{Z}}_{{\substack{\text{social}\\\text{inhibition}}}} \;-\; \underbrace{x_1(T)B_1}_{{\text{insecticides}}},\\ \frac{\mathrm{d}B_2}{\mathrm{d}t} &= \underbrace{\sigma_1 B_1}_{{\text{recruitment}}} \;\,-\,\; \underbrace{\sigma_2\frac{B_2^2}{Z}}_{{\substack{\text{social}\\\text{inhibition}}}} \;\;-\;\; \underbrace{\phi B_2}_{{\substack{\text{forager bee}\\\text{death}}}} \;\;-\;\; \underbrace{x_2(T)B_2}_{{\text{insecticides}}}, \end{align*} where $x_1(T)$ and $x_2(T)$ are time-averaged death rates due to insecticide exposure that depend on $T$. We then looked for “fixed point solutions”, points where all the derivatives are zero: this means, in theory, if the bee populations reach these values, they will stay there. If the fixed point is stable, it typically represents a long-time steady state of the system. However, if it is unstable, any small deviation from the exact fixed point will throw the equilibrium out of balance; thus these points are normally not realisable in practice. (In general, dynamical systems can be much more complicated, with periodic or chaotic behaviour for example, but luckily this approach works here.)
To simplify the algebra, it’s easier to express the equations in terms of the total population ($Z = B_1 + B_2$) and the proportion of forager bees ($W=B_2/Z$), which uniquely determine $B_1$ and $B_2$. We’re left with two quadratic equations to solve for any non-zero fixed points, \begin{equation*} \sigma_2 W^2 + (\phi + \sigma_1 + x_2)W -\sigma_1 = 0, \end{equation*} \begin{equation*} Z^2 – \frac{\beta}{(\phi+x_2-x_1)W+x_1}Z+\kappa^2 = 0. \end{equation*}
Below, there is a plot of the fixed point solutions, with stable solutions in solid lines and unstable solutions in dashed lines. As expected from the simulations, there is a non-zero stable solution in summer only when $T\geq42$, whereas the minimum value of $T$ in spring is much higher. This graph tells us that excluding other factors, for a given pesticide period $T$, if the initial bee total is below the unstable solution it will trend downwards until the bees are wiped out. If it is above the unstable solution, the population will trend upwards to the stable solution. While it tells us that the population can only be stable with $T \geq 42$ in summer, it also tells us that the greater the period, the lower the initial population would need to be for it to collapse, and the higher the population that can be sustained in the long term. (In other words: insecticides are bad news for the bees!)
To find the minimum value, $T_0$, for the hive to survive, ie the maximum frequency with which we can spray our pesticides without killing all our bees, we can take the discriminant of our second quadratic (since this will tell us whether we have any non-zero solutions), which simplifies to \begin{equation*} \frac{\beta}{\kappa} = 2[(\phi+x_2(T_0)-x_1(T_0))W(T_0)+x_1(T_0)]. \end{equation*}
Steady states in summer…
…and spring
This shows that $T_0$ depends on the ratio of biological parameters $\beta/\kappa$ (related to eclosion), whose values depend on the season. Inverting this to express $T_0$ in terms of $\beta/\kappa$ is messy but can be done graphically, using the plot below. The values in spring and summer are indicated on the plot. The minimum period increases drastically as $\beta/\kappa$ decreases, and by the time it has reached the spring value the minimum period is so large that assuming the toxicity is constant isn’t a good assumption. Of course honeybees don’t actually go extinct every spring; rather, this suggests that we should try playing with some other parameters, like the concentration of insecticides used.
Minimum period of pesticide sprays that the bees can survive
Can’t take the heat?
Summer heatwaves are no picnic for bees: the temperature of the hive needs to be regulated through fanning (workers using their wings to circulate air) and collecting water, which may be hard to find. We explored the consequences of a heatwave in summer, represented by a temporary increase in the death rate on days 131 to 134. Normally this would be catastrophic for our bee population, but given our previous results, we wondered if reducing the frequency of pesticide sprays in response could help.
In the simulation below, insecticides were used every 42 days in summer. The simulation over a two-year period shows that the bee population did not recover after the heatwave, and gradually declined to zero over the following winter and spring. Then, we ran another simulation where one insecticide spray just after the heatwave (on day 135) was skipped, and, as you can see in figure below, this time the bee population was able to recover!
While being sprayed with pesticides every 42 days, bee populations could not recover from a heat wave
Skipping a single pesticide spray just after the heatwave seemed to save the bees
Buzzing off
In this project we tried to illustrate how a wicked problem like honeybee populations can be modelled, by starting with a hive model and taking insecticides and heatwaves as examples to show how additional factors can be incorporated. The model suggests that bees are more susceptible to insecticides in spring than in summer, due to the seasonal change of the eclosion rate. With parameters specific to a given environment, the simulations could help to indicate whether insecticides are being used safely. Alternatively, intervention methods that target the eclosion rate, if feasible, might be able to improve resilience to insecticides.
It’s the nature of a wicked problem that no model can ever be constructed that would describe it in its entirety. The objective, rather, should be to capture key relationships of the system in a way that yields useful results.
No bees were harmed in the making of this model.
More from Chalkdust
In conversation with Anne Skeldon
Sam Kay and Ellen Jolley race against the body clock with the University of Surrey professor
On the cover: A new puzzle every day
Matthew Scroggs puts polyominoes inside polygons
Cheating at Tetris
Ffreuer Bristow works out which pieces can force someone to lose
The secret life of words
Barry Pawlik builds words and sentences from doubles
Coffee problems
Patrick Hedfeld has bean there, modelled that
Solving a cubic in 1646
Elinor Flavell takes squaring the circle to another level