Next Article in Journal
Stability Analysis of Simple Root Seeker for Nonlinear Equation
Next Article in Special Issue
An Inertial Subgradient Extragradient Method for Approximating Solutions to Equilibrium Problems in Hadamard Manifolds
Previous Article in Journal
Turbulent Free Convection and Thermal Radiation in an Air-Filled Cabinet with Partition on the Bottom Wall
Previous Article in Special Issue
Two Convergence Results for Inexact Infinite Products of Non-Expansive Mappings
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Parameter Estimation Analysis in a Model of Honey Production

by
Atanas Z. Atanasov
1,†,
Slavi G. Georgiev
2,3,*,‡ and
Lubin G. Vulkov
3,‡
1
Department of Agricultural Machinery, Agrarian Industrial Faculty, University of Ruse, 7004 Ruse, Bulgaria
2
Department of Informational Modeling, Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria
3
Department of Applied Mathematics and Statistics, Faculty of Natural Sciences and Education, University of Ruse, 7004 Ruse, Bulgaria
*
Author to whom correspondence should be addressed.
Current address: Department of Agricultural Machinery, University of Ruse, 8 Studentska Str., 7004 Ruse, Bulgaria.
Current address: Department of Applied Mathematics and Statistics, University of Ruse, 8 Studentska Str., 7004 Ruse, Bulgaria.
Axioms 2023, 12(2), 214; https://doi.org/10.3390/axioms12020214
Submission received: 30 December 2022 / Revised: 7 February 2023 / Accepted: 8 February 2023 / Published: 17 February 2023

Abstract

:
Honeybee losses are an extensive global problem. In this study, a new compartment model of honeybee population that mainly concerns honey production is developed. The model describes the interaction of the food stock with the brood (immature bees), adult bees and produced honey. In the present paper, the issue of an adequate model recovery is addressed and the parameter identification inverse problem is solved. An adjoint equation procedure to obtain the unknown parameter values by minimizing the functional error during a period of time is proposed. Numerical simulations with realistic data are discussed.

1. Introduction

Honeybee colonies are important for agriculture and the environment. They help plant reproduction by pollination, while beekeeping redounds to the development of rural areas. Unfortunately, in recent decades, a ubiquitous decline in both managed and unmanaged colonies has been observed. This is a global problem, since the bees contribute to the ecological equilibrium. If the bee population shrinks or disappears, plants would not get pollinated and would die off. Then, herbivorous animals would not have food and would go extinct, and they would be followed by carnivorous animals, including humans. Thus, preventing bee colonies from losses is of a paramount importance for preserving live on Earth in general.
Of the many species of bees, only a small number of them are eusocial; Apis mellifera is an example of eusocial behavior [1]. This species form colonies thus the survival, reproduction and honey production are directly dependent on the size and the structure of the colonies [2].
Honey, produced by honeybees, is a sweet natural substance, derived greatly from the nectar of flowers and transformed by a group of enzymes, which are present in the saliva of the worker bees. The honey is also airy and evaporates by its filtering, and is eventually stored inside the hives. Honey from Apis mellifera is one of the most essential zoo-agricultural goods for commercial trade in the world [3]. Regarding the honey trade, the USA is the global leader in imports. Concerning production, China is the global leader, following by Turkey, Iran, Ukraine and the Russian Federation. Finally, with respect to quality, Bulgarian honey is the most pure and sweet [4].
Beekeepers produce a variety of agriculture products, in addition to honey, including royal jelly, propolis and beeswax. This paper aims to develop mathematical modeling of the honeybee population dynamics and, therefore, honey production.
The most fundamental honeybee population model is suggested by [5], where only two compartments are explored—the young hive bees and the matured forager bees. This model is extended in [6,7], where the brood, the age of the foragers and the food are also included in the studies, accounting for the delay of maturing. Such investigation is done in [8], where a different form of the recruitment rate is used. In the study [9], exogenous stress is assumed to impact the recruitment process, social inhibition and the queen laying rate, causing a potential colony decline.
There are models developed as an effort to understand the decrease in colony numbers in recent decades. A survey in the USA suggests that treating against disease and mite infestation in the right way lowers the chance of colony loss [10]. Extensive study of the transition from hive to forager bees is performed in [11]. A comparison between the losses in different parts of the world is performed in [12].
The mysterious disease, whose causal factors are not entirely agreed on, is called Colony Collapse Disorder (CCD). It is characterized by rapid loss of forager bees but absence of dead bodies near the hive, lack of pest and mite invasion of the hive, and bees’ reluctance to consume food provided by the beekeeper. The first recorded massive colony loss is described as the ‘Isle of Wight Disease’ [13]. The effect of protein sources has been proposed as a potential cause for collapse [14]. A special CCD model is designed in [15], where the contagious adult bees are isolated from the others. A review of the suspected causal factors for the colony declines is summarized in [16].
Other models focus on particular parts of the surrounding environment such as food availability [17], age structure [18], seasonal effects [19], Varroa mites [20] and others [21,22], including the model memory property [23].
In [24], populations of adult and immature (brood) honeybees as well as their honeybee production are examined via mathematical and statistical modeling approaches. It is shown that, if a bee population is exposed to a stress factor (i.e., habitat destruction, Varroamites, climate variability, heavy metals, etc.), the number of individuals declines over time as well as the produced honey. The complex issue of the sustainability of honeybee colonies is important not only for the survival of the species but also for food security and the overall health of the environment. To ensure the sustainability of honeybee colonies, it is important to take measures such as providing adequate habitats, reducing pesticide exposure and promoting disease management practices. Aiming at the latter, the sophisticated processes of population dynamics have to be investigated via mathematical modeling.
In the present work we study the relationship between the population size of honeybees (Apis mellifera) and honey production if the bee colony is exposed to a number of stress factors that exogenously cause the death of individuals and therefore a possible reduction in honey production. Here lies the main originality of the study—suggesting a novel model for encountering the interaction between the bee castes and the amount of honey, stored in the hive.
Furthermore, in the investigation the inverse problem of identifying the food and honey consumption rates by the immature and adult bees is solved as well as the brood maturation rate. These quantities are of extreme importance for understanding the complex dynamics of the hive. It is done via the adjoint equations optimization approach. Such a study is performed in [25], where the contaminated bees are modeled as a separate compartment. Similar investigation is done in [26] but, for the coefficient identification, a trust-region reflexive algorithm is used.
This paper is organized as follows. In the next section, we extend the mathematical model, studied in [24], taking into account the food stock. What is more, we study the existence and non-negativity of the solutions. Section 3 is devoted to the parameter estimation analysis of the model. Section 4 is dedicated to numerical experiments regarding the direct and inverse problems. The paper is concluded in Section 5.

