Web Ecol., 22, 7–19, 2022
https://doi.org/10.5194/we-22-7-2022
Web Ecol., 22, 7–19, 2022
https://doi.org/10.5194/we-22-7-2022

Standard article 04 Mar 2022

Standard article | 04 Mar 2022 # An approach to the modeling of honey bee colonies

An approach to the modeling of honey bee colonies
Jhoana P. Romero-Leiton1, Alejandro Gutierrez2, Ivan Felipe Benavides3,5, Oscar E. Molina4, and Alejandra Pulgarín2 Jhoana P. Romero-Leiton et al.
• 2Facultad de Ciencias de la Educación, Universidad del Quindío, Armenia, Colombia
• 3Instituto de Estudios del Pacífico, Universidad Nacional, Pasto, Colombia
• 5Datambiente, Colombia

Correspondence: Jhoana P. Romero-Leiton (jpatirom3@gmail.com)

Abstract

In this work, populations of adult and immature honeybees and their honey production are studied through mathematical and statistical modeling approaches. Those models are complementary and are presented in disjunct form. They were used to show different modeling methods for honey bee population dynamics. The statistical approach consisted of a generalized linear model using data from the Department of Agriculture of the United States of America (USDA), which showed that the relationship between the number of colonies and the rate of honey production is not constant in time but decrease over the years. These models showed that when a bee population is subjected to a stress factor (i.e., habitat destruction, Varroa mite, climate variability, season, neonicotinoids, among others), the abundance of individuals decreases over time as well as the honey produced by the colonies. Finally, the mathematical approach consisted of two models: (1) a smooth model, in which conditions of existence and stability of the equilibrium solutions are determined by an ecological threshold value, and (2) a non-smooth model where the mortality rate of bees is included as a function of the number of adult bees in the population.

Share
1 Introduction

Bees are one of the most ecologically and commercially important insects in the world . There are approximately 20 000 species of bees , from which Apis mellifera is the most common for commercial purposes . Bees are of great importance not only for humans, but also for all plant species that they pollinate. The most important crops around the world are visited and pollinated by bees and rely on them for reproduction . In tropical crops, 70 % of the 1330 cultivated species are favored by these pollinators, while in European crops, 84 % of the 264 cultivated species depend on the animal pollination process, and it is known that the yields of some fruit, seed and nut crops decrease by more than 90 % without these pollinators .

Bees also contribute to the maintenance of native ecosystems, because pollination helps to regenerate trees, which help to conserve forest biodiversity, as well as preserving other ecosystem services . This makes bees important regulators of food production, forest balance and micro-climate dynamics (Winfree2010). From all species of bees, only a few are eusocial; Apis mellifera is an example of eusocial behavior . This species forms colonies so that the survival, reproduction and honey production depend on the structure and size of the colonies .

Honey is a sweet natural substance produced by honeybees using mainly the nectar of flowers, which is transformed by a group of enzymes present in the saliva of the worker bees, which is also airy and evaporates by its fluttering, and is finally stored inside the hives of the nest. Honey from Apis mellifera is one of the most important zoo-agricultural goods for commercial trade in the world . The utility of honey is diverse. For example, it works as an anti-flu agent . Honey also serves as a natural source of antioxidants which are effective in reducing the risk of heart and immune system diseases, cataracts and various inflammatory symptoms (Vit2004). Regarding the honey trade, the United States of America is the world leader in imports. Honey imports increased at a rate of 6956 t yr−1, making this country the main honey importer in the world (García2018). However, regarding production, China is the world leader (466.6 kt), followed by Turkey (94.7 kt), Iran (74.6 kt), Ukraine (73.7 kt) and the Russian Federation (68.4 kt) (Güngör2018).

Despite the great importance of bees in the world, an alarming concern has been raised in the last 20 years because of the decrease in the number of honey bees worldwide, which affects the quality of life of all human populations that depend (directly or indirectly) on them .

For these reasons, it is necessary to develop studies on honey bee population dynamics aiming to conserve and sustainably manage this species. Research in the area of mathematics and statistics has contributed to the development of models that allow us to understand the population dynamics of bee colonies . created a simple mathematical bee population model (HoPoMo) using differential equations to predict the population dynamics and resource fluctuations of a bee colony. Their model emphasized the importance of pollen supply and presented a mathematical model considering pollen and nectar as necessary foods for the colony. This model was developed based on the differential equations proposed by . built one of the most detailed population models of honey bee colony dynamics, consisting of 60 equations that track the life of a bee from an egg to an adult bee every day. The model considered the effect of seasonal changes in egg-laying rate, nurse bees on larval survival and pollen shortage on cannibalization. developed a model that focused on how internal demographic processes within a colony interact with food availability and brood rearing to alter growth in forager population size. The model was implemented as a series of differential equations based on the rate equations from the analytical models proposed by .

In this work, we describe from mathematical and statistical perspectives the relationship between the population size of honey bees (Apis mellifera) and honey production when the bee colony is subjected to some stress factor (habitat destruction, climatic variability, neonicotinoids, etc.) that externally causes the death of individuals and consequently the possible decrease in honey production. We present three different approaches. Two of them are purely mathematical and consider the interaction between immature and adult bees and the amount of honey produced using ordinary differential equations (ODEs).

First, we consider a smooth model in which a stability analysis of the equilibrium solutions is performed. The stability results are validated with numerical experiments using data obtained from the literature. Next, we consider a non-smooth model which is analyzed using Filippov's systems theory . Finally, a statistical model is developed considering honey production as the response variable, the number of colonies as a the main effect predictor, and year, latitude and area as covariates, using a generalized linear model (GLM).

2 The statistical model approach

