Food and environmental degradation as causative agents of honey bee colonies decline: Mathematical model approach

In this research, a new compartment model of honey bee population is developed to study the effects of gradual change of food availability as a result of environmental degradation on bee population growth and development. The model is proved to be mathematical well posed and a non-trivial equilibrium point is shown to exist and asymptotically stable under certain conditions. The model predicts a critical threshold environmental degradation rate above which the population size of bees decline and subsequently collapsed. Low environmental degradation and high food availability leads to stable bee population while the reverse leads to honey bee population collapse. Global sensitivity analysis is conducted to determine the most sensitive parameters of the model that can lead to colony collapse disorder. Numerical simulations are conducted to illustrate the results. DOI:10.46481/jnsps.2021.283


Introduction
Honey bees are social insects that live in colony and each member of the colony has its role. In each colony there is a queen, hundreds of male drones and thousands of female worker bees. Honey bees feed on flower pollen, and they also store nectar in their hives to feed young ones during winter. The life cycle of bee followed the incomplete metamorphosis starting from eggs. The eggs hatched to larvae which are sealed in the cell with wax cap by the worker bees where they developed into pupae. When fully developed, new adult worker bee emerged from the comb. It takes 21 days for a new worker bee to emerge The young worker bees are in the hives and part of their duty is to clean the hive, feed the queen, nurse the brood, convert nectar into honey and store it in the hive as food for members of the colony [6]. The old worker bees are the ones responsible for foraging but sometimes return to hive to take care of the brood when there are insufficient hive bees. Worker bees live up to 4-5 weeks in the summer and up to six month in the winter [5]. Drones are the male bees and their only duty in the colony is to mate with the queen that takes place outside the hive and die immediately. They are slightly smaller than the queen and larger than the worker bees. Honey bees play vital role in pollination of flowers in plant communities which is responsible for food production [7]. Studies have shown the key role played by honey bees in agricultural fields in seed and crop productions [8,9,10].
Since 2006 there is increase in reported number of colony collapse disorder (C.C.D) cases [11,12]. Colony collapse disorder is the sudden disappearance of honey bee colony. The most common symptom of C.C.D is the absence of worker bees from the hive, in which worker bees leave for foraging and they never returned [13]. Causes of C.C.D can be attributed to a large extent on environmental degradation factors which include excessive use of pesticides and insecticides (used to fight against honey bee pests), predators and diseases [14]. Some of these neonicotinoid family and Eutrophication of floral resources associated with anthropogenic, which have been established to have negative impact on honey bees including colony collapse [15,16,17]. Since floral diversity is required to ensure the supply of nectar and pollen to colony [18], bees travel long distance in search of food and in such case, they are exposed to dramatic temperature fluctuations which can cause stress and death [13]. Other factors of environmental degradation that contribute to C.C.D includes, radiation, climate change, global warming, geomagnetic disturbance [19,16], air pollution [20], predation by other insects [21] and electromagnetic radiation from the Sun [22].
Mathematical models play fundamental role in the study of biological phenomena and dynamics of such natural occurrences providing useful insight towards better understanding of the situation. DeGrandi-Hoffman in [23], developed a mathematical model that studied the laying rate of the queen bee. In their work, the model was divided into two components: the component that study the number of eggs laid by the queen each day and proportion of eggs that developed into drones and workers. The second component of the model track the development of eggs to adult bee. In [24], Schmickl and Crailsheim studied an age structured honey bee population model using difference equations. Their model was an extension of the one in [23]. Divison of labour and nutrients allocation in a colony is considered in their model. As a fellow up to causes of colony collapse disorders, of recent, some mathematical models were proposed, amongst which is the one by [25]. In [25], a model of honey bee colony collapse due to the contamination of forager bees in regions contaminated with pesticides was developed. The model attributed C.C.D. in terms of increased homing failure by contaminated forager bees. Furthermore, Petric et al. [26] investigated the interplay between Nosema ceranae infections and increased forager losses due to exposure to agricultural pesticides as leading factors for colony failure. On the side of diseases as causative agent of hive bee C.C.D., Dénes and Ibrahim [27] formulated a deterministic model that simultaneously modeled the spread of virus by Varroa mites among the hive bees population. They found three reproduction numbers as thresholds for global stability of four possible equilibria and later determined the key parameters that caused the collapse or survival of the bee colony. In Torres and Torres [28], a mathematical model was developed of a honey bee colony with Varroa destructor mites as parasitic disease. The model simulations revealed important parameters that determined the survival of colony. Comper and Eberl [29] formulated yet another mathematical model for a colony of honey bees infected with a parasite Nosema ceranae. Numerically, they analyzed effects of the parasite infection, food storage dynamics and implications towards colony survival and annual honey yield under treatment. On another cause of C.C.D. in honey bee colony, Khoury et al. in [30] presented a model that shows how death rate of forager bees leads to colony collapse. The model predicted a critical threshold forager death rate beneath which colonies regulate a stable population size. This model was improved in [31], by adding stored food in hive and brood compartments. It indicated how interactions of food and forager mortality affect the colony dynamics.
According to the model in [30], the quantity of hive stored food is proportional to the amount of foragers in the colony. This is true for a healthy bee colony with abundant food resources in its surrounding but for a colony in an environment where its resources are facing deleterious disturbance, the dynamics may not be the same. In our model we put into consideration pollen availability in the surrounding of honey bee colonies in modeling the food collection by foragers. It was reported in [32] that availability of pollen in agricultural landscapes is essential for the successful growth and reproduction of honey bee colonies. Field experiment also revealed that pollination of wild flowers and crops may be threatened by the widespread decline of pollinators (honey bees) [33]. In view of this, we extended the models in [30,31] by adding a compartment to represent the food collected by foragers from the surrounding environment and its subsequent degradation due to environmental factors.