2. Mathematical Model

In this section, we introduce a mathematical model that explains the interaction between the food stock and the brood (immature bees), adult bees and the amount of produced honey.
Following the results in [24,27] we establish a mathematical model that presents the interactions among brood B ( t ) at time t, adult bees A ( t ) and the amount of honey production M ( t ) , taking into account the weight of food stock F ( t ) .
We assume that the brood grows at a rate β , proportionally to the number of adult bees. This is given by term A / ( A + ν ) , where ν is the mean saturation rate (number of adult bees required for immature bees to achieve half of their maximal number). The number of bees surviving to the adult stage influences the number of immature bees.
The latter is modeled by the term ω B , where ω denotes the maturation rate to adult stage, and 1 / ω indicates the time spent before achieving the adult stage. The number of immature bees is decreased by natural death and it is modeled by the term μ B B , where μ B denotes the natural mortality rate of the immature stage. Following this discussion and those in [24,27] we consider the following system of ODEs:
d F d t = c A γ B ,
d B d t = β A A + ν ω B μ B B ,
d A d t = ω B μ A A σ A .
d M d t = ρ A A + u α M δ A M .
The model (1)–(4) is illustrated in the diagram of Figure 1.
It is assumed in the derivation of Equation (3) that the number of adult bees diminish naturally and it is demonstrated by the term μ A A , where μ A is the natural mortality rate of the adult stage. However, the bees can also die because of a stress factor. This is represented by the term σ A , where σ is the death rate due to a stressor (climate change, loss of habitat, heavy metals or pesticides, poor beekeeper’s management, etc.) acting on bees at the adult stage.
Equation (4) shows that the production of honey in hives increase at a rate ρ , which is influenced by the number of adult bees, given by the term A / ( A + u ) , where u is the mean saturation rate.
One important cause for decreasing of the honey is the feeding of immature bees, which is demonstrated by the term α M , where α is the honey loss rate.
The term δ A M represents the loss of honey production because of the consumption of adult bees, where δ is the adult bees’ honey consumption rate.
For more details on the specifications of the parameters in the model we refer to Table 2 in [24].
We solve the system of ordinary differential Equations (1)–(4) with initial conditions
F ( 0 ) = F 0 0 , B ( 0 ) = B 0 0 , A ( 0 ) = A 0 0 , M ( 0 ) = M 0 0 .
Using Theorem 7.1 in [28] one could easily prove that the subsystem (2)–(4) is positive (short for “non-negativity preserving”) in the sense that, if
B ( 0 ) 0 , A ( 0 ) 0 , M ( 0 ) 0 ,
then
B ( t ) 0 , A ( t ) 0 , M ( t ) 0 , t 0 .
This property is biologically relevant to the model.