For this section, we used the historical data of bee farming in the United States of America (USA) from 1985 to 2019, from the US Department of Agriculture (USDA), which are publicly available on USDA web page (https://usda.library.cornell.edu/concern/publications/hd76s004z?locale=en&page=3#release-items, last access: August 2022). Additionally, for a detailed revision, the R codes and the cleaned data set are available at the following GitHub link: https://github.com/jpatirom3/Honey-bees-modeling- (last access: March 2022). These data sets summarize the statistics of honey production and number of bee colonies by state and year in the USA. We used these data sets because they offer the most complete, comprehensive and publicly available data on honey production in the world. After checking the data sets for missing data, we decided to use the data from 39 states, which had complete values for all years.

A first exploration of these data sets revealed a significant decreasing trend in honey production over the years (Fig. 1), which is an expected pattern due to the worldwide decrease of working honeybees during the last 20 years . Since this kind of temporal trend can inflate the relationship between responses and predictors in statistical models , we fitted a model for each of seven 5-year sub-series (1985–1989, 1990–1994, 1995–1999, 2000–2004, 2005–2009, 2010–2014 and 2015–2019), in such a way that the temporal trend within each sub-series was significantly reduced in comparison to the full 1986–2019 series. Thus, every model was able to be built on n=194 data points, which is a reliable sample size for parameter estimation. This approach also allowed us to evaluate the temporal change of the parameters over time, which is important due to the already mentioned decline of bees. Figure 1Evolution of total honey production in the USA from 1985 to 2019. Units in the y axis are shown as thousand pounds. The solid line represents a temporal trend fitted by a locally weighted polynomial regression.

Thus, our modeling approach considered honey production (HP) (pounds) as the response variable for each of the seven 5-year sub-series; number of colonies (NC) as the main effect predictor; and year (Y), latitude (LAT) and area (AR) (hectares) of each state as covariates. Since NC and HP exhibited highly skewed distributions, we used an additive GLM log–log with a logged response and predictor in order to achieve normal and homoscedastic residuals . Equation (1) shows the specification of this model, where e is an error term with Gaussian distribution and Id is an identity link function.

$\begin{array}{}\text{(1)}& \mathrm{Id}\left(\mathrm{logHP}{\right)}_{i=\mathrm{1}}^{\mathrm{7}}={\mathit{\beta }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{2}}\left(\mathrm{logNC}\right)+{\mathit{\beta }}_{\mathrm{3}}Y+{\mathit{\beta }}_{\mathrm{4}}\mathrm{LAT}+{\mathit{\beta }}_{\mathrm{5}}\mathrm{AR}+e\end{array}$

Estimation of β parameters was performed trough the iterative reweighted least squares method (IWLS) (Stirling1984). Different combinations of covariates were fitted in a top-down stepwise modeling approach, which started by fitting a saturated model with all covariates and ended up with a reduced final model with a minimum number of covariates that captured most of the variance in the response variable (Table 1). The Akaike information criterion (AIC) and the root mean square error (RMSE) were used as metrics of performance and prediction accuracy to compare among models and define the final models for each year's sub-series. β2, of main interest for this section, was interpreted as a percentage unit change in HP caused by a percentage unit change in NC, while keeping the other covariates parameters constant. All these procedures were done in R studio (R Team2019).

Table 1Results of analysis of deviance for the statistical models fitted to the data for every 5-year sub-series. Residual deviance represents the degree in which the response variable can be predicted by model predictors. The higher the value, the better the model is able to make predictions. The F value is the rate of model variance; it provides a way to compare the fitted model against a model with only intercept and no predictors. Higher F values represent better fits of the data to the model. F values provide the significance of the effect of each predictor in the models. Values equal to or lower than 0.05 represent a significant effect. R2 is the coefficient of determination; it represents the percentage of variance in the response variable that is explained by the fitted model. Table 1 shows the results of the statistical modeling processes in terms of analysis of deviance. Covariates in the final reduced models were Y and LAT. AR was dropped because AIC values did not decrease in any of the models that included it. Figure 2 shows the predicted values of the models as a function of NC and HP. In general, the performance and prediction accuracy of all models were satisfactory, since R2 ranged between 60 % and 93 %. Parameter estimation showed that β2, of main interest for our study, reduced over the years from 1.14±(0.06) in 1985–1989 to 1.05±(0.02) in 2014–2019 (Fig. 3) (values in brackets represent 95 % confidence intervals). This is a reduction of 8 % in the effect of NC on HP. These values can be interpreted in such a way that during 1985–1989, a 1 % increase in NC caused an increase of 1.1 % to 1.14 % in HP, whilst in 2014–2019, a 1 % increase in NC caused an increase of 1.03 % to 1.07 % in HP. Figure 2Real HP (honey production) values (black solid dots) and predicted HP values (red open circles) from GLMs plotted against real NC values in logarithmic scale for the seven sub-series of years. Figure 3Decay of β2 over the years. Vertical bars represent 95 % confidence intervals estimated from GLMs. Figure 4Honey bee and food dynamics represented in the model (Eq. 5). Dashed lines represent food consumption and mortality.

3 The smooth mathematical model approach

In this section we formulate a smooth mathematical model that describes the interaction among immature bees at time t(B(t)), adult bees at time t(T(t)) and the amount of honey produced at time t(M(t)). We assumed that immature bees grow at a rate β, proportionally to the number of adult bees. This is represented by term $T/\left(T+\mathit{\nu }\right)$, where ν is the average saturation rate (number of adult bees needed for immature bees to reach the half of its maximum number). The number of bees reaching the adult stage affects the number of immature bees. This is modeled by term ωB, where ω represents the maturation rate to the adult stage, and $\mathrm{1}/\mathit{\omega }$ represents the time spent before reaching the adult state. Additionally, the number of immature bees is reduced by natural mortality. This is modeled by term μBB, where μB represents the natural death rate of the immature stage. In this way, the following ODE is proposed to model the variation in the number of immature bees at time t:

$\begin{array}{}\text{(2)}& \frac{\mathrm{d}B}{\mathrm{d}t}=\mathit{\beta }\frac{T}{T+\mathit{\nu }}-\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{B}}B.\end{array}$

Similarly, we assumed that the number of adult bees decreases naturally. This is represented by term μTT, where μT is the natural death rate of the adult stage. Additionally, the bees can also die due to a stress factor. This was modeled by the term σT, where σ is the death rate due to a stress factor (loss of habitat, pesticides, climate change, mismanagement by the beekeeper, among others) on bees at the adult stage. Thus, the number of adult bees at time t can be represented by the following ODE:

$\begin{array}{}\text{(3)}& \frac{\mathrm{d}T}{\mathrm{d}t}=\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{T}}T-\mathit{\sigma }T.\end{array}$

Finally, the production of honey in the hives increases at a rate ρ, which depends on the number of adult bees. This is represented by the term $T/\left(T+u\right)$, where u is the average saturation rate. The amount of produced honey decreases due to different factors; one of them is the loss of honey due to natural causes such as support for the hive and feeding of immature bees, which is represented by the term αM, where α is the rate of honey loss. Another cause is the consumption by adult bees, which is modeled by term δTM, where δ is the rate of honey consumption by adult bees. Therefore, the following ODE models the variation in the proportion of existing honey at time t:

$\begin{array}{}\text{(4)}& \frac{\mathrm{d}M}{\mathrm{d}t}=\mathit{\rho }\frac{T}{T+u}-\mathit{\alpha }M-\mathit{\delta }TM.\end{array}$

From Eqs. (2)–(4) we obtain the following system of ODEs:

$\begin{array}{}\text{(5)}& \left\{\begin{array}{ll}\frac{\mathrm{d}B}{\mathrm{d}t}& =\mathit{\beta }\frac{T}{T+\mathit{\nu }}-\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{B}}B\\ \frac{\mathrm{d}T}{\mathrm{d}t}& =\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{T}}T-\mathit{\sigma }T\\ \frac{\mathrm{d}M}{\mathrm{d}t}& =\mathit{\rho }\frac{T}{T+u}-\mathit{\alpha }M-\mathit{\delta }TM\end{array}\right\.\end{array}$

Model (Eq. 5) can be represented in Fig. 4.

The specifications of parameters involved in the model (Eq. 5) are shown in Table 2.

Table 2Parameters involved in the model (Eq. 5): description and dimension. Let us define $\mathit{ϵ}=\mathrm{min}\mathit{\left\{}{\mathit{\mu }}_{\mathrm{B}},{\mathit{\mu }}_{\mathrm{T}}\mathit{\right\}}$. Then, the set of biological interest of the model (Eq. 5) is given by

$\begin{array}{}\text{(6)}& \mathrm{\Omega }=\left\{\left(B,T,M\right)\in {\mathbb{R}}_{+}^{\mathrm{3}}:\mathrm{0}\le B+T\le \frac{\mathit{\beta }}{\mathit{ϵ}},\phantom{\rule{0.25em}{0ex}}\mathrm{0}\le M\le \frac{\mathit{\rho }}{\mathit{\alpha }}\right\}.\end{array}$

The following lemma ensures that set Ω has biological sense; that is, all solutions starting in Ω remain there for all t≥0.

Lema 3.1. The set Ω defined in Eq. (6) is positively invariant for the solutions of the system (Eq. 5).

The proof of the above lema can be seen in Appendix A.

## 3.1 Equilibrium points

For calculating the equilibrium points and stability of the model (Eq. 5), it should be clarified that all parameters given in Table 2 are strictly greater than zero. The equilibrium points of the system (Eq. 5) are given for the solution of the following system of algebraic equations

$\begin{array}{}\text{(7)}& \left\{\begin{array}{l}\mathrm{0}=\mathit{\beta }\frac{T}{T+\mathit{\nu }}-\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{B}}B\\ \mathrm{0}=\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{T}}T-\mathit{\sigma }T\\ \mathrm{0}=\mathit{\rho }\frac{T}{T+u}-\mathit{\alpha }M-\mathit{\delta }TM\end{array}\right\.\end{array}$