Model Formulation
Our model is extension of the ones in [30,31]. We add a compartment that models the food in bee colony surrounding to investigate how gradual change in the environment will affect the dynamics of bee colony. The components of the model are divided into five compartments. These are the brood class (B) which refers to the bee in either egg, larva or pupa stage of development. Broods that developed into adult bees living in the hive are denoted by (H). Bees that transit into foraging from the hive bees formed the forager class given as (F). The food stored in the hive collected by foragers is denoted by class ( f ). The nectar and pollen in the environment constitute the food in the colony surrounding is represented by (E). Therefore, the 447 equations of the model as illustrated in the schematic diagram in Figure 1, are given below: The model satisfies the nonnegative initial conditions Equation (1) describes the rate of change of brood with time.
The parameter L represents the laying rate of eggs by the queen.
Holling type III function is used in equation (1) to show how the survival of brood depends on the amount of stored food in the hive. Mechelis-Menten function is used to show the impact of hive bees to the brood survival as it was used in [31]. The parameter b controls the rate at which recruitment into foraging decreases as stored food increase. The parameter v determine how rapidly the survival rate of brood tends to one as H increases. The parameter φ in equation (1) is the rate at which the eggs are developing into young adult bees. Equation (2) describes the rate of change of hive bees (H) with time. The hive bee population is increasing by number of adult bees that emerged pupation and is decreasing by number of young bees that are recruited into foraging. Hive bees become foragers at rate α min while α max governs the strength of the effect that low food stores has in the transition to foragers. Parameter σ is the rate at which forager bees are returning back into the hive when there are insufficient hive bees to take care of the brood. Equation (3) represents the rate of change of forager bees with time. The population of foragers is increasing by the number of hive bees that are recruited into foraging and it is decreasing by the number of foragers that transit back into the hive due to social inhibition of the bee colony and foragers that die while foraging at rate m. Equation (4) represents the rate of change of stored food in the hive with time. The quantity of stored food is increasing by the amount of food collected per forager from the surrounding. The stored food compartment is modelled in such a way that food gathering by foragers depends on availability of food in the surrounding. Saturation function is used to model the effect of food density in the colony surrounding on rate of food collection by foragers. The parameter c is the amount of food collected per forager and d is the half saturation constant. The hive stored food is decreasing by the quantity of food consumed by foragers, hive bees and broad at rates γ F , γ H and γ B , respectively. Here, we assume that the food consumed by foragers and hive bees by γ A , to be the average of γ F and γ H , as their ratio in the hive equilibrates quickly [31]. Equation (5) represents the rate of change of food in the bee colony surrounding with time. We use a Mechelis-Menten function to show the impact of foragers in the production of food through pollination. The parameter λ is the food collection rate by foragers in the surrounding while ρ is the half saturation constant of the forager in the environment. The food in colony surrounding is decreasing due to environmental degradation at rate µ.