3. Parameter Identification

In this section, the parameter inverse problem is defined. Such problems appear very often in practice. The problem (1)–(5), where the values of the parameters are known, is well-posed and it is called a direct problem. However, in the real world, the values of some of the coefficients are not directly measurable but they are very important for professional honeybee management. Their reconstruction, provided that additional information is given, is referred to an inverse problem. Inverse problems are ill-posed and harder to solve. We employ the adjoint equation optimization approach [29,30].
The parameters to be reconstructed are p = ( p 1 , p 2 , p 3 , p 4 , p 5 ) , p 1 = α , p 2 = γ , p 3 = δ , p 4 = σ , p 5 = ω , and
p S adm = p R 5 : 0 < p i < P i , i = 1 , , 5 .
The admissible set S adm is defined by the biology of the honeybee [31]. To find the parameters p , though, some new information must be brought. In many cases it is possible to measure the model functions at some discrete times. In reality, counting the brood B is a difficult task, so we adopt measurements of the functions
F obs ( t k ) = X k , A obs ( t k ) = Y k , M obs ( t k ) = Z k
for k = 1 , , K . We assume all functions are measured at some predefined time instances. The observation times for every function may be different.
In practice, the observations are obtained from electronic devices equipping the hive. In a quasi-real setting, first the direct problem is solved and then the observations are extracted from the solution to the direct problem.
To solve the inverse problem, the least-square function
Φ ( p ) = Φ ( α , γ , δ , σ , ω ) = Φ F ( α , γ , δ , σ , ω ) + Φ A ( α , γ , δ , σ , ω ) + Φ M ( α , γ , δ , σ , ω ) = k = 1 K ( F ( t k ; p ) X k ) 2 + k = 1 K ( A ( t k ; p ) Y k ) 2 + k = 1 K ( M ( t k ; p ) Z k ) 2
is minimized, e.g., by a gradient method [32], where Ψ ( t k ; p ) , Ψ { F , A , M } are the theoretical quantities from the model and Ξ k , Ξ { X , Y , Z } are the observed values in practice.
Now we state an expression for the gradient of the function Φ p : = Φ ( p ) .
Theorem 1.
The gradient Φ p ( Φ α , Φ γ , Φ δ , Φ σ , Φ ω ) is given by
Φ α = 0 T φ M ( t ) M ( t ) d t ,
Φ γ = 0 T φ F ( t ) B ( t ) d t ,
Φ δ = 0 T φ M ( t ) A ( t ) M ( t ) d t ,
Φ σ = 0 T φ A ( t ) A ( t ) d t ,
Φ ω = 0 T φ A ( t ) B ( t ) d t ,
where the triple { φ M , φ F , φ A } is the unique solution of the adjoint system
d φ F d t = 2 k = 1 K ( F X ) δ ( t t k ) ,
d φ A d t = c φ F + ( μ A + σ ) φ A + δ · M ρ u ( A + u ) 2 φ M + 2 k = 1 K ( A Y ) δ ( t t k ) ,
d φ M d t = ( α + δ · A ) φ M + 2 k = 1 K ( M Z ) δ ( t t k ) ,
φ F ( T ) = φ A ( T ) = φ M ( T ) = 0 .
Proof. 
We denote δ p = ( δ α , δ γ , δ δ , δ σ , δ ω ) and δ α = ε h 1 , δ γ = ε h 2 , δ δ = ε h 3 , δ σ = ε h 4 , δ ω = ε h 5 .
If δ F ( t ; p ) = F ( t ; p + δ p ) F ( t ; p ) , δ A ( t ; p ) = A ( t ; p + δ p ) A ( t ; p ) and δ M ( t ; p ) = M ( t ; p + δ p ) M ( t ; p ) , write the ODE system for F ( t ; p + δ p ) , A ( t ; p + δ p ) and M ( t ; p + δ p ) as (1), (3) and (4) with initial conditions F 0 , A 0 and M 0 (5).
Then, calculate the differences of the corresponding equations to obtain an ODE system for δ F , δ A and δ M with zero initial conditions.
d d t δ F = c δ A δ γ B ,
d d t δ A = ( μ A + σ ) δ A δ σ A + δ w B ,
d d t δ M = ρ u δ A ( A + u ) 2 δ · M δ A ( α + δ · A ) δ M δ α M δ δ A M .
We find the increment of the functional Φ ( p ) :
Φ ( p + δ p ) Φ ( p ) = 2 k = 1 K δ F ( t k ; p ) F ( t k ; p ) X k + 2 k = 1 K δ A ( t k ; p ) A ( t k ; p ) Y k + 2 k = 1 K δ M ( t k ; p ) M ( t k ; p ) Z k = 2 k = 1 K 0 T δ F ( t k ; p ) F ( t k ; p ) X k δ ( t t k ) d t + 2 k = 1 K 0 T δ A ( t k ; p ) A ( t k ; p ) Y k δ ( t t k ) d t + 2 k = 1 K 0 T δ M ( t k ; p ) M ( t k ; p ) Z k δ ( t t k ) d t .
Let us multiply Equations (18)–(20) by smooth functions φ F ( t ) , φ A ( t ) and φ M ( t ) s.t. φ F ( T ) = φ A ( T ) = φ M ( T ) = 0 and integrate both sides of the results from 0 to T:
0 T φ F d d t δ F + φ A d d t δ A + φ M d d t δ M d t = c 0 T φ F δ A d t δ γ 0 T φ F B d t ( μ A + σ ) 0 T φ A δ A d t δ σ 0 T φ A A d t + δ w 0 T φ A B d t + ρ u 0 T φ M δ A ( A + u ) 2 d t δ 0 T φ M M δ A d t ( α + δ · A ) 0 T φ M δ M d t δ α 0 T φ M M d t δ δ 0 T φ M A M d t .
On the other hand, integrating by parts and using the facts that φ F ( T ) = φ A ( T ) = φ M ( T ) = 0 and δ F ( 0 ) = δ A ( 0 ) = δ M ( 0 ) = 0 , we obtain
0 T φ F d d t δ F d t + 0 T φ A d d t δ A d t + 0 T φ M d d t δ M d t = 0 T δ F d φ F d t d t 0 T δ A d φ A d t d t 0 T δ M d φ M d t d t .
Let us place the expressions for d φ F d t , d φ A d t and d φ M d t from (14)–(16) in (22):
0 T φ F d d t δ F + φ A d d t δ A + φ M d d t δ M d t = c 0 T φ F δ A d t ( μ A + σ ) 0 T φ A δ A d t δ 0 T φ M M δ A d t + ρ u 0 T φ M 1 ( A + u ) 2 δ A d t ( α + δ · A ) 0 T φ M δ M d t 2 0 T δ F k = 1 K ( F X ) δ ( t t k ) d t 2 0 T δ A k = 1 K ( A Y ) δ ( t t k ) d t 2 0 T δ M k = 1 K ( M Z ) δ ( t t k ) d t .
Equating (21) and (23) yields
2 0 T δ F k = 1 K ( F X ) δ ( t t k ) d t + 2 0 T δ A k = 1 K ( A Y ) δ ( t t k ) d t + 2 0 T δ M k = 1 K ( M Z ) δ ( t t k ) d t = δ α 0 T φ M M d t + δ γ 0 T φ F B d t + δ δ 0 T φ M A M d t + δ σ 0 T φ A A d t δ ω 0 T φ A B d t .
Rewriting the last expression give
Φ ( α + ε h 1 , γ + ε h 2 , δ + ε h 3 , σ + ε h 4 , w + ε h 5 ) Φ ( α , γ , δ , σ , ω ) = h 1 0 τ φ M M d t + h 2 0 τ φ F B d t + h 3 0 τ φ M A M d t + h 4 0 τ φ A A d t h 5 0 τ φ A B d t ε .
Now, taking h 2 = h 3 = h 4 = h 5 = 0 , dividing both sides by ε h 1 and taking the limit ε 0 we find the formula for Φ α in the theorem.
Analogously, we obtain the formulae for Φ γ , Φ δ , Φ σ and Φ ω (10)–(13). □
Employing the fundamental property of the Dirac-delta function 0 T f ( t ) δ ( t t k ) d t = f ( t k ) , t k ( 0 , T ) , where f ( t ) is a continuous function, (14)–(17) could be rewritten in its equivalent form:
d φ F d t = 0 , t t k , k = 1 , , K , d φ A d t = c φ F + ( μ A + σ ) φ A + δ · M ρ u ( A + u ) 2 φ M , t t k , k = 1 , , K , d φ M d t = ( α + δ · A ) φ M , t t k , k = 1 , , K , φ F t = t k = 2 F ( t k ; p ) X k , k = 1 , , K , φ A t = t k = 2 A ( t k ; p ) Y k , k = 1 , , K , φ M t = t k = 2 M ( t k ; p ) Z k , k = 1 , , K , φ F ( T ) = φ A ( T ) = φ M ( T ) = 0 .
Having obtained the gradient, we employ an iterative procedure as follows, where the new approximation p s + 1 is defined by
p s + 1 = p s r Φ ( p s ) ,
where r R 5 + are gradient multipliers. The iterations start at chosen p 0 and end if p s : = p s + 1 p s < ε p , where ε p is a tolerance quantity, else increase s : = s + 1 and start a new iteration. The final approximation is denoted with p ˇ and it is called a nonlinear estimator.

