the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An approach to the modeling of honey bee colonies
Jhoana P. RomeroLeiton
Alejandro Gutierrez
Ivan Felipe Benavides
Oscar E. Molina
Alejandra Pulgarín
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 nonsmooth model where the mortality rate of bees is included as a function of the number of adult bees in the population.
Bees are one of the most ecologically and commercially important insects in the world (Brown and Paxton, 2009). There are approximately 20 000 species of bees (NatesParra, 2011), from which Apis mellifera is the most common for commercial purposes (Hardstone and Scott, 2010). 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 (Patel et al., 2021). 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 (García et al., 2016; Klein et al., 2007; Garibaldi et al., 2012).
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 (NatesParra, 2016). This makes bees important regulators of food production, forest balance and microclimate dynamics (Winfree, 2010). From all species of bees, only a few are eusocial; Apis mellifera is an example of eusocial behavior (Woodard et al., 2011). This species forms colonies so that the survival, reproduction and honey production depend on the structure and size of the colonies (Mattila and Seeley, 2007).
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 zooagricultural goods for commercial trade in the world (Campos and Leyva, 2018; Ulloa et al., 2010; Farouk et al., 2014). The utility of honey is diverse. For example, it works as an antiflu agent (Farouk et al., 2014). 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 (Vit, 2004). 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ía, 2018). 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ör, 2018).
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 (Ratnieks and Carreck, 2010; CoxFoster and VanEngelsdorp, 2009; NatesParra, 2016).
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 (Belić et al., 1986; Schmickl and Crailsheim, 2008; Khoury et al., 2011; Russell et al., 2013; Jatulan et al., 2015; Dennis and Kemp, 2016; Booton et al., 2017). Schmickl and Crailsheim (2007) 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 (Bagheri and Mirzaie, 2019) 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 Khoury et al. (2013). Schmickl and Crailsheim (2007) 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 egglaying rate, nurse bees on larval survival and pollen shortage on cannibalization. Russell et al. (2013) 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 Khoury et al. (2013).
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 nonsmooth model which is analyzed using Filippov's systems theory (Dieci and Lopez, 2009). 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).
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#releaseitems, 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/Honeybeesmodeling (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 (Bretagnolle and Gaba, 2015). Since this kind of temporal trend can inflate the relationship between responses and predictors in statistical models (Wooldridge, 2015), we fitted a model for each of seven 5year subseries (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 subseries 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.
Thus, our modeling approach considered honey production (HP) (pounds) as the response variable for each of the seven 5year subseries; 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 (Quinn and Keough, 2002). Equation (1) shows the specification of this model, where e is an error term with Gaussian distribution and Id is an identity link function.
Estimation of β parameters was performed trough the iterative reweighted least squares method (IWLS) (Stirling, 1984). Different combinations of covariates were fitted in a topdown 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 subseries. β_{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 Team, 2019).
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 R^{2} 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.
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/(T+\mathit{\nu})$, 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 μ_{B}B, 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:
Similarly, we assumed that the number of adult bees decreases naturally. This is represented by term μ_{T}T, 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:
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/(T+u)$, 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:
From Eqs. (2)–(4) we obtain the following system of ODEs:
Model (Eq. 5) can be represented in Fig. 4.
The specifications of parameters involved in the model (Eq. 5) are shown in Table 2.
Let us define $\mathit{\u03f5}=\mathrm{min}\mathit{\{}{\mathit{\mu}}_{\mathrm{B}},{\mathit{\mu}}_{\mathrm{T}}\mathit{\}}$. Then, the set of biological interest of the model (Eq. 5) is given by
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
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}}=(\mathrm{0},\mathrm{0},\mathrm{0})$. Besides, if ${\mathit{\lambda}}^{*}>\mathit{\nu}$ there is a nontrivial equilibrium ${\mathbf{E}}_{\mathrm{1}}=({b}^{*},{\mathit{\lambda}}^{*}\mathit{\nu},{m}^{*})$, 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 righthand side of the system (Eq. 5), which is given by the following Jacobian matrix:
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}}=(\mathrm{0},\mathrm{0},\mathrm{0})$ is locally and asymptotically stable, whereas if ${\mathit{\lambda}}^{*}>\mathit{\nu}$, it is unstable. If ${\mathit{\lambda}}^{*}=\mathit{\nu}$, E_{0} is a nonhyperbolic equilibrium point. Similarly, If ${\mathit{\lambda}}^{*}>\mathit{\nu}$, the coexistence equilibrium ${\mathbf{E}}_{\mathrm{1}}=({b}^{*},{\mathit{\lambda}}^{*}\mathit{\nu},{m}^{*})$ is locally and asymptotically stable, whereas if ${\mathit{\lambda}}^{*}<\mathit{\nu}$, it is unstable. If ${\mathit{\lambda}}^{*}=\mathit{\nu}$, E_{1} is a nonhyperbolic equilibrium point.
A summary of the results from Proposition 3.3 is given in Table 3.
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
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.
Russell et al. (2013)Torres et al. (2015)Martínez et al. (2018)Russell et al. (2013)Russell et al. (2013)Winfree et al. (2009)Khoury et al. (2013)Brodschneider and Crailsheim (2010)In Fig. 5 we present the stability of the equilibrium E_{0}, 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 6 numerically represents the stability of the equilibrium E_{1}. 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.
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 nonsmooth 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 (Dieci and Lopez, 2009) to study the system (Eq. 5), where R_{1} and R_{2} are considered as death thresholds. If the number T of adult bees is big enough from a R_{2}, the death of bees due to a stress factor (σ) will be maximum (σ_{max}). Similarly, from a R_{1}, 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 R_{1} and R_{2}, the death rate due to a stress factor will be stable (σ_{1}). The previous consideration can be mathematically written as follows:
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 (Amador et al., 2021). The limits between the different regions are called discontinuity surfaces, which are usually denoted as Σ_{i} (Hogan et al., 2016).
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 pseudoequilibrium 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/Honeybeesmodeling (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 8 shows the switching surfaces Σ_{1} (left side) and Σ_{2} (right side), which divide the first octant into the following regions: S_{1} (left side of Σ_{1}), S_{2} (region between Σ_{1} and Σ_{2}), and S_{3} (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 S_{2} and S_{3}, 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 breakeven point, as shown by Fig. 7.
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 nonsmooth 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 R_{1} and R_{2} were considered for the stress factor. The Filippov analysis for the nonsmooth model showed that there is always a crossover in the switching surfaces. After 100 weeks, the trajectories switched between regions S_{2} and S_{3}, 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 S_{2}. We have taken into account the fact that the equilibrium point E_{1} is stable, and, therefore, the equilibrium value is reached when t→∞.
For the statistical approach, generalized linear models were built to study seven 5year subseries 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 monocultures (Belsky and Joshi, 2019) or the infestation by the ectoparasitic Varroa mite destructor (GuzmánNovoa et al., 2010; Martin et al., 2012).
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 nonsmooth 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.
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
Thus, $B+T\le \mathit{\beta}/\mathit{\u03f5}$. Similarly, for M we get that
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),
Replacing B defined in Eq. (A1) in the first equation of the system (Eq. 7) we obtain
From the previous equation we have T=0 or T≠0. First, let us suppose that T≠0, thus
where
Replacing T defined on Eq. (A2) into the third equation of the system, Eq. (7), we have that
from where
Now, replacing T defined on Eq. (A2) into the Eq. (A1), we obtain that
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}}=(\mathrm{0},\mathrm{0},\mathrm{0})$.
A2 Local stability analysis
The stability of E_{0} can be found by evaluating the above matrix in E_{0}. Thus,
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:
From the Routh–Hurwitz criterion (DeJesus and Kaufman, 1987), we have that the roots of the above polynomial have negative real parts, if and only if the following conditions are met:
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 E_{1}. Evaluating matrix (Eq. 8) in E_{1} we obtain
An eigenvalue of the previous matrix is $\mathit{\eta}=\mathit{\alpha}\mathit{\delta}({\mathit{\lambda}}^{*}\mathit{\nu})$, which is negative if ${\mathit{\lambda}}^{*}\ge \mathit{\nu}$, while the other two eigenvalues are given by the roots of the following characteristic polynomial:
Similarly, from the Routh–Hurwitz criterion (DeJesus and Kaufman, 1987) the roots of the above polynomial have negative real parts, if and only if the following conditions are met:
It is clear that Θ_{1}>0 while Θ_{2}>0, if and only if ${\mathit{\lambda}}^{*}>\mathit{\nu}.$
We assume that the state space consists of the three regions S_{1}, S_{2} and S_{3}⊂ℝ^{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}}(B,T,M)$, ${F}_{\mathrm{2}}(B,T,M)$ and ${F}_{\mathrm{3}}(B,T,M)$ determined by the value of T according to Eq. (9). If T<R_{1}, we have
When ${R}_{\mathrm{1}}\le T\le {R}_{\mathrm{2}}$, we have that
Finally, at T>R_{2}
These three functions are also assumed to be defined through the state space, even though they are only used in their respective regions (Amador et al., 2013). Thus, the regions S_{1}, S_{2} and S_{3} can be defined by the following sets:
The regions are separated by a switching surface: limit of discontinuity represented by the scalar function ${H}^{\left(i\right)}:{\mathbb{R}}^{\mathrm{3}}\u27f6\mathbb{R}$; ${H}^{\left(i\right)}\left(x\right)=T{R}_{i}$ with $i=\mathrm{1},\mathrm{2}$, so that
As $x=(B,T,M)\in {R}^{\mathrm{3}}$, the gradient of H^{(i)}(x) can be defined as follows:
Thus, taking the regions Σ_{i} in which the gradient does not null in the set ${\mathrm{\Gamma}}_{i}=\mathit{\{}x\in {\mathbb{R}}^{\mathrm{3}}:{H}^{\left(i\right)}(x)=\mathrm{0},\phantom{\rule{0.25em}{0ex}}i=\mathrm{1},\mathrm{2}\mathit{\}}$, for x∈Γ_{1}, we define
and, for x∈Γ_{2},
Filippov's systems theory (Amador et al., 2017) states that at points x, where
With the vector field
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<{\displaystyle \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>{\displaystyle \frac{\left({\mathit{\mu}}_{\mathrm{T}}+{\mathit{\sigma}}_{\mathrm{1}}\right)T}{\mathit{\omega}}}.$$ 
Tangential point when
$$B={\displaystyle \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={\displaystyle \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}}}<B<{\displaystyle \frac{\left({\mathit{\mu}}_{\mathrm{T}}+{\mathit{\sigma}}_{\mathrm{1}}\right)T}{\mathit{\omega}}$$with the equation $\dot{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}}={\displaystyle \frac{\mathit{\omega}}{{\mathit{\sigma}}_{\mathrm{min}}{\mathit{\sigma}}_{\mathrm{1}}}}\left({\displaystyle \frac{B}{T}}\right){\displaystyle \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<{\displaystyle \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>{\displaystyle \frac{\left({\mathit{\mu}}_{\mathrm{T}}+{\mathit{\sigma}}_{\mathrm{max}}\right)T}{\mathit{\omega}}}.$$ 
Tangential point when
$$B={\displaystyle \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={\displaystyle \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}}}<B<{\displaystyle \frac{\left({\mathit{\mu}}_{\mathrm{T}}+{\mathit{\sigma}}_{\mathrm{max}}\right)T}{\mathit{\omega}}},$$with the equation $\dot{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}}={\displaystyle \frac{\mathit{\omega}}{{\mathit{\sigma}}_{\mathrm{1}}{\mathit{\sigma}}_{\mathrm{max}}}}\left({\displaystyle \frac{B}{T}}\right){\displaystyle \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}{\displaystyle \frac{T}{T+u}}\mathit{\alpha}M\mathit{\delta}TM\end{array}\right).$$
The R and Matlab codes and the cleaned data set are available at the following GitHub link: https://github.com/jpatirom3/Honeybeesmodeling (RomeroLeiton, 2022).
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#releaseitems, National Agricultural Statistics Service, 2021).
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 nonsmooth 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.
The contact author has declared that neither they nor their coauthors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Jhoana P. RomeroLeiton 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.
This paper was edited by John M. Halley and reviewed by Sergio Estay and one anonymous referee.
Amador, J., Granada, H., Redondo, J., and Olivar, G.: Nonlinear and Nonsmooth Dynamics in StressSickness 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/s1259101201382, 2013. a
Amador, J. A., Redondo, J. M., OlivarTost, G., and Erazo, C.: CooperationBased 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.: Stressmediated 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/s1359301503025, 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
CoxFoster, D. and VanEngelsdorp, D.: Saving the honeybee, Scient. Am., 300, 40–47, https://doi.org/10.1038/scientificamerican040940, 2009. a
DeJesus, E. and Kaufman, C.: RouthHurwitz 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/S071834292016000300008, 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ánNovoa, E., Eccles, L., Calvete, Y., Mcgowan, J., Kelly, P., and CorreaBení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 intracolony 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., SteffanDewenter, 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
NatesParra, G.: Genética del comportamiento: abejas como modelo, Acta Biológica Colombiana, 16, 213–229, 2011. a
NatesParra, 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#releaseitems, 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/s13280020013339, 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.Rproject.org (last access: 12 July 2021), 2019. a
RomeroLeiton, J. P.: jpatirom3/Honeybeesmodeling, GitHub [code], https://github.com/jpatirom3/Honeybeesmodeling, 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 selforganization 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.17496632.2010.05449.x, 2010. a
Winfree, R., Aguilar, R., Vázquez, D., LeBuhn, G., and Aizen, M.: A metaanalysis of bees' responses to anthropogenic disturbance, Ecology, 90, 2068–2076, https://doi.org/10.1890/081245.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
 Abstract
 Introduction
 The statistical model approach
 The smooth mathematical model approach
 The nonsmooth mathematical model approach
 Discussion
 Appendix A
 Appendix B
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References
 Abstract
 Introduction
 The statistical model approach
 The smooth mathematical model approach
 The nonsmooth mathematical model approach
 Discussion
 Appendix A
 Appendix B
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References