Basic properties of the model
The basic properties of the model, which include existence, uniqueness, boundedness and stabilities of equilibria can now be studied. Since the model is representing living organisms, without loss of generality, we assume that all variables and parameters are nonnegative.
To prove the existence, uniqueness and boundedness of solution, we claim the following result:

Proof.
For existence and uniqueness of solution, the system of equations (1)- (5), can be expressed as where g is Lipschitz continuous function while . It follows from the existence and uniqueness theorem in [34] that the system of equations (1) Using Standard comparison theorem [35], we get From equation (4) we have From equation (5), Hence the solution is bounded in Ω. This ends the proof of Theorem 1.

Existence and stability of equilibrium
When the food in both hive and surroundings increase with bounds, these approaches the non-trivial equilibrium which can be obtained by setting the right hand sides of system (1)-(5) to zero. Thus, after cumbersome computations, we obtained the equilibrium as: The non-trivial equilibrium is real and positive if, λ > µ and µ < λL mρ+L . Combining these conditions, we obtained the critical or threshold environmental degradation value µ crit as The bee colony will therefore collapse if µ crit > min λ, λL mρ+L .
For the local stability of equilibrium point (7), we need to compute all the eigenvalues ω for the characteristics polynomial P(ω), with respect to the jacobian matrix J * of the system evaluated at (7). Computing the Jacobian of the model equations (1)-(5) at (7) gives,  Table 1 when all conditions are satisfied (µ < λ and µ < λL mρ+L , µ = 0.044, λ = 0.07 and λL mρ+L = 0.046.) Let 2 ) 2 , c 9 = λρE * (ρ+F * ) 2 , c 10 = λF * ρ+F * . Thus, the characteristic polynomial P(ω) is defined as where I is a 5 × 5 identity matrix. Evaluating the determinant in (9), we get Hence to show LAS for (7), we need to establish that all the five roots of P(ω) are strictly less than 0. However, due to high dimensionality of P(ω), and the algebraic complexity of parameters involved, the analysis of stability by linearization becomes cumbersome analytically. Thus Numerical approach is used to 450 show the stability of the non-trivial equilibrium point (7) in Section 3.

Numerical simulations
In this section, numerical simulations are conducted to illustrate analytical and other results. Furthermore, we present sensitivities for effects of varying certain parameters in the model as they affect the bee colony and food in the hive and environment. Figure (2) shows numerical solutions for the non-trivial equilibrium using the parameter values in Table 1. In particular, λL mρ+L = 0.046 > µ = 0.044 to satisfy the conditions of Theorem 2. Using the initial conditions as used in [31], one can observe that bees and food are asymptotically converging to the non-trivial equilibrium point (7) as stated in Theorem 2.
In Figure 3, numerical simulations of the model are displayed using three different values of µ = 0.043, 0.046 and 0.047 to show the effect of environmental degradation on honey bee colony. The simulation is run for 400 days using 33, 0000 initial total number of bees in the colony. It can be observed from figure (3) that when µ = 0.043 (mild effect), the bees survived throughout the 400 days. However, when the environment degradation parameter µ is increased to 0.046, the bees population start to decline after 200 days and the colony collapsed in about 360 days. Furthermore, when the environment degradation parameter µ is increased to 0.047, the bees population start to decline in just 100 days and the colony collapsed (C.C.D. occurred) in about 270 days. Figure 4, depicts the numerical simulations of the model using three different values of µ = 0.043, 0.044, 0.045 to show the effect of environmental degradation on food stored in bee hive and environment. The simulation is run for 500 days with 1.8 × 10 4 kg of initial total food. It can be observed from the figure that when µ = 0.043 (mild effect), the food grows in abundance. However, when the environment degradation parameter µ is increased to 0.044, the food decline drastically. Similarly, when µ is increased to 0.045, the food start to decline and exhausted after about 270 days. The fluctuations of total food observed in Figure 4 translate the seasonal variation in the environment.
In order to confirm the effect of forager death on bee colony as found in [31], in Figure 5, we illustrates the numerical simulation of the model using different forager death rate. As can be observed from the figure, when m = 0.1 and 0.2 (low death rates), the bees population increase. However, when the death rate of forager m is increased to 0.3, at the initial stage there is fluctuations in the total number of bees. Subsequently, after 100 days, the number turns to zero (C.C.D. occurred). This result is in conformity with the one presented in [31].