Details about the solution of the previous system can be found in Appendix A. We get the following proposition.

Proposition 3.2. Model (Eq. 5) always has a trivial equilibrium ${\mathbf{E}}_{\mathrm{0}}=\left(\mathrm{0},\mathrm{0},\mathrm{0}\right)$. Besides, if ${\mathit{\lambda }}^{*}>\mathit{\nu }$ there is a non-trivial equilibrium ${\mathbf{E}}_{\mathrm{1}}=\left({b}^{*},{\mathit{\lambda }}^{*}-\mathit{\nu },{m}^{*}\right)$, where b* and m* are defined in Eqs. (A4) and (A5), respectively.

## 3.2 Local stability analysis

Here, we determine the conditions for the local stability of the equilibrium solutions given in Proposition 3.2. We use the linearization of the vector field given by the right-hand side of the system (Eq. 5), which is given by the following Jacobian matrix:

$\begin{array}{}\text{(8)}& \mathbf{J}\left(\mathbf{E}\right)=\left(\begin{array}{lll}-\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)& \frac{\mathit{\beta }\mathit{\nu }}{\left(T+\mathit{\nu }{\right)}^{\mathrm{2}}}& \mathrm{0}\\ \mathit{\omega }& -\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)& \mathrm{0}\\ \mathrm{0}& \frac{\mathit{\rho }u}{\left(T+u{\right)}^{\mathrm{2}}}-\mathit{\delta }M& -\left(\mathit{\alpha }+\mathit{\delta }T\right)\end{array}\right).\end{array}$

The stability of each equilibrium point given on Proposition 3.2 can be found by evaluating the above matrix in each equilibrium point. To see the mathematical details see Appendix A. We obtain the following proposition.

Proposition 3.3. If ${\mathit{\lambda }}^{*}<\mathit{\nu }$, the trivial equilibrium ${\mathbf{E}}_{\mathrm{0}}=\left(\mathrm{0},\mathrm{0},\mathrm{0}\right)$ is locally and asymptotically stable, whereas if ${\mathit{\lambda }}^{*}>\mathit{\nu }$, it is unstable. If ${\mathit{\lambda }}^{*}=\mathit{\nu }$, E0 is a non-hyperbolic equilibrium point. Similarly, If ${\mathit{\lambda }}^{*}>\mathit{\nu }$, the coexistence equilibrium ${\mathbf{E}}_{\mathrm{1}}=\left({b}^{*},{\mathit{\lambda }}^{*}-\mathit{\nu },{m}^{*}\right)$ is locally and asymptotically stable, whereas if ${\mathit{\lambda }}^{*}<\mathit{\nu }$, it is unstable. If ${\mathit{\lambda }}^{*}=\mathit{\nu }$, E1 is a non-hyperbolic equilibrium point.

A summary of the results from Proposition 3.3 is given in Table 3.

Table 3Condition of the existence and stability of the equilibrium points of the model (Eq. 5). Remark 3.4. The existence and stability conditions of the equilibrium points of the system (Eq. 5) are given in terms of an ecological threshold λ* defined in Eq. (A3), which can be rewritten as

${\mathit{\lambda }}^{*}=\frac{\mathit{\omega }}{{\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }}*\frac{\mathit{\beta }}{\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}}.$

Since μT+σ represents the total death rate of adult bees and ω+μB the total death rate of immature bees, the expression $\mathit{\omega }/{\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }$ can be interpreted as the proportion of immature bees that equals the number of adult bees when these remain constant. Similarly, $\mathit{\beta }/\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}$ is the proportion of adult bees that equals the number of immature bees when these remain constant. Thus, λ* is the total net number of bees produced by the colony during the lifetime of the hive. If the total death rate of adult bees is large enough, it would lead to the collapse of the hive because the term λ* would approach to zero. This is because the existence condition of the coexistence equilibrium states that ${\mathit{\lambda }}^{*}>\mathit{\nu }$, where the constant ν is the number of adult bees required for the growth rate of immature bees to reach half its maximum value.

## 3.3 Numerical simulations

In this section we present some numerical experiments that illustrate the behavior of the hive and honey subjected to different parameter values. Some values were taken from the literature and others were assumed. In Table 4 we show the rank of the parameter values included in the model and their respective references.

In Fig. 5 we present the stability of the equilibrium E0, which is obtained with a higher value for the reproduction rate of immature bees. Here, σ=0.8 (considered a high value), ρ=0.13, μB=0.11, μT=0.29, ω=0.95, α=0.95, δ=0.571, u=1, ν=1 and β=0.92. With these values the number of bees and their honey will constantly decrease until a collapse. Figure 5Numerical representation of the equilibrium E0. Here, ${\mathit{\lambda }}^{*}=\mathrm{0.7564}$ and ν=1.

Figure 6 numerically represents the stability of the equilibrium E1. Here, the mortality rate due to a stress factor was reduced to an average value of σ=0.4, and the honey production rate was increased to ρ=0.23. The loss of honey due to natural causes has also increased to α=0.1. By including the above changes in parameter values, the amounts of immature adults and honey production remain constant in B=0.14, T=0.2 and M=0.289. However, this is less likely to occur because, due to multiple factors, the hives tend to collapse after a certain time. Figure 6Numerical representation of the equilibrium E1. Here, ${\mathit{\lambda }}^{*}=\mathrm{1.195}$ and ν=1.

4 The non-smooth mathematical model approach

Some mathematical models whose arguments are described by piecewise smooth functions have numerous dynamics and applications, as well as an extensive mathematical structure, even though their behavior is not simply described in terms of modern qualitative theory of dynamical systems. However, this classical theory of dynamic systems does not include a significant number of phenomena that occur in practice. These phenomena are defined with functions that are non-smooth in their arguments, such as electrical circuits in which there are commutations (switches), control systems, models in different sciences and population models in biology, where some continuous changes can generate discrete actions. Such is the case of bee populations, in which different stress factors can cause some thresholds of mortality and affect both the number of bees in the immature state, as well as in the adult state, and therefore honey production.

Thus, in this section we use the Filippov systems to study the system (Eq. 5), where R1 and R2 are considered as death thresholds. If the number T of adult bees is big enough from a R2, the death of bees due to a stress factor (σ) will be maximum (σmax). Similarly, from a R1, if the number of adult bees is small enough, the death rate due to a stress factor will be minimal (σmin). Finally, if we consider that the number of adult bees is in equilibrium, that is between R1 and R2, the death rate due to a stress factor will be stable (σ1). The previous consideration can be mathematically written as follows:

$\begin{array}{}\text{(9)}& \mathit{\sigma }=\left\{\begin{array}{ll}{\mathit{\sigma }}_{\mathrm{max}}& \mathrm{if}\phantom{\rule{0.25em}{0ex}}T>{R}_{\mathrm{2}},\\ {\mathit{\sigma }}_{\mathrm{1}}& \mathrm{if}\phantom{\rule{0.25em}{0ex}}{R}_{\mathrm{1}}\le T\le {R}_{\mathrm{2}},\\ {\mathit{\sigma }}_{\mathrm{min}}& \mathrm{if}\phantom{\rule{0.25em}{0ex}}T<{R}_{\mathrm{1}}.\end{array}\right\\end{array}$

The main characteristic of the Filippov systems is the subdivision of the state space into disjoint subregions, and a smooth vector field is defined within each one of them . The limits between the different regions are called discontinuity surfaces, which are usually denoted as Σi .

Mathematical details about this approach can be seen in Appendix B. However, this presents only a mathematical understanding of Filippov's systems; a biological sense of this will be determined in future work, as well as the study of the presence of pseudo-equilibrium points.

## 4.1 Numerical simulations

The numerical simulations of this section were performed using the MATLAB interface, and the codes can be consulted at the following GitHub link: https://github.com/jpatirom3/Honey-bees-modeling- (last access: March 2022). Numerical simulations for the described systems are presented below. Figure 7 shows the trajectory of the solutions of the states B(t), T(t) and M(t). In this figure, we take B(0)=1, T(0)=1 and M(0)=1 as initial conditions. Table 5 summarizes the values used for each parameter. Figure 7Trajectories for B, T and M. Here the initial condition is $\left(B\left(\mathrm{0}\right),T\left(\mathrm{0}\right),M\left(\mathrm{0}\right)\right)=\left(\mathrm{1},\mathrm{1},\mathrm{1}\right)$ and β=0.9, ω=0.95, ρ=0.23, μB=0.1, μT=0.29, α=0.018, u=1, ν=1 and δ=0.5. The ecological threshold is ${\mathit{\lambda }}^{*}=\mathrm{1.2063}$. Figure 8Numerical simulation of the model (Eq. 9) with initial conditions – number of immature bees: 1; number of adult bees: 0.15, 0.35 and 0.7; and amount of honey: 1. The trajectories are represented by the red, blue and green curves, respectively. The parameter values are given in Table 5.

Figure 8 shows the switching surfaces Σ1 (left side) and Σ2 (right side), which divide the first octant into the following regions: S1 (left side of Σ1), S2 (region between Σ1 and Σ2), and S3 (right side of Σ2), as well as the phase portrait. The trajectories show that regardless of the initial conditions or the region of start, there is always a crossover in the switching planes, which is proven analytically with the Filippov analysis as shown above. After some time, the trajectories remain in regions S2 and S3, which means that the number of adult bees is switching within these regions limited by the switching surface Σ2. Finally, the number of adult bees seeks stability at the break-even point, as shown by Fig. 7.

5 Discussion

In this work, we attempted to contrast three different modeling approaches to understand the population dynamics of immature bees, adult bees and the honey produced in the hive. We first formulated a smooth mathematical model in which conditions of existence and stability of equilibrium solutions were found. Afterwards, we adapted the previous model to a non-smooth model which was studied using Filippov's systems theory. Finally, we analyzed a honey bee data set obtained from USDA, using a GLM.

From the smooth model, we obtained two equilibrium points: (1) a coexistence between the number of bees and honey production and (2) a collapse in both the number of bees and honey production. An important aspect to highlight here is that the existence and stability of these equilibrium points are in terms of the ecological threshold λ*, given in Remark 3.4, which represents the net number of bees produced within the colony over the lifetime of the hive. Additionally, the numerical experiments of these two points suggested that a high death rate due to a stress factor (σ=0.8), which was used in Fig. 5, leads to a collapse of the hive in 50 weeks.

Subsequently, a smooth model was proposed in which two death thresholds R1 and R2 were considered for the stress factor. The Filippov analysis for the non-smooth model showed that there is always a crossover in the switching surfaces. After 100 weeks, the trajectories switched between regions S2 and S3, which defined a normal and high stress rate respectively, as it was shown in Fig. 8. Finally, the orbits stabilized at the equilibrium point located in region S2. We have taken into account the fact that the equilibrium point E1 is stable, and, therefore, the equilibrium value is reached when t→∞.

For the statistical approach, generalized linear models were built to study seven 5-year sub-series of colony abundance and honey production data (1985 to 2019) from the USDA. We observed that the number of colonies required to produce a same amount of honey is not constant but reduced over the years as shown by Fig. 3. This points to the influence of a stress factor, such as the change of natural habitats to mono-cultures or the infestation by the ectoparasitic Varroa mite destructor .

From these three modeling approaches, we conclude that honey production is reduced due to the decrease in the total number of bees within the colony and due to an external stress factor, which could be climate variability, pesticides, parasites such as the Varroa mite, or habitat loss, among other things. This can be observed in Fig. 6, in the simulations of the non-smooth model (Eqs. 7 and 8), and mainly in the statistical results in Fig. 3. From numerical experiments it was observed that the amount of honey produced is considerably reduced in time; however, it does not fully collapse but stay constant.

Appendix A

Proof of Lema 3.1.

Proof. It is clear that $\frac{T}{T+\mathit{\nu }}\le \mathrm{1}$. Adding the first two equations of the system (Eq. 5) we obtain

$\begin{array}{ll}\frac{\mathrm{d}\left(B+T\right)}{\mathrm{d}t}& =\mathit{\beta }\frac{T}{T+\mathit{\nu }}-{\mathit{\mu }}_{\mathrm{B}}B-{\mathit{\mu }}_{\mathrm{T}}T-\mathit{\sigma }T\\ & \le \mathit{\beta }-{\mathit{\mu }}_{\mathrm{B}}B-{\mathit{\mu }}_{\mathrm{T}}T\\ & \le \mathit{\beta }-\mathit{ϵ}\left(B+T\right).\end{array}$

Thus, $B+T\le \mathit{\beta }/\mathit{ϵ}$. Similarly, for M we get that

$\begin{array}{ll}\frac{\mathrm{d}M}{\mathrm{d}t}& =\mathit{\rho }\frac{T}{T+u}-\mathit{\alpha }M-\mathit{\delta }TM.\\ & \le \mathit{\rho }-\mathit{\alpha }M-\mathit{\delta }TM\\ & \le \mathit{\rho }-\mathit{\alpha }M.\end{array}$

Therefore, $M\le \mathit{\rho }/\mathit{\alpha }$, which concludes the proof.

## A1 Equilibrium points

To solve the algebraic system given on Eq. (7), we proceed as follows.

From the second equation in Eq. (7),

$\begin{array}{}\text{(A1)}& B=\left(\frac{{\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }}{\mathit{\omega }}\right)T.\end{array}$

Replacing B defined in Eq. (A1) in the first equation of the system (Eq. 7) we obtain

$T\left(\frac{\mathit{\beta }}{T+\mathit{\nu }}-\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)\left(\mathrm{1}+\frac{{\mathit{\mu }}_{\mathrm{B}}}{\mathit{\omega }}\right)\right)=\mathrm{0}.$

From the previous equation we have T=0 or T≠0. First, let us suppose that T≠0, thus

$\begin{array}{}\text{(A2)}& T={\mathit{\lambda }}^{*}-\mathit{\nu },\end{array}$

where

$\begin{array}{}\text{(A3)}& {\mathit{\lambda }}^{*}=\frac{\mathit{\omega }\mathit{\beta }}{\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)}.\end{array}$

Replacing T defined on Eq. (A2) into the third equation of the system, Eq. (7), we have that

$\mathit{\rho }\frac{{\mathit{\lambda }}^{*}-\mathit{\nu }}{{\mathit{\lambda }}^{*}-\mathit{\nu }+u}-\mathit{\alpha }M-\left({\mathit{\lambda }}^{*}-\mathit{\nu }\right)\mathit{\delta }M=\mathrm{0},$

from where

$\begin{array}{}\text{(A4)}& M=\frac{\mathit{\rho }\left({\mathit{\lambda }}^{*}-\mathit{\nu }\right)}{\left(\mathit{\alpha }+\left({\mathit{\lambda }}^{*}-\mathit{\nu }\right)\mathit{\delta }\right)\left({\mathit{\lambda }}^{*}-\mathit{\nu }-u\right)}=:{m}^{*}.\end{array}$

Now, replacing T defined on Eq. (A2) into the Eq. (A1), we obtain that

$\begin{array}{}\text{(A5)}& B=\frac{\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)\left({\mathit{\lambda }}^{*}-\mathit{\nu }\right)}{\mathit{\omega }}=:{b}^{*}.\end{array}$