4. Numerical Experiments

This section is devoted to presenting numerical tests which demonstrate the algorithm application. Firstly, the numerical algorithm is summarized. Then, the direct problem is solved and its solution is used to obtain measurements for the inverse problem.

4.1. Numerical Procedure

All the programming code is implemented in the MATLAB® environment. For solving the ODE systems (1)–(5) and (14)–(17), a Runge–Kutta-type method is used. The algorithm for solving the inverse problem could be described as follows:
  • Choose initial approximation p 0 .
  • Set s : = 0 .
  • Until p s < ε p do
    3.1.
    Solve system (1)–(5) with p s to obtain F, B, A and M.
    3.2.
    Solve system (14)–(17) to obtain φ F , φ A and φ M .
    3.3.
    Compute the gradient Φ p (9)–(13).
    3.4.
    Calculate p s + 1 by (24) and set s : = s + 1 .
  • The estimator is set to p ˇ : = p s .

4.2. Direct Problem

Let us first solve the direct problem (1)–(5) with realistic data given in [24,27]. The adult food collection rate is assumed to be c = 0.04 g/bee/day. The larval consumption rate is γ = 0.12 g/bee/day. The brood reproduction rate is β = 0.92 bee/day. The adult maturation rate is ω = 0.95 day−1. The brood natural mortality rate is μ B = 0.11 day−1. The adult bee natural mortality rate is μ A = 0.29 day−1. The adult bee stressor mortality rate is σ = 0.1 day−1. The honey production rate is ρ = 0.23 bees/day. The rate of natural honey loss is α = 0.018 g/day. The honey consumption rate is δ = 0.571 g/bee/day. The half saturation rates are ν = u = 1 thousand bees.
We simulate the hive development for a typical foraging season, lasting T = 100 days. At the beginning of the season, there are F 0 = 10 kilograms of food stores, B 0 = 2000 larvae, A 0 = 10,000 adult bees and M 0 = 1 kilogram honey. The outcome is plotted in Figure 2.
It could be observed that the hive approaches its equilibrium state relatively fast. It is characterized by a small amount of honey as well as a small number of larvae and adult bees. This is approved by the phase space diagram for a fixed F 0 (Figure 3), which shows no dependence on the initial conditions. Only in case of B 0 = A 0 = 0 , then the extinction equilibrium is approached.
Of course, it is not always true. If there is a hazard present in the environment, i.e., the stress death rate is as high as σ = 0.5 , then the extinction equilibrium is the only attractor, see Figure 4. This unarguably means that the hive would eventually collapse unless something is drastically changed.