Global sensitivity analysis
From previous numerical simulations, we tested the variation of two parameters (µ, m) as they affect bee colony and total food in the surrounding. In order to investigate further the influence of other parameters, in this section, we conduct global sensitivity analysis on seven parameters that are associated with either foraging or food consumption to determine the extent of contributions on this matter. These parameters are c (rate of food collection per foragers), α min (Rate of transition from hive to foraging), α max ( Rate of transition from hive to foraging in the absence of food), m (foragers death rate), σ ( rate of inhibition), γ A (rate of stored food consumption by foragers and hive bees), γ B (rate of stored food consumption by the brood), λ (food collection factor by foragers from surrounding) and µ (rate at which food within the hive environment is degrading). Using the method of partial rank correlation coefficients (PRCC), as presented by [37] and successfully implemented in [38,39], we carry out the global sensitivity analysis of mentioned parameters on food stored in the colony and that in the environment. The main objective is to determine the most influential parameters that contribute to colony failure as a result of food scarcity in the hive and surrounding environment. To compute the PRCC values, we use the MatLab R2017b with ranges of parameters involved as given in Table 1, varied with ±50%, divided into 1000 sample sizes. The parameter with PRCC value far away from zero, on absolute term indicates is more influential compared to one closer to zero. The results of global sensitivity analysis are presented graphically and summarized quantitatively in Figures 6, 7 and Table 2, respectively.
In Figure 6, global sensitivity analysis of aforementioned parameters on food stored in hive is displayed. It can be observed that µ and λ are the most influential parameters followed by m and c in that order. Figure 7 presents the PRCC simulations of parameters on food in the colony surrounding. The results indicate that food resources in colony surrounding is also most sensitive to µ and λ, followed by m and c respectively.
It is worth noting that in all the results of sensitivity analysis shown in Figures 6 and 7, as summarized in Table 2, the parameter µ, which is the environment degradation rate, is the 451   overall most sensitive parameter in reciprocal manner.

Discussion and concluding remarks
In this research, we extended the models in [30,31] by incorporating a compartment of food in the surrounding bee colony environment to determine the effect of environmental degradation on bee colony collapse. Mathematical analysis shows that the model is well posed and has a non-trivial equilibrium which is determined to be locally asymptotically stable under certain conditions. The model suggests that the availability of food in the colony surrounding and also within the hive have strong influence on colony growth and development. When environmental degradation rate µ is high, the quantity of food in the surrounding and hive reduces asymptotically to zero. Consequently, the number of bees also reduce to zero at higher environmental degradation. However, at lower rates of µ, the amount of food in the hive and surrounding become abundant so also the number of bees. This agree with the assertion in [31] that in an environment with abundant of food, a healthy bee colony will continue to accumulate food into the hive until they become limited by storage space, seasonal change or dearth in forage.
Furthermore, it is observed that sudden change of environmental degradation rate has a complex effect on the bee population. When the environment degradation is low, bees will accumulate enough food into their hives and any sudden increase will not affect the bee colony since there is enough stored food. When the rate is suddenly decreased it takes a while before the colony population start to increase and so also the stored food. This has been shown by experiment in [40], when food is provided to a starved colony it takes time for the bees to return to their normal way of foraging.
It is worth reporting here, the implications and impacts of our results as par the ecological and other issues are concerned. As shown from the results, high rate of environmental degradation will leads to colony collapse due to non-availability of food in the hive and surroundings. This will directly affects food/fruits production as a results of absence of pollinators (bees). In addition, this will affects the ecosystem of the surrounding due to absence of shrubs and trees. From our results, again, we can address the colony collapse by choosing appropriate critical parameter values responsible for C.C.D. To further promote the bee colony development, from agricultural point of view, by planting abundant flowering plants with minimum proximity to the colonies so as to reduce the risk of foragers going far in search of food. From Governments sides, sponsored bills and laws should be enacted to stop or reduce the usage of chemicals in controlling pests or parasites in farming. Appropriate measures should also be taken by government to monitor the transportation of bees across various locations to avoid spread of diseases amongst the colonies. As a further work, our model can be extended by incorporating spread of diseases/pests in the colony or seasonal variation in food supply or climatic changes.