Note that a sufficient condition for T defined in Eq. (A2), M defined in Eq. (A4) and B defined in Eq. (A5) to be positive is that ${\mathit{\lambda }}^{*}>\mathit{\nu }.$

Now, if T=0 in Eq. (A2), we obtain the trivial equilibrium solution ${E}_{\mathrm{0}}=\left(\mathrm{0},\mathrm{0},\mathrm{0}\right)$.

## A2 Local stability analysis

The stability of E0 can be found by evaluating the above matrix in E0. Thus,

$\mathbf{J}\left({\mathbf{E}}_{\mathrm{0}}\right)=\left(\begin{array}{lll}-\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)& \frac{\mathit{\beta }}{\mathit{\nu }}& \mathrm{0}\\ \mathit{\omega }& -\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)& \mathrm{0}\\ \mathrm{0}& \frac{\mathit{\rho }}{u}& -\mathit{\alpha }\end{array}\right).$

An eigenvalue of the previous matrix is $\mathit{\xi }=-\mathit{\alpha }$, whereas the other two eigenvalues are given by the roots of the following characteristic polynomial:

$\begin{array}{ll}p\left(\mathit{\xi }\right)& ={\mathit{\xi }}^{\mathrm{2}}+\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}+{\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }\right)\mathit{\xi }\\ & +\left[\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }\right)-\frac{\mathit{\omega }\mathit{\beta }}{\mathit{\nu }}\right]\\ & ={\mathit{\xi }}^{\mathrm{2}}+\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}+{\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }\right)\mathit{\xi }\\ & +\frac{\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }\right)}{\mathit{\nu }}\left(\mathit{\nu }-{\mathit{\lambda }}^{*}\right)\\ & ={\mathit{\xi }}^{\mathrm{2}}+{a}_{\mathrm{1}}\mathit{\xi }+{a}_{\mathrm{0}}.\end{array}$

From the Routh–Hurwitz criterion , we have that the roots of the above polynomial have negative real parts, if and only if the following conditions are met:

$\begin{array}{ll}& {\mathrm{\Delta }}_{\mathrm{1}}=\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}+{\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }\right)>\mathrm{0}\\ & {\mathrm{\Delta }}_{\mathrm{2}}=\left|\begin{array}{cc}\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}+{\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }& \mathrm{1}\\ \mathrm{0}& \frac{\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }\right)}{\mathit{\nu }}\left(\mathit{\nu }-{\mathit{\lambda }}^{*}\right)\end{array}\right|>\mathrm{0}.\end{array}$

It is clear that Λ1>0 and Λ2>0, if and only if ${\mathit{\lambda }}^{*}<\mathit{\nu }.$

Now, we will determine the stability conditions for the equilibrium E1. Evaluating matrix (Eq. 8) in E1 we obtain

$\begin{array}{ll}& \mathbf{J}\left({\mathbf{E}}_{\mathrm{1}}\right)=\\ & \left(\begin{array}{lll}-\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)& \frac{\mathit{\beta }\mathit{\nu }}{{\left({\mathit{\lambda }}^{*}\right)}^{\mathrm{2}}}& \mathrm{0}\\ \mathit{\omega }& -\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)& \mathrm{0}\\ \mathrm{0}& \frac{\mathit{\rho }u}{{\left({\mathit{\lambda }}^{*}-\mathit{\nu }+u\right)}^{\mathrm{2}}}-\mathit{\delta }{m}^{*}& -\left(\mathit{\alpha }+\mathit{\delta }\left({\mathit{\lambda }}^{*}-\mathit{\nu }\right)\right)\end{array}\right).\end{array}$

An eigenvalue of the previous matrix is $\mathit{\eta }=-\mathit{\alpha }-\mathit{\delta }\left({\mathit{\lambda }}^{*}-\mathit{\nu }\right)$, which is negative if ${\mathit{\lambda }}^{*}\ge \mathit{\nu }$, while the other two eigenvalues are given by the roots of the following characteristic polynomial:

$\begin{array}{ll}p\left(\mathit{\eta }\right)=& {\mathit{\eta }}^{\mathrm{2}}+\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}+{\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)\mathit{\eta }\\ & +\left[\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)-\frac{\mathit{\beta }\mathit{\nu }\mathit{\omega }}{{\left({\mathit{\lambda }}^{*}\right)}^{\mathrm{2}}}\right]\\ & ={\mathit{\eta }}^{\mathrm{2}}+\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}+{\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }\right)\mathit{\eta }\\ & +\frac{\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)}{\mathit{\nu }}\left[\mathit{\nu }-\frac{{\mathit{\nu }}^{\mathrm{2}}}{{\mathit{\lambda }}^{*}}\right]\\ & ={\mathit{\eta }}^{\mathrm{2}}+\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}+{\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }\right)\mathit{\eta }\\ & +\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)\left[\mathrm{1}-\frac{\mathit{\nu }}{{\mathit{\lambda }}^{*}}\right]\\ & ={\mathit{\eta }}^{\mathrm{2}}+{b}_{\mathrm{1}}\mathit{\eta }+{b}_{\mathrm{0}}.\end{array}$

Similarly, from the Routh–Hurwitz criterion the roots of the above polynomial have negative real parts, if and only if the following conditions are met:

$\begin{array}{ll}& {\mathrm{\Theta }}_{\mathrm{1}}=\left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}+{\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }\right)>\mathrm{0}\\ & {\mathrm{\Theta }}_{\mathrm{2}}=\left|\begin{array}{cc}\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}+{\mathit{\mu }}_{\mathrm{T}}+\mathit{\delta }& \mathrm{1}\\ \mathrm{0}& \left(\mathit{\omega }+{\mathit{\mu }}_{\mathrm{B}}\right)\left({\mathit{\mu }}_{\mathrm{T}}+\mathit{\sigma }\right)\left(\mathrm{1}-\frac{\mathit{\nu }}{{\mathit{\lambda }}^{*}}\right)\end{array}\right|>\mathrm{0}.\end{array}$

It is clear that Θ1>0 while Θ2>0, if and only if ${\mathit{\lambda }}^{*}>\mathit{\nu }.$

Appendix B

We assume that the state space consists of the three regions S1, S2 and S3⊂ℝ3, which are separated by two discontinuity surfaces Σi with $i=\mathrm{1},\mathrm{2}$, defined as the null set of a smooth scalar function $H\left(x\right):{\mathbb{R}}^{\mathrm{3}}\to \mathbb{R}$, in such a way that the equations (Eq. 5) under the conditions (Eq. 9) are the following smooth functions: ${F}_{\mathrm{1}}\left(B,T,M\right)$, ${F}_{\mathrm{2}}\left(B,T,M\right)$ and ${F}_{\mathrm{3}}\left(B,T,M\right)$ determined by the value of T according to Eq. (9). If T<R1, we have

${F}_{\mathrm{1}}\left(B,T,M\right)=\left\{\begin{array}{l}\mathit{\beta }\frac{T}{T+\mathit{\nu }}-\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{B}}B\\ \mathit{\omega }B-{\mathit{\mu }}_{\mathrm{T}}T-{\mathit{\sigma }}_{\mathrm{min}}T\\ \mathit{\rho }\frac{T}{T+u}-\mathit{\alpha }M-\mathit{\delta }TM\end{array}\right\.$

When ${R}_{\mathrm{1}}\le T\le {R}_{\mathrm{2}}$, we have that

${F}_{\mathrm{2}}\left(B,T,M\right)=\left\{\begin{array}{l}\mathit{\beta }\frac{T}{T+\mathit{\nu }}-\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{B}}B\\ \mathit{\omega }B-{\mathit{\mu }}_{\mathrm{T}}T-{\mathit{\sigma }}_{\mathrm{1}}T\\ \mathit{\rho }\frac{T}{T+u}-\mathit{\alpha }M-\mathit{\delta }TM\end{array}\right\.$