4.3. Inverse Problem

Let us solve the inverse problem of identifying the parameters (6) p = ( α , γ , δ , σ , ω ) = ( 0.018 , 0.12 , 0.571 , 0.1 , 0.95 ) . The values of the other parameters and initial conditions remain the same as in the direct problem setting.
We define K = 19 equidistantly distributed observations of type (7), i.e., one observation in every 5 days. The admissible set is set to S adm ( 0 , 1 ) 5 . The values r are tuned empirically and they are given in Table 1.
The respective values (8) are Φ F ( p ˇ ) = 0.1815 , Φ A ( p ˇ ) = 3.7902 and Φ M ( p ˇ ) = 0.3807 . The parameters are recovered with moderate precision, but the honeybee dynamics are reconstructed in an accurate manner. The root mean squared errors are small as RMSE F = 0.0977 , RMSE A = 0.4466 and RMSE M = 0.1416 .
Finally, we perform a test with perturbed measurements to explore the impact of the observation error on the parameter identification. Every electronic device has its instrumental error, so testing with noisy observation is meaningful. We add Gaussian noise to the observations (7), in particular the error in a single observation is not greater than 1 % with 95 % confidence. The results, following the same steps, are given in Table 2.
The outcomes are similar as the functional values are Φ F ( p ˇ ) = 0.2005 , Φ A ( p ˇ ) = 3.6929 and Φ M ( p ˇ ) = 0.3403 . The root mean squared errors are again small as RMSE F = 0.1027 , RMSE A = 0.4409 and RMSE M = 0.1338 . All these demonstrate the robustness and the applicability of the suggested approach with realistic data.