Finally, at T>R2

${F}_{\mathrm{3}}\left(B,T,M\right)=\left\{\begin{array}{l}\mathit{\beta }\frac{T}{T+\mathit{\nu }}-\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{B}}B\\ \mathit{\omega }B-{\mathit{\mu }}_{\mathrm{T}}T-{\mathit{\sigma }}_{\mathrm{max}}T\\ \mathit{\rho }\frac{T}{T+u}-\mathit{\alpha }M-\mathit{\delta }TM\end{array}\right\.$

These three functions are also assumed to be defined through the state space, even though they are only used in their respective regions . Thus, the regions S1, S2 and S3 can be defined by the following sets:

$\begin{array}{ll}& {S}_{\mathrm{1}}=\left\{\left(B,T,M\right)\in {\mathbb{R}}^{\mathrm{3}}:T<{R}_{\mathrm{1}}\right\}\\ & {S}_{\mathrm{2}}=\left\{\left(B,T,M\right)\in {\mathbb{R}}^{\mathrm{3}}:{R}_{\mathrm{1}}\le T\le {R}_{\mathrm{2}}\right\}\\ & {S}_{\mathrm{3}}=\left\{\left(B,T,M\right)\in {\mathbb{R}}^{\mathrm{3}}:T>{R}_{\mathrm{2}}\right\}.\end{array}$

The regions are separated by a switching surface: limit of discontinuity represented by the scalar function ${H}^{\left(i\right)}:{\mathbb{R}}^{\mathrm{3}}⟶\mathbb{R}$; ${H}^{\left(i\right)}\left(x\right)=T-{R}_{i}$ with $i=\mathrm{1},\mathrm{2}$, so that

${\mathrm{\Sigma }}_{i}=\left\{x\in {\mathbb{R}}^{\mathrm{3}}:T={R}_{i},\phantom{\rule{0.25em}{0ex}}i=\mathrm{1},\mathrm{2}\right\}.$

As $x=\left(B,T,M\right)\in {R}^{\mathrm{3}}$, the gradient of H(i)(x) can be defined as follows:

${H}_{x}^{\left(i\right)}\left(x\right)=\left(\mathrm{0},\mathrm{1},\mathrm{0}\right),\phantom{\rule{0.25em}{0ex}}i=\mathrm{1},\mathrm{2}.$

Thus, taking the regions Σi in which the gradient does not null in the set ${\mathrm{\Gamma }}_{i}=\mathit{\left\{}x\in {\mathbb{R}}^{\mathrm{3}}:{H}^{\left(i\right)}\left(x\right)=\mathrm{0},\phantom{\rule{0.25em}{0ex}}i=\mathrm{1},\mathrm{2}\mathit{\right\}}$, for x∈Γ1, we define

$\begin{array}{ll}{\mathit{\phi }}^{\left(\mathrm{1}\right)}\left(x\right)& =〈{H}_{x}^{\left(\mathrm{1}\right)}\left(x\right),{F}_{\mathrm{1}}\left(x\right)〉〈{H}_{x}^{\left(\mathrm{1}\right)}\left(x\right),{F}_{\mathrm{2}}\left(x\right)〉\\ & =\left(\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{min}}\right)T\right)\left(\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T\right),\end{array}$

and, for x∈Γ2,

$\begin{array}{ll}{\mathit{\phi }}^{\left(\mathrm{2}\right)}\left(x\right)& =〈{H}_{x}^{\left(\mathrm{2}\right)}\left(x\right),{F}_{\mathrm{2}}\left(x\right)〉〈{H}_{x}^{\left(\mathrm{2}\right)}\left(x\right),{F}_{\mathrm{3}}\left(x\right)〉\\ & =\left(\mathit{\omega }B-\left({\mathit{\mu }}_{T}+{\mathit{\sigma }}_{\mathrm{1}}\right)T\right)\left(\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{max}}\right)T\right).\end{array}$

Filippov's systems theory states that at points x, where

$\begin{array}{ll}{\mathit{\phi }}^{\left(i\right)}\left(x\right)>\mathrm{0},& i=\mathrm{1},\mathrm{2},\phantom{\rule{0.25em}{0ex}}\mathrm{there}\phantom{\rule{0.25em}{0ex}}\mathrm{is}\phantom{\rule{0.25em}{0ex}}\mathrm{a}\phantom{\rule{0.25em}{0ex}}\mathrm{crossing}\phantom{\rule{0.25em}{0ex}}\mathrm{point},\\ {\mathit{\phi }}^{\left(i\right)}\left(x\right)=\mathrm{0},& i=\mathrm{1},\mathrm{2},\phantom{\rule{0.25em}{0ex}}\mathrm{there}\phantom{\rule{0.25em}{0ex}}\mathrm{are}\phantom{\rule{0.25em}{0ex}}\mathrm{tangential}\phantom{\rule{0.25em}{0ex}}\mathrm{points},\\ {\mathit{\phi }}^{\left(i\right)}\left(x\right)<\mathrm{0},& i=\mathrm{1},\mathrm{2},\phantom{\rule{0.25em}{0ex}}\mathrm{there}\phantom{\rule{0.25em}{0ex}}\mathrm{is}\phantom{\rule{0.25em}{0ex}}\mathrm{a}\phantom{\rule{0.25em}{0ex}}\mathrm{sliding}\phantom{\rule{0.25em}{0ex}}\mathrm{point},\\ & \mathrm{with}\phantom{\rule{0.25em}{0ex}}\mathrm{the}\phantom{\rule{0.25em}{0ex}}\mathrm{equation}\phantom{\rule{0.25em}{0ex}}\stackrel{\mathrm{˙}}{x}={g}_{i}\left(x\right),\phantom{\rule{0.25em}{0ex}}x\in {\mathrm{\Sigma }}_{i}.\end{array}$

With the vector field

$\begin{array}{ll}& {g}_{i}\left(x\right)={\mathit{\lambda }}_{i}{F}_{i}\left(x\right)+\left(\mathrm{1}-{\mathit{\lambda }}_{i}\right){F}_{i+\mathrm{1}}\left(x\right)\phantom{\rule{0.25em}{0ex}}\mathrm{with}\\ & {\mathit{\lambda }}_{i}=\frac{〈{H}_{x}^{\left(i\right)}\left(x\right),{F}_{i+\mathrm{1}}\left(x\right)〉}{〈{H}_{x}^{\left(i\right)}\left(x\right),{F}_{i+\mathrm{1}}\left(x\right)-{F}_{i}\left(x\right)〉},\phantom{\rule{0.25em}{0ex}}i=\mathrm{1},\mathrm{2}.\end{array}$

For x∈Γ1, we have the following cases.

• Crossing point, when φ(1)(x)>0, i.e.,

$\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{min}}\right)T>\mathrm{0}\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T>\mathrm{0},$

or

$\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{min}}\right)T<\mathrm{0}\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T<\mathrm{0},$

as σmin<σ1, we have

$B<\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{min}}\right)T}{\mathit{\omega }},\phantom{\rule{0.25em}{0ex}}\mathrm{or}\phantom{\rule{0.25em}{0ex}}B>\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T}{\mathit{\omega }}.$
• Tangential point when

$B=\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{min}}\right)T}{\mathit{\omega }}\phantom{\rule{0.25em}{0ex}}\mathrm{or}\phantom{\rule{0.25em}{0ex}}B=\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T}{\mathit{\omega }}.$
• Sliding point at

$\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{min}}\right)T}{\mathit{\omega }}

with the equation $\stackrel{\mathrm{˙}}{x}={g}_{\mathrm{1}}\left(x\right)$, where the vector field is given by

${g}_{\mathrm{1}}\left(x\right)={\mathit{\lambda }}_{\mathrm{1}}{F}_{\mathrm{1}}\left(x\right)+\left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{1}}\right){F}_{\mathrm{1}}\left(x\right),$

with

${\mathit{\lambda }}_{\mathrm{1}}=\frac{\mathit{\omega }}{{\mathit{\sigma }}_{\mathrm{min}}-{\mathit{\sigma }}_{\mathrm{1}}}\left(\frac{B}{T}\right)-\frac{{\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{min}}-{\mathit{\sigma }}_{\mathrm{1}}},$

i.e.,

${g}_{\mathrm{1}}\left(x\right)=\left\{\begin{array}{l}\mathit{\beta }\frac{T}{T+\mathit{\nu }}-\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{B}}B\\ \left[\frac{\mathit{\omega }}{{\mathit{\sigma }}_{\mathrm{min}}-{\mathit{\sigma }}_{\mathrm{1}}}\left(\frac{B}{T}\right)-\frac{{\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{min}}-{\mathit{\sigma }}_{\mathrm{1}}}\right]\left({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{min}}\right)T\\ +\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T\\ \mathit{\rho }\frac{T}{T+u}-\mathit{\alpha }M-\mathit{\delta }TM\end{array}.\right\$

For x∈Γ2, we have the following cases.

• Crossing point, when φ(2)(x)>0, i.e.,

$\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T>\mathrm{0}\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{max}}\right)T>\mathrm{0},$

or

$\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T<\mathrm{0}\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{max}}\right)T<\mathrm{0},$

as σ1<σmax, we have

$B<\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T}{\mathit{\omega }},\phantom{\rule{0.25em}{0ex}}\mathrm{or}\phantom{\rule{0.25em}{0ex}}B>\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{max}}\right)T}{\mathit{\omega }}.$
• Tangential point when

$B=\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T}{\mathit{\omega }}\phantom{\rule{0.25em}{0ex}}\mathrm{or}\phantom{\rule{0.25em}{0ex}}B=\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{max}}\right)T}{\mathit{\omega }}.$
• Sliding point at

$\frac{\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{1}}\right)T}{\mathit{\omega }}

with the equation $\stackrel{\mathrm{˙}}{x}={g}_{\mathrm{2}}\left(x\right)$, where the vector field is given by

${g}_{\mathrm{2}}\left(x\right)={\mathit{\lambda }}_{\mathrm{2}}{F}_{\mathrm{2}}\left(x\right)+\left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{2}}\right){F}_{\mathrm{3}}\left(x\right),$

with

${\mathit{\lambda }}_{\mathrm{2}}=\frac{\mathit{\omega }}{{\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{max}}}\left(\frac{B}{T}\right)-\frac{{\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{max}}}{{\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{1}}},$

i.e.,

${g}_{\mathrm{2}}\left(x\right)=\left\{\begin{array}{l}\mathit{\beta }\frac{T}{T+\mathit{\nu }}-\mathit{\omega }B-{\mathit{\mu }}_{\mathrm{B}}B\\ \left[\frac{\mathit{\omega }}{{\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{max}}}\left(\frac{B}{T}\right)-\frac{{\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{max}}}{{\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{max}}}\right]\left({\mathit{\sigma }}_{\mathrm{max}}-{\mathit{\sigma }}_{\mathrm{1}}\right)T\\ +\mathit{\omega }B-\left({\mathit{\mu }}_{\mathrm{T}}+{\mathit{\sigma }}_{\mathrm{max}}\right)T\\ \mathit{\rho }\frac{T}{T+u}-\mathit{\alpha }M-\mathit{\delta }TM\end{array}\right\.$
Code availability

The R and Matlab codes and the cleaned data set are available at the following GitHub link: https://github.com/jpatirom3/Honey-bees-modeling- .

Data availability

The historical data of bee farming in the United States of America (USA) from 1985 to 2019, from the US Department of Agriculture (USDA), used in this paper are publicly available on the USDA web page (https://usda.library.cornell.edu/concern/publications/hd76s004z?locale=en&page=3#release-items, ).

Author contributions

JPR performed the conceptualization, model formulation, formal analysis, supervision, validation, and writing of the original draft and developed the methodology and software. AG performed the investigation and data collection, developed the methodology, and participated in the review and editing. IFB performed the data curation, visualization, and statistical analysis; developed the R code; and participated in the review and editing. OM performed the non-smooth model analysis, developed the MATLAB code, and participated in the review and editing. AP performed the model formulation, carried out the funding and resource acquisition, participated in the project administration and in the review and editing.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

Jhoana P. Romero-Leiton and Ivan Felipe Benavides thank the support of Fundación Ceiba (Colombia). All authors are grateful to Victor Solarte from Blue Note Data Analysis SAS and to the anonymous referees for the recommendations that helped to improve this article.

Review statement

This paper was edited by John M. Halley and reviewed by Sergio Estay and one anonymous referee.

References

Amador, J., Granada, H., Redondo, J., and Olivar, G.: Nonlinear and Nonsmooth Dynamics in Stress-Sickness Processes, Ciencia en Desarrollo, 8, 9–19, 2017. a

Amador, J. A., Olivar, G., and Angulo, F.: Smooth and Filippov models of sustainable development: Bifurcations and numerical computations, Different. Equat. Dynam. Syst., 21, 173–184, https://doi.org/10.1007/s12591-012-0138-2, 2013. a

Amador, J. A., Redondo, J. M., Olivar-Tost, G., and Erazo, C.: Cooperation-Based Modeling of Sustainable Development: An Approach from Filippov's Systems, Complexity, https://doi.org/10.1155/2021/4249106, 2021. a

Bagheri, S. and Mirzaie, M.: A mathematical model of honey bee colony dynamics to predict the effect of pollen on colony failure, PloS One, 14, 1–19, https://doi.org/10.1371/journal.pone.0225632, 2019. a

Belić, M., Škarka, V., Deneubourg, J., and Lax, M.: Mathematical model of honeycomb construction, J. Math. Biol., 24, 437–449, https://doi.org/10.1007/BF01236891, 1986. a

Belsky, J. and Joshi, N.: Impact of biotic and abiotic stressors on managed and feral bees, Insects, 10, 233, https://doi.org/10.3390/insects10080233, 2019. a

Booton, R., Iwasa, Y., Marshall, J., and Childs, D.: Stress-mediated Allee effects can cause the sudden collapse of honey bee colonies, J. Theor. Biol., 420, 213–219, https://doi.org/10.1016/j.jtbi.2017.03.009, 2017. a

Bretagnolle, V. and Gaba, S.: Weeds for bees? A review, Agron. Sustain. Dev., 35, 891–909, https://doi.org/10.1007/s13593-015-0302-5, 2015. a

Brodschneider, R. and Crailsheim, K.: Nutrition and health in honey bees, Apidologie, 41, 278–294, https://doi.org/10.1051/apido/2010012, 2010. a

Brown, M. and Paxton, R.: The conservation of bees: a global perspective, Apidologie, 40, 410–416, https://doi.org/10.1051/apido/2009019, 2009. a

Campos, M. and Leyva, C.: El mercado internacional de la miel de abeja y la competitividad de México, Revista de Economía, 35, 87–123, https://doi.org/10.33937/reveco.2018.92, 2018. a

Cox-Foster, D. and VanEngelsdorp, D.: Saving the honeybee, Scient. Am., 300, 40–47, https://doi.org/10.1038/scientificamerican0409-40, 2009. a

DeJesus, E. and Kaufman, C.: Routh-Hurwitz criterion in the examination of eigenvalues of a system of nonlinear ordinary differential equations, Phys. Rev. A, 35, 5288–5290, https://doi.org/10.1103/PhysRevA.35.5288, 1987. a, b

Dennis, B. and Kemp, W.: How hives collapse: Allee effects, ecological resilience, and the honey bee, PloS One, 11, 1–17, https://doi.org/10.1371/journal.pone.0150055, 2016. a

Dieci, L. and Lopez, L.: Sliding motion in Filippov differential systems: theoretical results and a computational approach, SIAM J. Numer. Anal., 47, 2023–2051, 2009. a, b

Farouk, K., Palmera, K., and Sepúlveda, P.: Abejas, InfoZoa Boletín de zoología, 6, 1–12, 2014. a, b

García, M., Ríos, L., and Álvarez, J.: La polinización en los sistemas de producción agrícola, Revisión sistemática de la literatura, Idesia (Arica), 34, 53–68, https://doi.org/10.4067/S0718-34292016000300008, 2016. a

García, N.: The current situation on the international honey market, Bee World, 95, 89–94, https://doi.org/10.1080/0005772X.2018.1483814, 2018. a

Garibaldi, A., Morales, C., Ashworth, L., Chacoff, N., and Aizen, M.: Los polinizadores en la agricultura, Ciencia Hoy, 21, 34–43, 2012. a

Güngör, E.: Determination of optimun management strategy for honey production forest lands using a'wot and conjoint analysis: a case study in Turkey, Appl. Ecol. Environ. Res., 16, 3437–3459, https://doi.org/10.15666/aeer/1603_34373459, 2018. a

Guzmán-Novoa, E., Eccles, L., Calvete, Y., Mcgowan, J., Kelly, P., and Correa-Benítez, A.: Varroa destructor is the main culprit for the death and reduced populations of overwintered honey bee (Apis mellifera) colonies in Ontario, Canada, Apidologie, 41, 443–450, https://doi.org/10.1051/apido/2009076, 2010. a

Hardstone, M. and Scott, J.: Is Apis mellifera more sensitive to insecticides than other insects?, Pest Manage. Sci., 66, 1171–1180, https://doi.org/10.1002/ps.2001, 2010. a

Hogan, S. J., Homer, M. E., Jeffrey, M. R., and Szalai, R.: Piecewise smooth dynamical systems theory: the case of the missing boundary equilibrium bifurcations, J. Nonlin. Sci., 26, 1161–1173, 2016. a

Jatulan, E., Rabajante, J., Banaay, C., Fajardo, A., and Jose, E.: A mathematical model of intra-colony spread of American foulbrood in European honeybees (Apis mellifera), PloS One, 10, 1–13, https://doi.org/10.1371/journal.pone.0143805, 2015. a

Khoury, D., Myerscough, M., and Barron, A.: A quantitative model of honey bee colony population dynamics, PloS One, 6, 1–6, https://doi.org/10.1371/journal.pone.0018491, 2011. a

Khoury, D., Barron, A., and Myerscough, M.: Modelling food and population dynamics in honey bee colonies, PloS One, 8, 1–7, https://doi.org/10.1371/journal.pone.0059084, 2013. a, b, c

Klein, A., Vaissiere, B., Cane, J., Steffan-Dewenter, I., Cunningham, S., Kremen, C., and Tscharntke, T.: Importance of pollinators in changing landscapes for world crops, P. Roy. Soc. B, 274, 303–313, https://doi.org/10.1098/rspb.2006.3721, 2007. a

Martin, S., Highfield, A., Brettell, L., Villalobos, E., Budge, G., Powell, M., Nikaido, S., and Schroeder, D.: Global honey bee viral landscape altered by a parasitic mite, Science, 336, 1304–1306, https://doi.org/10.1126/science.1220941, 2012. a

Martínez, E., Arroyo, H., Aguilar, N., García, J., Santoyo, V., and Aguilar, J.: Dinámica de adopción de buenas prácticas de producción de miel en la península de Yucatán, México, Revista mexicana de ciencias pecuarias, 9, 48–67, https://doi.org/10.22319/rmcp.v9i1.4366, 2018. a

Mattila, H. and Seeley, T.: Genetic diversity in honey bee colonies enhances productivity and fitness, Science, 317, 362–364, https://doi.org/10.1126/science.1143046, 2007. a

Nates-Parra, G.: Genética del comportamiento: abejas como modelo, Acta Biológica Colombiana, 16, 213–229, 2011. a

Nates-Parra, G.: Iniciativa Colombiana de Polinizadores: Capítulo Abejas, ICPA, http://hdl.handle.net/11438/8800 (last access: 25 February 2022), 2016. a, b

National Agricultural Statistics Service: Honey production, USDA [data set], https://usda.library.cornell.edu/concern/publications/hd76s004z?locale=en&page=3#release-items, last access: August 2021. a

Patel, V., Pauli, N., Biggs, E., Barbour, L., and Boruff, B.: Why bees are critical for achieving sustainable development, Ambio, 50, 49–59, https://doi.org/10.1007/s13280-020-01333-9, 2021. a

Quinn, G. and Keough, M.: Experimental Design and Data Analysis for Biologists, Cambridge University Press, https://doi.org/10.1017/CBO9780511806384, 2002. a

Ratnieks, F. and Carreck, N.: Clarity on honey bee collapse?, Science, 327, 152–153, https://doi.org/10.1126/science.1185563, 2010. a

R Team: R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, https://www.R-project.org (last access: 12 July 2021), 2019. a

Romero-Leiton, J. P.: jpatirom3/Honey-bees-modeling-, GitHub [code], https://github.com/jpatirom3/Honey-bees-modeling-, last access: March 2022. a

Russell, S., Barron, A., and Harris, D.: Dynamic modelling of honey bee (Apis mellifera) colony growth and failure, Ecol. Model., 265, 158–169, https://doi.org/10.1016/j.ecolmodel.2013.06.005, 2013. a, b, c, d, e

Schmickl, T. and Crailsheim, K.: HoPoMo: A model of honeybee intracolonial population dynamics and resource management, Ecol. Model., 204, 219–245, https://doi.org/10.1016/j.ecolmodel.2007.01.001, 2007.  a, b

Schmickl, T. and Crailsheim, K.: TaskSelSim: a model of the self-organization of the division of labour in honeybees, Math. Comput. Model. Dynam. Syst., 14, 101–125, https://doi.org/10.1080/13873950701846662, 2008. a

Stirling, W.: Iteratively reweighted least squares for models with a linear part, J. Roy. Stat. Soc. Ser. C, 33, 7–17, https://doi.org/10.2307/2347657, 1984. a

Torres, D., Ricoy, U., and Roybal, S.: Modeling honey bee populations, PloS One, 10, 1–28, https://doi.org/10.1371/journal.pone.0130966, 2015. a

Ulloa, J., Mondragón, P., Rodríguez, R., Reséndis, J., and Rosas, P.: La miel de abeja y su importancia, Revista Fuente, 2, 11–18, 2010. a

Vit, P.: Productos de la colmena recolectados y procesados por las abejas: Miel, polen y propóleos, Revista del Instituto Nacional de Higiene Rafael Rangel, 35, 32–39, 2004. a

Winfree, R.: The conservation and restoration of wild bees, Ann. N.Y. Acad. Sci., 1195, 169–197, https://doi.org/10.1111/j.1749-6632.2010.05449.x, 2010. a

Winfree, R., Aguilar, R., Vázquez, D., LeBuhn, G., and Aizen, M.: A meta-analysis of bees' responses to anthropogenic disturbance, Ecology, 90, 2068–2076, https://doi.org/10.1890/08-1245.1, 2009. a

Woodard, S., Fischman, B., Venkat, A., Hudson, M., Varala, K., Cameron, S., Clark, A., and Robinson, G.: Genes involved in convergent evolution of eusociality in bees, P. Natl. Acad. Sci. USA, 108, 7472–7477, https://doi.org/10.1073/pnas.1103457108, 2011. a

Wooldridge, J.: Introductory econometrics: A modern approach, Cengage learning, J. Appl. Econometr., 17, 191–193, 2015. a