5. Conclusions

Honeybees are one of the most important species on Earth. Their steady colony number decline is a major global problem. To fight this issue, professional honeybee management must take well-designed precautionary measures. The obtained results in this study help beekeepers to foresee the forward colony dynamics. It is crucial to have the ability to simulate the future development and it is here where mathematical modeling comes to the rescue. Then, adequate measures could be undertaken in order to prevent or to revert a colony collapse.
The novelty of the paper is twofold. To begin with, we proposed a new mathematical approach for modeling of honeybee colonies. We analyzed populations of immature and adult bees as well as their honey production. In the context of honeybee colony dynamics, we model the interaction between the different compartments, focusing on parameter recovery. Secondly, the defined ill-posed problem is solved by means of the adjoint equation optimization method. The reconstructed parameters are unobservable in reality but vital for the colony population dynamics. The computational examples with realistic data demonstrate how to apply the approach in practice.
There are many ways to further develop this research. The considered model could be extended to account for mites, viruses and other hazards. Temperature and seasonal effects also worth considering. What is more, activating the hereditary property via fractional-order derivatives almost always results in a better fit. A broader qualitative analysis to better understand the complex phenomena, processing in the hive, is on the agenda as well.

Author Contributions

Conceptualization, A.Z.A. and L.G.V.; methodology, S.G.G. and L.G.V.; software, S.G.G.; validation, S.G.G.; formal analysis, L.G.V.; investigation, S.G.G. and L.G.V.; resources, A.Z.A. and S.G.G.; data curation, A.Z.A. and S.G.G.; writing—original draft preparation, S.G.G.; writing—review and editing, S.G.G. and L.G.V.; visualization, S.G.G.; supervision, L.G.V.; project administration, A.Z.A.; funding acquisition, A.Z.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Bulgarian National Science Fund under Project KP-06-PN 46-7 “Design and research of fundamental technologies and methods for precision apiculture”.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Woodard, S.; Fischman, B.; Venkat, A.; Hudson, M.; Varala, K.; Cameron, S.; Clark, A.; Robinson, G. Genes involved in convergent evolution of eusociality in bees. Proc. Natl. Acad. Sci. USA 2011, 108, 7472–7477. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Matilla, H.; Seeley, T. Genetic diversity in honey bee colonies enhances productivity and fitness. Science 2007, 317, 362–364. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Farouk, K.; Palmera, K.; Sepúlveda, P. Abejas. In InfoZoa Boletín de Zoología; Universidad del Magdalena: Santa Marta, Colombia, 2014; Volume 6. [Google Scholar]
  4. Bulgarian Honey. 2023. Available online: https://www.bulgarianhoney.com/quality.htm (accessed on 7 February 2023).
  5. Khoury, D.S.; Myerscough, M.R.; Barron, A.B. A quantitative model of honey bee colony population dynamics. PLoS ONE 2011, 6, e18491. [Google Scholar] [CrossRef] [Green Version]
  6. Khoury, D.S.; Barron, A.B.; Meyerscough, M.R. Modelling food and population dynamics honey bee colonies. PLoS ONE 2013, 8, e0059084. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Meyerscough, M.R.; Khoury, D.S.; Ronzani, S.; Barron, A.B. Why do hives die? Using mathematics to solve the problem of honey bee colony collapse. In The Role and Importance of Mathematics in Innovation: Proceedings of the Forum “Math-for-Industry”; Anderssen, B., Ed.; Springer: Singapore, 2017; Volume 25, pp. 35–50. [Google Scholar]
  8. Russel, S.; Barron, A.B.; Harris, D. Dynamics modelling of honeybee (Apis mellifera) colony growth and failure. Ecol. Model. 2013, 265, 138–169. [Google Scholar]
  9. Booton, R.D.; Iwasa, Y.; Marshall, J.A.R.; Childs, D.Z. Stress-mediated Allee effects can cause the sudden collapse of honey bee colonies. J. Theor. Biol. 2017, 420, 213–219. [Google Scholar] [CrossRef] [Green Version]
  10. Finley, J.; Camazine, S.; Frazier, M. The epidemic of honey bee colony losses during the 1995–1996 season. Am. Bee J. 1996, 136, 805–808. [Google Scholar]
  11. Amdam, G.V.; Omholt, S.W. The hive bee to forager transition in honeybee colonies: The double repressor hypothesis. J. Theor. Biol. 2003, 223, 451–464. [Google Scholar] [CrossRef]
  12. van der Zee, R.; Pisa, L.; Andronov, S.; Brodschneider, R.; Charriere, J.D.; Chlebo, R.; Coffey, M.F.; Cralisheim, K.; Dahle, B.; Gajda, A.; et al. Managed honey bee colony losses in Canada, China, Europe, Israel and Turkey for the winters of 2008–2009 and 2009–2010. J. Apic. Res. 2012, 51, 100–114. [Google Scholar] [CrossRef] [Green Version]
  13. Bailey, L. The ‘Isle of Wight Disease’: The Origin and Significance of the Myth. Bee World 1964, 45, 32–37. [Google Scholar] [CrossRef]
  14. Kulincevic, J.M.; Rothenbuhler, W.C.; Rinderer, T.E. Disappearing disease. Part 1—Effects of certain protein sources given to honey bee colonies in Florida. Am. Bee J. 1982, 122, 189–191. [Google Scholar]
  15. Dornberger, L.; Mitchell, C.; Hull, B.; Ventura, W.; Shopp, H.; Kribs-Zaleta, C.; Kojouharov, H.; Grover, J. Death of the Bees: A Mathematical Model of Colony Collapse Disorder; Technical Report 2012-12, Mathematics Preprint Series; University of Texas at Arlington Mathematics Department: Arlington, TX, USA, 2012. [Google Scholar]
  16. Hristov, P.; Shumkova, R.; Palova, N.; Neov, B. Factors associated with honey bee colony losses: A mini-review. Vet. Sci. 2020, 7, 166. [Google Scholar] [CrossRef] [PubMed]
  17. Bagheri, S.; Mirzaie, M. A mathematical model of honey bee colony dynamics to predict the effect of pollen on colony failure. PLoS ONE 2019, 14, e0225632. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Betti, M.I.; Wahl, L.M.; Zamir, M. Reproduction number and asymptotic stability for the dynamics of a honey bee colony with continuous age structure. Bull. Math. Biol. 2017, 79, 1586–1611. [Google Scholar] [CrossRef] [Green Version]
  19. Switanek, M.; Crailsheim, K.; Truhetz, H.; Brodschneider, R. Modelling seasonal effects of temperature and precipitation on honey bee winter mortality in a temperate climate. Sci. Total Environ. 2017, 579, 1581–1587. [Google Scholar] [CrossRef]
  20. Ratti, V.; Kevan, P.G.; Eberl, H.J. A mathematical model of forager loss in honeybee colonies infested with Varroa destructor and the acute bee paralysis virus. Bull. Math. Biol. 2017, 79, 1218–1253. [Google Scholar] [CrossRef]
  21. Becher, M.A.; Osborne, J.L.; Thorbek, P.; Kennedy, P.J.; Grimm, V. Review: Towards a systems approach for understanding honeybee decline: A stocktaking and synthesis of existing models. J. Appl. Ecol. 2013, 50, 868–880. [Google Scholar] [CrossRef] [Green Version]
  22. Torres, D.J.; Ricoy, V.M.; Roybal, S. Modelling honey bee populations. PLoS ONE 2015, 10, e0130966. [Google Scholar] [CrossRef]
  23. Yıldız, T.A. A fractional dynamical model for honeybee colony population. Int. J. Biomath. 2018, 11, 1850063. [Google Scholar] [CrossRef]
  24. Romero-Leiton, J.P.; Gutierrez, A.; Benavides, I.F.; Molina, O.E.; Pulgarín, A. An approach to the modeling of honey bee colonies. Web Ecol. 2013, 22, 7–19. [Google Scholar] [CrossRef]
  25. Atanasov, A.Z.; Georgiev, S.G.; Vulkov, L.G. Reconstruction analysis of honeybee colony collapse disorder modeling. Optim. Eng. 2021, 22, 2481–2503. [Google Scholar] [CrossRef]
  26. Atanasov, A.Z.; Georgiev, S.G.; Vulkov, L.G. Parameter identification of Colony Collapse Disorder in honeybees as a contagion. In Modelling and Development of Intelligent Systems; Simian, D., Stoica, L.F., Eds.; Springer: Dordrecht, The Netherlands, 2021; Volume 1341, pp. 363–377. [Google Scholar]
  27. Hong, W.; Chen, B.; Lu, Y.; Lu, C.; Liu, S. Using system equalization principle to study the effects of multiple factors to the development of bee colony. Ecol. Model. 2022, 470, 110002. [Google Scholar] [CrossRef]
  28. Hundsdorfer, W.; Vermer, J. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations; Springer: Berlin/Heidelberg, Germany, 2003. [Google Scholar]
  29. Marchuk, G.I. Adjoint Equations and Analysis of Complex Systems; Kluwer: Dordrecht, The Netherlands, 1995. [Google Scholar]
  30. Marchuk, G.I.; Agoshkov, V.I.; Shutyaev, V.P. Adjoint Equations and Perturbation Algorithms in Nonlinear Problems; CRC Press: Boca Raton, FL, USA, 1996. [Google Scholar]
  31. Winston, W.L. The Biology of the Honey Bee; Harvard University Press: Cambridge, MA, USA, 1991. [Google Scholar]
  32. Ma, C.; Jiang, L. Some research on Levenberg–Marquardt method for the nonlinear equations. Appl. Math. Comp. 2007, 184, 1032–1040. [Google Scholar] [CrossRef]
Figure 1. Schematic representation of model (1)–(4).
Figure 1. Schematic representation of model (1)–(4).
Axioms 12 00214 g001
Figure 2. Solution to the direct problem (1)–(5).
Figure 2. Solution to the direct problem (1)–(5).
Axioms 12 00214 g002
Figure 3. Phase space diagram: non-trivial equilibrium.
Figure 3. Phase space diagram: non-trivial equilibrium.
Axioms 12 00214 g003
Figure 4. Phase space diagram: extinction equilibrium.
Figure 4. Phase space diagram: extinction equilibrium.
Axioms 12 00214 g004
Table 1. Simulation with ε p = 8 × 10 4 .
Table 1. Simulation with ε p = 8 × 10 4 .
Parameter p i p 0 i p ˇ i | p i p ˇ i | | p i p ˇ i | p i r i
α 0.0180.020.02740.00940.5234 4 × 10 23
γ 0.120.100.09910.02090.1746 7 × 10 5
δ 0.5710.500.36330.20770.3638 5 × 10 23
σ 0.10.200.18860.08860.8862 1 × 10 25
ω 0.951.000.89920.05080.0535 1 × 10 23
Table 2. Simulation with perturbed observations and ε p = 8 × 10 4 .
Table 2. Simulation with perturbed observations and ε p = 8 × 10 4 .
Parameter p i p 0 i p ˇ i | p i p ˇ i | | p i p ˇ i | p i r i
α 0.0180.020.02690.00890.4942 4 × 10 23
γ 0.120.100.09890.02110.1758 7 × 10 5
δ 0.5710.500.37310.19790.3465 5 × 10 23
σ 0.10.200.18690.08690.8689 1 × 10 25
ω 0.951.000.89270.05730.0603 1 × 10 23
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Atanasov, A.Z.; Georgiev, S.G.; Vulkov, L.G. Parameter Estimation Analysis in a Model of Honey Production. Axioms 2023, 12, 214. https://doi.org/10.3390/axioms12020214

AMA Style

Atanasov AZ, Georgiev SG, Vulkov LG. Parameter Estimation Analysis in a Model of Honey Production. Axioms. 2023; 12(2):214. https://doi.org/10.3390/axioms12020214

Chicago/Turabian Style

Atanasov, Atanas Z., Slavi G. Georgiev, and Lubin G. Vulkov. 2023. "Parameter Estimation Analysis in a Model of Honey Production" Axioms 12, no. 2: 214. https://doi.org/10.3390/axioms12020214

APA Style

Atanasov, A. Z., Georgiev, S. G., & Vulkov, L. G. (2023). Parameter Estimation Analysis in a Model of Honey Production. Axioms, 12(2), 214. https://doi.org/10.3390/axioms12020214

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop