Analysis of the Influence of Brood Deaths on Honeybee Population

Abstract

Featured Application

This model can be utilized to predict the long-term impacts of insecticide exposure on honeybee populations, aiding in the development of effective conservation strategies and regulatory policies to protect the honeybees as crucial pollinators.

Abstract

Many mathematical models using ordinary differential equations (ODEs) have been used to investigate what type of stressors cause honeybee colonies collapse. We propose a simple model of a delayed differential equation system (DDE) to describe the effect of insecticides over brood death rate and its influence over honeybee population dynamics. First, we remember some basic facts for the model with no delay. To analyze our model, we study the equilibria and perform stability and sensitivity analysis of the DDE system. Next, by using the delay time $\tau$ as a bifurcation parameter, we find that no Hopf bifurcation could arise as the time lag $\tau$ varies within biologically plausible ranges. Numerical simulations with real data are studied for the biological significance of the model.

1. Introduction

Honeybee colonies consist of distinct castes, including broods (encompassing eggs and larvae), sealed broods (pupae), queens, male drones, and female worker bees. Drones primarily exist to fertilize a queen during seasonal mating flights, after which they promptly perish. The queen’s main role is egg-laying. When a virgin queen hatches, she actively seeks out and eliminates rival virgin queens before finding mates to begin her reproductive duties.

The transformation from egg to adult bee takes roughly 21 days. Upon reaching adulthood, worker bees gradually transition from performing tasks within the hive to those outside as they age. These roles range from brood care and hive upkeep to foraging. Honeybees produce food not only for themselves but also for the queen and the developing brood (eggs and larvae). Beyond honey production for sustenance, they play a crucial role in pollinating a wide array of plant species.

Honeybee colonies exhibit a highly organized social structure, where each caste performs specific functions essential for the colony’s survival and growth. Worker bees, the most numerous caste, ensure the smooth functioning of the hive by performing duties that change as they age. Initially, young workers, also known as nurse bees, tend to the brood by feeding and grooming them. As they grow older, they transition to tasks such as cleaning the hive, processing food, and eventually foraging for nectar and pollen.

Queens are central to the colony’s reproduction and genetic diversity. A newly emerged virgin queen embarks on a series of orientation flights, followed by mating flights where she mates with multiple drones. This genetic diversity is crucial for the colony’s resilience to diseases and environmental changes. After mating, the queen returns to the hive to lay eggs, ensuring the colony’s continuity.

Drones, while having a singular focus on reproduction, contribute significantly to the genetic health of the bee population. Their presence during the mating season ensures that queens mate with drones from different colonies, thus promoting genetic variation.

The colony’s dynamic lifecycle and division of labor are critical for its success. Worker bees’ ability to adapt to different roles based on the colony’s needs demonstrates their remarkable flexibility and the complex social structure of the hive. This adaptability is key to the colony’s ability to thrive in varying environmental conditions.

Understanding the sophisticated dynamics of honeybee colonies, including the roles and lifecycle of different castes, is fundamental for developing accurate models that predict their behavior and response to stressors. These models are essential for devising strategies to protect and sustain honeybee populations, which are vital for biodiversity and agriculture through their pollination activities.

Honeybee populations typically follow a seasonal cycle influenced by the availability of resources [1]. In spring (April through June), the lengthening days and increased pollen supply lead to a birth rate that exceeds the death rate, causing the population to expand. The population peaks toward the end of June and stays elevated until mid-summer, after which it starts to decline. This downward trend continues into autumn and winter, reaching its lowest point in February.

Regarding brood dynamics, the queen’s egg-laying activity is minimal during late winter but rises sharply in spring to take advantage of the flowering season. In summer, the egg-laying rate remains elevated but slightly lower than spring levels. In autumn, certain plant species like ironweed bloom, triggering a second, smaller peak in egg-laying. After this brief increase, the queen’s egg-laying rate drops to zero by December [2,3].

The seasonal fluctuations in honeybee populations are primarily driven by environmental factors and the availability of resources. In the spring, longer days and abundant pollen from blooming flowers provide optimal conditions for the queen to increase her egg-laying rate. This results in a surge in the number of brood and, subsequently, adult bees. The peak population during late June supports the hive’s activities, including foraging and storing food for the coming months.

As summer progresses, the initial abundance of resources begins to wane. The reduction in pollen and nectar sources causes the population growth to slow and eventually decline. By mid-summer, the hive experiences a gradual decrease in population as the death rate begins to exceed the birth rate. This trend continues into the fall, with the population steadily declining as fewer resources are available and the queen’s egg-laying rate diminishes.

During the late fall and winter, the queen’s egg-laying rate significantly drops, reaching minimal levels. The colony’s focus shifts to conserving energy and resources, resulting in the lowest population levels in February. The winter months are a period of survival, with the colony relying on stored food and reducing activity to maintain warmth and protect the remaining bees.

The spring resurgence in egg-laying and population growth is a critical period for the hive. The increased egg-laying rate in response to the blooming period ensures that the colony can replenish its numbers and prepare for the next cycle of growth and decline. The potential second peak in the fall, triggered by the blooming of specific plant species, provides an additional boost to the colony’s numbers before the onset of winter.

Understanding these seasonal dynamics is essential for managing and protecting honeybee populations. Accurate models that incorporate these patterns can help predict population changes and identify critical periods for intervention. Such models are invaluable for beekeepers and researchers in developing strategies to sustain healthy colonies and mitigate the effects of environmental stressors.

Extensive research has been carried out to explore the factors contributing to the collapse of honeybee colonies [4,5]. Potential causes of colony decline include parasites, diseases, agricultural pesticides, chemicals used by beekeepers, genetically modified crops, insecticides, and shifts in cultural practices [6,7].

To tackle this issue and address the difficulty of directly monitoring honeybee population dynamics, the development of mathematical models becomes essential. These models help identify key factors that affect bee populations. Various models have been formulated to simulate honeybee population dynamics and to assess the influence of human-induced stressors, such as pathogens, parasites, and nutritional challenges, on the colony’s long-term survival [7,8,9,10,11,12,13,14,15,16].

Among these previous studies, the work of Khoury et al. (2011) stands out as highly relevant. These studies provide valuable insights into how various stressors and factors affect honeybee populations. Khoury et al.’s model highlighted the critical role that forager death rates play in colony dynamics, demonstrating how increased mortality among foragers can lead to significant declines in colony population and growth [13,14]. Russell et al. expanded on this by incorporating seasonal variations and resource availability, providing a more comprehensive view of the factors that influence colony health throughout the year [16].

Perry et al.’s case study focused on the practical implications of food shortages and rapid maturation, illustrating how these factors can directly impact colony success. Their findings emphasize the importance of adequate food supply and optimal maturation rates for maintaining healthy bee populations [7].

The development of such mathematical models, including ones with memory [17,18] is essential for predicting and managing the health of honeybee colonies. By understanding the dynamics of honeybee populations and the impact of various stressors, researchers and beekeepers can develop strategies to mitigate the effects of these stressors and promote the long-term viability of colonies. This research underscores the need for continued investigation into the complex interplay of factors affecting honeybee populations and the development of robust models that can inform effective management practices.

These models serve as crucial tools for beekeepers, allowing them to anticipate potential issues and take proactive measures to protect their colonies. With accurate predictions and a deeper understanding of colony dynamics, beekeepers can implement targeted interventions to address specific stressors, ultimately contributing to the sustainability and resilience of honeybee populations.

Despite their contributions, these models overlook the mortality rate of broods, concentrating instead on the death rates of foraging bees. As a result, they fail to directly investigate how brood mortality might contribute to colony collapse disorder (CCD) [19]. The study in [20] introduces a mathematical model that specifically evaluates the impact of brood death rates, resulting from insecticide exposure, on honeybee population dynamics. This model also accounts for the seasonal variations in the queen’s egg-laying rate.

The novelty of our study lies in explicitly modeling brood mortality as a function of insecticide exposure and analyzing its delayed effects on the colony’s stability. This provides new insights into how external stressors propagate through the colony over time. By demonstrating the biological significance of these delays, our work offers a an approach to the development of more sophisticated predictive tools in apiculture. This advancement situates the proposed model as a step forward in understanding and mitigating risks to honeybee populations, which are crucial for global agriculture and ecosystem sustainability.

The theoretical part of the paper [20] studies the stability of the brood honeybee model described and investigates how the brood delay $\tau$ affects population dynamics within the hive.

We first study the local stability of the equilibrium points of the ODEs. Next, using the delay $\tau$ as a bifurcation parameter, we demonstrate that the positive equilibrium loses its stability and the equation cannot exhibit a Hopf bifurcation.

Our goal is to use the theoretical results from this model to obtain realistic predictions for honeybee population levels and to identify critical periods of growth for honeybee colonies in villages such as Brestovitsa, Yuper, etc. [21].

By examining the brood death rate, this model addresses a crucial aspect that has been overlooked in previous studies. The incorporation of brood mortality into the model allows for a more comprehensive understanding of the factors influencing honeybee population dynamics. This approach is particularly important for predicting potential colony collapse due to increased brood mortality. The stability analysis and investigation of brood delay $\tau$ provide insights into how time delays in brood development affect overall colony stability [22].

These theoretical insights are not merely academic; they have practical applications for beekeepers and researchers. By identifying critical periods in the growth and stability of honeybee colonies, this model can help beekeepers implement timely interventions to mitigate adverse effects and support colony health.

Applying these theoretical results to real-world scenarios, such as the honeybee colonies in villages, allows for more accurate predictions and effective management strategies. Understanding when colonies are most vulnerable enables beekeepers to take proactive measures, such as adjusting beekeeping practices or applying treatments to protect the brood.

The remaining part of the paper is organized as follows. In the next section, we introduce the ODE and DDE models. Some results, obtained in [20], and new ones for the ODE model are discussed in Section 3. Results concerning the existence, uniqueness, as well as the stability of the solutions for the DDE model are reported in Section 4. Computational experiments with synthetic and real input data are performed in Section 5. The last section provides conclusions for the modelers and beekeepers on how to select the specific parameters and select the best model for honeybee populations in the hive.

2. Model Formulation

Among previous studies, such as the review by [23], one of the most relevant is the work by [14], which focuses primarily on the modeling of food availability and its impact on honeybee populations.

In their study, Khoury et al. extended the basic model of [13] to explicitly include both food and brood. It is assumed that food availability influences both the number of brood successfully reared to adulthood and the rate at which bees transition from hive duties to foraging. In this section, we introduce two models. The first one is the ODE model proposed in [20]. Following the concept of [13,14], we then extend the ODE model of [20] by incorporating a time delay in the number of brood.

Here, we focus on a simpler model to analyze in detail the effect of the delay in food availability on hive population dynamics.

2.1. Biological Assumptions

Let $B\left(t\right)$ represent the number of broods (eggs, larvae, and pupae) in a hive at time t, and $H\left(t\right)$ denote the number of adult bees at the same time. The time t is measured in 12-day intervals, reflecting the data collection period. The following ecological assumptions describe the population dynamics of western honeybees (Apis mellifera):

• The queen lays eggs at a rate r, adding to the brood compartment. It is assumed that the queen is consistently fertile and well nourished.

• The brood receives care from adult bees, following a Holling Type 3 functional response. This model assumes that a threshold number of adults, K, is required for adequate brood care. Without sufficient adults, the colony faces collapse, while a sufficient number ensures its prosperity.

• Brood mortality follows an exponential decay pattern, proportional to the brood population at time t. Similarly, adult bee mortality is also modeled as an exponential decay, proportional to the adult population at time t.

• Brood develop into adult bees at a constant rate $\varphi$ is termed the eclosion rate.

• The model excludes drones, assuming that their primary role is mating with queens, which is not considered in this system.

We consider a Holling Type 3 functional response for modeling the “birth” of brood, as proposed by Holling in 1975. This response accounts for the capture rate, which increases in direct proportion to the adult population. At a low adult population, the rate of brood capture accelerates. Let $t_s$ denote the available capture time, which is the period during which adults can search for and care for the brood.

These assumptions lead to the following system of nonlinear ordinary differential equations (ODEs) modeling the honeybee population dynamics.

2.2. The ODE Model

In this model, r represents the egg-laying rate, and K denotes the critical population of adult bees over time. K is considered the critical number of bees required to care for the brood.

The available adult bees must locate and tend to more brood for the hive to succeed. Following the aforementioned assumptions, the authors of [20] proposed the following ordinary differential equation (ODE) model:

(1)$\begin{array}{c}{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}=\frac{r{H}^2}{K^2+H^2}-{\mu }_{b}B-\varphi B\equiv f(B,H),}\hfill \\ {\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}=\varphi B-{\mu }_{h}H\equiv g(B,H)}\hfill \end{array}$

A mathematical analysis of the model in Equation (1) is presented in [20, Section III].

It is shown that the solutions of the system in Equation (1) are bounded and positive, and Equation (1) has no limit cycles in ${\mathbb{R}}^2$.

Thus, the system in Equation (1) is biologically well defined, ensuring that the population does not grow unbounded over time. The maximum value of each compartment corresponds to the maximum population of broods and adults that can exist [24].

For the brood population, the upper limit is proportional to the egg-laying rate r. A lower value of r results in a smaller bound. Similarly, if the combined rate of brood mortality and maturation, ${\mu }_b+\varphi$, is excessively high, the upper limit will also be reduced.

Therefore, for practical applications that assist beekeepers, it is necessary to estimate the parameters of the system in Equation (1) based on concrete measurements.

The model in Equation (1) follows the studies performed in papers [13,14].

By utilizing these insights and mathematical formulations, we aim to provide a robust framework for understanding and predicting honeybee population dynamics. This model serves as a valuable tool for researchers and beekeepers, enabling them to anticipate changes in colony health and take appropriate measures to support and sustain healthy bee populations. Through continuous refinement and validation against real-world data, the model can be further enhanced to address additional factors affecting bee populations, ultimately contributing to the preservation of these essential pollinators.

2.3. The DDE Model

Following the basic models from [14,25], we extend Equation (1) to

(2)${\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}={\displaystyle \frac{r{H}^2}{K^2+H^2}}-{\mu }_{b}B-\varphi B,$

(3)${\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}=\varphi B(t-\tau )-{\mu }_{h}H,$

where the term $\varphi B(t-\tau )$ represents the rate at which adult bees emerge from the brood population. The bees that emerge at time t are those that entered the population at time $t-\tau$.

The parameters $\{r,{\mu }_b,{\mu }_h,\varphi ,K\}\in S_{\mathrm{adm}}$, where the admissible set comes from biological significance and mathematical reasoning, and $S_{\mathrm{adm}}\equiv (0,50000)\times {(0,2)}^3\times (0,50000)$ [14,20].

The models in [14,25] represent the interactions between hive bees, foragers, food, and brood. The dependent variables include the stored food in the hive f, the uncapped brood items B, the number of hive bees H, and the number of foragers F. Hive bees emerge $\tau$ days after pupation. Bees transition from the hive bee class to become foragers, a process modeled by the recruitment function $R(H,F,f)$, where f is a measure of the amount of food stored in the hive and available for colony use [14].

As noted in [25], these models are complex and cannot be easily analyzed mathematically. Therefore, in Equation (3) above, the function $R(H,F,f)$ is approximated by the constant ${\mu }_h$.

This delay differential equation model (DDE) captures the essential dynamics of honeybee populations by incorporating a time delay $\tau$ in the transition from brood to adult bees. The inclusion of this time delay is critical, as it reflects the biological reality that there is a lag between when brood are initially cared for and when they emerge as adult bees.

The analysis of such a delayed model provides insights into the stability and behavior of honeybee populations under different conditions. By examining how the delay $\tau$ affects population dynamics, we can better understand the potential for periodic oscillations or other complex behaviors that might arise in the colony.

Through this approach, we aim to extend the applicability of previous models by incorporating realistic biological delays, thus providing a more accurate and comprehensive framework for studying honeybee population dynamics. This model not only aids in theoretical research but also offers practical tools for beekeepers to manage their colonies more effectively by understanding the critical periods and factors that influence bee emergence and overall colony health.

2.4. Sensitivity Analysis

In this subsection, we will conduct a sensitivity analysis of the basic coexistence number ${\mathcal{R}}_d$ with respect to its parameters [20]:

${\displaystyle {\mathcal{R}}_d=\frac{r\varphi }{2K{\mu }_h({\mu }_b+\varphi )}.}$

Essentially, the sensitivity analysis measures the impact of each one of the parameters in the model. The sensitivity of the basic coexistence number with respect to a particular parameter is the amount of relative change of ${\mathcal{R}}_d$ to relative change in p. If ${\mathcal{R}}_d$ is influenced by and differentiable with respect to all parameters, then the sensitivity indices could be denoted by means of partial derivatives in the following way:

${\displaystyle {\mathcal{S}}_p^{{\mathcal{R}}_d}=\frac{p}{{\mathcal{R}}_d}·\frac{\partial {\mathcal{R}}_d}{\partial p}.}$

In the case of the models in Equations (1)–(3),

$\begin{array}{cc}{\displaystyle {\mathcal{S}}_r^{{\mathcal{R}}_d}=1>0,}\hfill & {\displaystyle {\mathcal{S}}_{{\mu }_h}^{{\mathcal{R}}_d}=-1<0,}\hfill \\ {\displaystyle {\mathcal{S}}_{\varphi }^{{\mathcal{R}}_d}=\frac{{\mu }_b}{{\mu }_b+\varphi }>0,}\hfill & {\displaystyle {\mathcal{S}}_K^{{\mathcal{R}}_d}=-1<0,}\hfill \\ \hfill & {\displaystyle {\mathcal{S}}_{{\mu }_b}^{{\mathcal{R}}_d}=-\frac{{\mu }_b}{{\mu }_b+\varphi }<0.}\hfill \end{array}$

3. Well-Posedness and Stability Analysis for the ODE Model

The starting points of the stability analysis of the system in Equations (1)–(3) is the study of steady states, i. e., equilibria.

The point $E_0(0,0)$ is the trivial one and solving the system $f(u,v)=0,g(u,v)=0$ of algebraic equations provides the so-called endemic equilibria for the model Equation (1).

They were found in [20] and the results are as follows:

$E_2(e_{21}, e_{22}),E_3(e_{31}, e_{32}),$

$e_{21}={\displaystyle \frac{r\varphi +\xi }{2\varphi ({\mu }_b+\varphi )}},e_{22}={\displaystyle \frac{r\varphi +\xi }{2{\mu }_h({\mu }_b+\varphi )}},$

$e_{31}={\displaystyle \frac{r\varphi -\xi }{2\varphi ({\mu }_b+\varphi )}},e_{32}={\displaystyle \frac{r\varphi -\xi }{2{\mu }_h({\mu }_b+\varphi )}},$

where

$\xi =\sqrt{r^2{\varphi }^2-4K^2{\mu }_h^2{({\mu }_b+\varphi )}^2}.$

The Jacobian to the system Equation (1) is as follows: in the trivial equilibrium,

$J(0,0)=\left|\begin{array}{cc}-({\mu }_b+\varphi )& \varphi \\ \varphi & -{\mu }_h\end{array}\right|,det\left(J\right)={\mu }_h({\mu }_b+\varphi )>0,\mathrm{tr}\left(J\right)=-({\mu }_b+{\mu }_h+\varphi );$

the second equilibrium,

$J(B^{*},H^{*})=\left|\begin{array}{cc}-({\mu }_b+\varphi )& {\displaystyle \frac{{\mu }_h({\mu }_b+\varphi )(r\varphi -\xi )}{r{\varphi }^2}}\\ \varphi & -{\mu }_h\end{array}\right|,$

$det\left(J\right)={\displaystyle \frac{{\mu }_h({\mu }_b+\varphi )\xi }{r\varphi }},\mathrm{tr}\left(J\right)=-({\mu }_b+{\mu }_h+\varphi );$

and the third equilibrium,

$J(B^{*},H^{*})=\left|\begin{array}{cc}-({\mu }_b+\varphi )& {\displaystyle \frac{{\mu }_h({\mu }_b+\varphi )(r\varphi +\xi )}{r{\varphi }^2}}\\ \varphi & -{\mu }_h\end{array}\right|,$

$det\left(J\right)=-{\displaystyle \frac{{\mu }_h({\mu }_b+\varphi )\xi }{r\varphi }},\mathrm{tr}\left(J\right)=-({\mu }_b+{\mu }_h+\varphi ).$

Therefore, the trace will always be negative, independent of $B^{*},H^{*}$.

So, if $det\left(J(B^{*},H^{*})\right)>0$, the equilibrium point $(B^{*},H^{*})$ will be asymptotically stable and, if $det\left(J(B^{*},H^{*})\right)<0$, the equilibrium point $(B^{*},H^{*})$ will be unstable (for non-zero positive parameters).

Due to $J(0,0)$, the trivial equilibrium ($B^{*} = 0,H^{*} = 0$) is always stable.

4. Stability Analysis for the DDE Model

Let

$u_1\left(t\right)=B\left(t\right)-B^{*},u_2\left(t\right)=H\left(t\right)-H^{*},$

where the point $(B^{*}, H^{*})$ is one of the equilibria points above.

The linearization of Equations (2) and (3) at $E_1(0,0)$ is

(4)${\displaystyle \frac{\mathrm{d}{u}_1}{\mathrm{d}t}}=-{\mu }_b{u}_1\left(t\right)-\varphi u_1\left(t\right)+{\displaystyle \frac{2r{K}^2{H}^{*}}{{({\left(H^{*}\right)}^2+K^2)}^2}}{u}_2\left(t\right),$

(5)${\displaystyle \frac{\mathrm{d}{u}_2}{\mathrm{d}t}}=\varphi u_1(t-\tau )-{\mu }_h{u}_2\left(t\right).$

Its characteristic equation is

(6)${\lambda }^2+({\mu }_h+q)\lambda +{\mu }_{h}q-s\varphi e^{-\lambda \tau }=0,$

where

$q={\mu }_b+\varphi ,s={\displaystyle \frac{2r{K}^2{H}^{*}}{{\left({\left(H^{*}\right)}^2+K^2\right)}^2}}.$

When there is no delay, that is, $\tau =0$, the corresponding Equation (6) reduces to

(7)${\lambda }^2+({\mu }_h+q)\lambda +{\mu }_{h}q-s\varphi =0,$

and the corresponding eigenvalues were reported in the previous section.

By the Hartman–Grobman theorem [26], since the eigenvalues of the linearized system in Equations (4) and (5) have non-zero real parts, the quantitative behavior of solutions for the nonlinear system in Equations (2) and (3) is the same as the linearized system in Equations (4) and (5) in a neighborhood of the equilibrium point $E(B^{*},H^{*})$.

Now, we investigate the distribution of roots of the transcendental Equation (6) since the stability of the point $(0,0)$ of the linear system in Equations (4) and (5) is governed by the locations of the roots of the characteristic Equation (6).

By the roots of Equation (6), and the results of Section 3, there exists ${\tau }_0$ such that $\mathfrak{Re}\left(\lambda \right(\tau \left)\right) < 0$ for $\tau \in [0,{\tau }_0)$.

Since a loss of asymptotic stability arises when $\mathfrak{Re}\left(\lambda \left({\tau }^{*}\right)\right) = 0$, we examine whether there exists a ${\tau }^{*} > 0$ for which $\mathfrak{Re}\left(\lambda \left({\tau }^{*}\right)\right) = 0$, that is, we would like to know when Equation (6) has purely imaginary roots.

We suppose that, for $\tau = {\tau }^{*}$, we have $\lambda = iw$ with $w > 0$; then, we have the following result.

Lemma 1.

For the system in Equations (4) and (5), the transcendental Equation (6) has one purely imaginary root.

Proof. 

For $\tau = {\tau }^{*}$, let $\lambda = iw$ be the root of Equation (6) with w real and without loss of generality $w > 0$.

Then,

$i^2{w}^2+({\mu }_h+q)iw+{\mu }_{h}q-s\varphi e^{-iw\tau }=0,$

that is,

$-w^2+({\mu }_h+q)iw+{\mu }_{h}q-s\varphi (cosw\tau -isinw\tau )=0.$

Separating real and imaginary parts,

(8)$w^2-{\mu }_{h}q=-s\varphi ·cosw\tau$

and

(9)$-({\mu }_h+q)w=s\varphi ·sinw\tau ,$

which is equivalent to

$w^4+({\mu }_h^2+q^2)w^2-s^2{\varphi }^2+{\left({\mu }_{h}q\right)}^2=0.$

Letting $w^2 = t$, we have

(10)$t^4+({\mu }_h^2+q^2)t-s^2{\varphi }^2+{\left({\mu }_{h}q\right)}^2=0.$

which implies that this equation governs the possible values of $\tau$ and w, for which Equation (6) can have purely imaginary roots.

Solving Equation (10) with respect to $w = \sqrt{t}$, we find

(11)$w=\sqrt{{\displaystyle \frac{-({\mu }_h^2+q^2)+\sqrt{{({\mu }_h^2+q^2)}^2+4s^2{\varphi }^2-4{\left({\mu }_{h}q\right)}^2}}{2}}}$

and, from Equation (8), we obtain

(12)${\tau }_k={\displaystyle \frac{2}{w}}arctan\left(\sqrt{-{\displaystyle \frac{w^2-{\mu }_{h}q+s}{w^2-{\mu }_{h}q-s}}}\right)+{\displaystyle \frac{2\pi k}{w}},k=0,1,2,\dots ,$

while from Equation (9), we derive

(13)${\tau }_k={\displaystyle \frac{2}{w}}arctan\left({\displaystyle \frac{s-\sqrt{s^2-{({\mu }_h+q)}^2{w}^2}}{({\mu }_h+q)w}}\right)+{\displaystyle \frac{2\pi k}{w}},k=0,1,2,\dots ,$

where formulae (12) and (13) are equivalent. This finalizes the proof of the lemma. □

Theorem 1.

For the DDEs modeling the population dynamics of Apis Mellifera, Equations (2) and (3), the following assertions (properties) are fulfilled:

• (a) 

  If $\tau \in [0,{\tau }_0)$, the equilibrium point $(0,0)$ of the system in Equations (2) and (3) is asymptotically stable.

• (b) 

  If $\tau > {\tau }_0$, then the equilibrium point $(B^{*}, H^{*})$ of the system in Equations (2) and (3) is unstable.

• (c) 

  The system in Equations (2) and (3) does not undergo a Hopf bifurcation at the equilibrium point $(B^{*}, H^{*})$ with respect to the time lag τ.

Proof. 

The assertions (a) and (b) follow from Lemma 1.

In order to prove point (c), we differentiate the characteristic Equation (6) with respect to $\tau$:

$2\lambda {\displaystyle \frac{d\lambda }{d\tau }}+({\mu }_h+q){\displaystyle \frac{d\lambda }{d\tau }}+s\varphi e^{-\lambda \tau }\left(\tau {\displaystyle \frac{d\lambda }{d\tau }}+\lambda \right)=0.$

Hence,

${\left({\displaystyle \frac{d\lambda }{d\tau }}\right)}^{-1}=-{\displaystyle \frac{\tau }{\lambda }}-{\displaystyle \frac{2\lambda +{\mu }_h+q}{\lambda s\varphi e^{-\lambda \tau }}}.$

Therefore,

$\mathfrak{Re}{\left.\left({\left({\displaystyle \frac{d\lambda }{d\tau }}\right)}^{-1}\right)\right|}_{\lambda =iw}=\mathfrak{Re}\left[{\displaystyle \frac{{\tau }_{k}i}{w}}\right]+\mathfrak{Re}\left[-{\displaystyle \frac{2iw+q+{\mu }_h}{s\varphi iw(cosw{\tau }_k-isinw{\tau }_k)}}\right]$

$={\displaystyle \frac{{\mu }_h^2+q^2+2w^2}{s^2{\varphi }^2}}>0.$

Applying Rouche’s theorem, see, e.g., [27], we find that the transversality condition is satisfied, so the only condition for exhibiting a Hopf bifurcation is that w is real, or

$s\varphi >{\mu }_{h}q$

from Equation (11). This condition is never satisfied for values of parameters belonging to $S_{\mathrm{adm}}$. □

5. Computational Simulations

In this section, we present implicit numerical methods for ODEs and DDEs. To begin, we give the algorithms, then we conduct computational experiments with real data [28]. First, we show the results with the classical model in Equation (1). Further, we display the outcome with the delayed model in Equations (2) and (3).

5.1. Implicit Numerical Methods

The choice of a numerical method is very important for the accuracy and computational cost of the solution. Usually, for complex nonlinear systems, the implicit methods perform better in terms of computational resources if a predefined accuracy is required.

The problem with the classical implicit methods is that they are very computationally demanding, since, at each step, for each stage, we need to solve a system of nonlinear algebraic equations. It is done by a suitable iteration method, which makes the entire approach computationally intensive. A possible way to try to diminish the load is to apply a special choice of the coefficients in Equation (14). We employ the Mono-Implicit Runge–Kutta (MIRK) method [29].

If the ODE to solve has the form

${\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}=f\left(t,y\left(t\right)\right),y\left(0\right)=y_0,$

then the MIRK method could be defined as follows:

(14)$y_{n+1}=y_n+h_n\sum _{r=1}^s{b}_r{K}_{n+1}^r,$

where s is the number of stages, $h_n=t_{n+1}-t_n$ is the size of the step, and

$K_{n+1}^r=f(t_n+c_r{h}_n,(1-v_r)y_n+v_r{y}_{n+1}+h_n\sum _{j=1}^{r-1}{a}_{rj}{K}_{n+1}^j).$

The other coefficients could be defined by means of the special form of the Butcher tableau in Equation (15):

(15)$\begin{array}{ccccccc}{c}_1& v_1& 0& 0& \cdots & 0& 0\\ c_2& v_2& a_{2,1}& 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ c_s& v_s& a_{s,1}& a_{s,2}& \cdots & a_{s,s-1}& 0\\ & & b_1& b_2& \cdots & \cdots & b_s\end{array}\equiv \begin{array}{ccc}C& V& A\\ & & B\end{array}.$

For consistency, it must be true that $C = V + A\ast e$, where the elements of vector $C = {(c_1,c_2,\dots ,c_s)}^{\top }$ are called abscissae, $B = (b_1,b_2,\dots ,b_s)$, ${\displaystyle \sum _{i=1}^s{b}_i = 1}$, $V = {(v_1,v_2,\dots ,v_s)}^{\top }$, $e = {(1,1,\dots ,1)}^{\top }$, and $A = {\left\{a_{ij}\right\}}_{i,j=1}^s$ is the Runge–Kutta matrix. The abscissae fulfill ${\displaystyle c_r=v_r + \sum _{i=1}^s{a}_{ri}}$.

The benefit of the MIRK over the other implicit methods is that the stages explicitly relate to each other, and they appear implicit only with respect to $y_{n+1}$. In this manner, the method retains the high accuracy of the implicit class, while the computational resources needed are significantly reduced.

In this case, our DDE looks like

${\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}} = f\left(t,y\left(t\right),y(t-\tau )\right),y\left(t\right) = y_0\left(t\right) for -\tau \le t \le 0,$

then it is possible to adapt the MIRK method for delayed equations [30].

Similar to Equation (14), the delayed MIRK is defined by

$y_{n+1}=y_n + h_n\sum _{r=1}^s{b}_r{K}_{n+1}^r,$

but now the stages are defined as

$K_{n+1}^r = f(t_n+c_r{h}_n, (1-v_r)y_n + v_r{y}_{n+1} + h_n\sum _{j=1}^{r-1}{a}_{rj}{K}_{n+1}^j, y(t_n+c_r{h}_n-\tau )).$

This modification is enough if the delay $\tau$ is multiple of the step $h_n$, since, then, the delayed term $y(t_n + c_r{h}_n - \tau )$ would always coincide with a grid point.

It is possible to obtain different schemes, if the coefficients A, V, and B vary. Of course, they need to meet certain conditions [31]. In our experiments, we used the third-order scheme [30].

$\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 1& 1& 0& 0& 0& 0\\ {\displaystyle \frac{3}{4}}& {\displaystyle \frac{27}{32}}& {\displaystyle \frac{3}{64}}& -{\displaystyle \frac{9}{64}}& 0& 0\\ {\displaystyle \frac{1}{5}}& {\displaystyle \frac{13}{125}}& {\displaystyle \frac{16}{125}}& -{\displaystyle \frac{4}{125}}& 0& 0\\ & & {\displaystyle \frac{5}{9}}& -{\displaystyle \frac{3}{8}}& {\displaystyle \frac{128}{99}}& -{\displaystyle \frac{125}{264}}\end{array}.$

5.2. Data Description

The data we used were collected and described by Harris [28]. They concern the foraging seasons during two years—1975 and 1976. The colonies were initiated with approximately 7000 adult bees, and were used for commercial purposes. The data were recorded on a 12-day basis. To estimate the bee numbers, some colonies were killed at various times throughout the season and the corpses were counted. Then, time series for eggs, larvae, pupae, and adult bees <19 days and >19 days of age were provided for both years. In addition, the collected honey was also measured.

Although old, this piece of data is often used in the literature due to its reliability. Furthermore, since many studies have employed it, different models could be compared, based on the same data.

The values of the parameters and initial conditions were suggested by [20], while the optimal value of the lag $\tau$ is proposed by [14].

5.3. Numerical Simulations for ODEs (Equation (1))

In this subsection, we test the model in Equation (1) with real data from [28] and set as follows. For the 1975 data, we used r = 25,061.572, ${\mu }_B = 0.575$, ${\mu }_H = 0.302$, $\varphi = 0.571$, and $K = 3805.467$ and the initial conditions were $B_0 = 6125$ and $H_0 = 5362$ [20].

For the 1976 data, we employed r = 27,666.177, ${\mu }_B = 0.583$, ${\mu }_H = 0.246$, $\varphi = 0.571$, and $K = 4732.389$ and the initial conditions were $B_0 = 5982$ and $H_0 = 5362$ [20]. The results are plotted in Figure 1.

It is observed that the hive is attracted to its non-trivial equilibrium $(B^{*}, H^{*})$ relatively fast.

The phase plane portraits are given in Figure 2.

5.4. Numerical Simulations for DDE (Equations (2) and (3))

In this subsection, we continue with the aforementioned parameters and initial values, but also set $\tau = 12$ [14]. The results are plotted in Figure 3.

In this case, the hive also achieves the equilibrium point $(B^{*},H^{*})$, but slower than the modelling without delay. This outcome bears much more resemblance to the reality.

The phase plane portraits are given in Figure 4.

5.5. Discussion

The inclusion of time delays allows for the identification of critical thresholds and the complex behavior of the quantities under study. As our results demonstrate, using the delay time $\tau$ enables modeling that non-delayed models fail to express adequately. Understanding these phenomena provides crucial insights into the stability of honeybee populations and their resilience to environmental stressors.

Incorporating time delays into our model offers a more realistic representation of honeybee colony dynamics by accounting for the inherent lag between key biological processes, such as egg-laying, brood development, and adult maturation. This temporal component aligns closely with the actual lifecycle of honeybees and their response to stressors like insecticides, which may not have immediate effects but influence population dynamics over time. In the context of precision apiculture, this enhanced realism enables beekeepers and researchers to predict colony responses to environmental changes or stressors with greater accuracy, thereby supporting proactive management strategies to prevent colony collapse. Compared to traditional models without delays, our approach captures critical thresholds and time-dependent interactions that could otherwise be overlooked, such as delayed population declines due to brood death. This finding situates our model within the broader context of mathematical models by addressing a key limitation in existing frameworks and demonstrating how time lag effects can refine predictions, improve intervention planning, and ultimately contribute to the sustainability of honeybee populations.

6. Conclusions

In this study, we presented a delayed ordinary differential equation system (DDE) to model the influence of brood deaths on honeybee population dynamics under the impact of insecticides. Our analysis has shown that incorporating time delays into the model provides significant advantages over traditional models without time lags.

One of the primary advantages of delayed modeling is its ability to more accurately capture the biological processes and temporal dynamics of honeybee populations. In real-world scenarios, the effects of environmental stressors such as insecticides do not manifest instantaneously. There is a natural delay between the exposure to the stressor, its impact on the brood, and the subsequent effect on the overall population. By including this time lag, our model offers a more realistic representation of these delayed effects, leading to more accurate predictions and insights.

Additionally, the delayed model enhances the sensitivity analysis, offering a deeper understanding of how different parameters influence the system’s behavior over time. This is particularly important for developing targeted intervention strategies. By simulating various scenarios with both synthetic and real data, we have illustrated the biological significance of the model, highlighting how delayed responses can lead to different population outcomes under identical conditions.

In the future, we plan to further study more advanced variants of the model with time-dependent egg-laying rate and more sophisticated variants of the recruitment function. Also, an important question that arises is the following: How does the coefficient parameters’ uncertainty affect the model predictions? It leads to the solution to an inverse coefficient DDE problem upon measurements of $B\left(t\right)$ and $H\left(t\right)$.

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.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding authors.

Acknowledgments

The authors are grateful to the anonymous referees for the useful suggestions and comments.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:

CCD	Colony Collapse Disorder
DDE	Delayed Differential Equation
MIRK	Mono-Implicit Runge–Kutta
ODE	Ordinary Differential Equation

Footnotes

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.

Figures

Figure 1. The solution to Equation (1) for 1975 data (left) and 1976 data (right).

Figure 2. The phase plane portraits of Equation (1) for 1975 data (left) and 1976 data (right).

Figure 3. The solution to Equations (2) and (3) for 1975 data (left) and 1976 data (right).

Figure 4. The phase plane portraits of Equations (2) and (3) for 1975 data (left) and 1976 data (right).

References

1. Chen, J.; Messan, K.; Messan, M.R.; DeGrandi-Hoffman, G.; Bai, D.; Kang, Y. How to model honeybee population dynamics: Stage structure and seasonality. Math. Appl. Sci. Eng.; 2020; 1, pp. 91-125. [DOI: https://dx.doi.org/10.5206/mase/10559]

2. Mid-Atlantic Apiculture Research and Extension Consortium (MAAREC). The Colony and Its Organization; 2017. Available online: https://canr.udel.edu/maarec/honey-bee-biology/the-colony-and-its-organization/ (accessed on 5 December 2024).

3. Mid-Atlantic Apiculture Research and Extension Consortium (MAAREC). Basic Bee Biology for Beekeepers; MAAREC, 2004; 1.4. Available online: https://bee-health.extension.org/basic-bee-biology-for-beekeepers/ (accessed on 5 December 2024).

4. Atanasov, A.Z.; Georgiev, S.G.; Vulkov, L.G. Reconstruction analysis of honeybee colony collapse disorder modeling. Optim. Eng.; 2021; 22, pp. 2481-2503. [DOI: https://dx.doi.org/10.1007/s11081-021-09678-0]

5. van Engelsdorp, D.; Evans, J.D.; Saegerman, C.; Mullin, C.; Haubruge, E.; Nguyen, B.K.; Frazier, M.; Frazier, J.; Cox-Foster, D.; Chen, Y. et al. Colony Collapse Disorder: A descriptive study. PLoS ONE; 2009; 4, e6481.

6. Oldroyd, B.P. What’s killing american honey bees?. PLoS Biol.; 2007; 5, e168. [DOI: https://dx.doi.org/10.1371/journal.pbio.0050168] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17564497]

7. Perry, C.J.; Søvik, E.; Myerscough, M.R.; Barron, A.B. Rapid behavioral maturation accelerates failure of stressed honey bee colonies. Proc. Natl. Acad. Sci. USA; 2015; 112, pp. 3427-3432. [DOI: https://dx.doi.org/10.1073/pnas.1422089112]

8. Aronstein, K.A.; Saldivar, E.; Vega, R.; Westmiller, S.; Douglas, A.E. How varroa parasitism affects the immunological and nutritional status of the honey bee, apis mellifera. Insects; 2012; 3, pp. 601-615. [DOI: https://dx.doi.org/10.3390/insects3030601]

9. Betti, M.I.; Wahl, L.M.; Zamir, M. Effects of infection on honey bee population dynamics: A model. PLoS ONE; 2014; 9, E110237. [DOI: https://dx.doi.org/10.1371/journal.pone.0110237]

10. Camazine, S.; Sneyd, J.; Jenkins, M.J.; Murray, J.D. A mathematical model of self-organized pattern formation on the combs of honeybee colonies. J. Theor. Biol.; 1990; 147, pp. 553-571. [DOI: https://dx.doi.org/10.1016/S0022-5193(05)80264-4]

11. Eberl, H.J.; Frederick, M.R.; Kevan, P.G. Importance of brood maintaince items in simple models of the honeybee varroa destructive-acute bee paralysis virus complex. Electron. J. Differ. Equ. Conf.; 2010; 19, pp. 85-98.

12. Kang, Y.; Rodriguez-Rodriguez, M.; Evilsizor, S. Ecological and evolutionary dynamics of two-stage models of social insects with egg cannibalism. J. Math. Anal. And Appl.; 2015; 430, pp. 324-353. [DOI: https://dx.doi.org/10.1016/j.jmaa.2015.04.079]

13. Khoury, D.S.; Myerscough, M.R.; Barron, A.B. A quantitative model of honey bee colony population dynamics. PLoS ONE; 2011; 6, e18491. [DOI: https://dx.doi.org/10.1371/journal.pone.0018491] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21533156]

14. Khoury, D.S.; Barron, A.B.; Meyerscough, M.R. Modelling food and population dynamics honey bee colonies. PLoS ONE; 2013; 8, e0059084. [DOI: https://dx.doi.org/10.1371/journal.pone.0059084] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23667418]

15. Kribs-Zaleta, C.M.; Mitchell, C. Modeling Colony Collapse Disorder in honeybees as a contagion. Math. Biosci. Eng.; 2014; 11, pp. 1275-1294. [DOI: https://dx.doi.org/10.3934/mbe.2014.11.1275] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25365602]

16. Russel, S.; Barron, A.B.; Harris, D. Dynamics modelling of honeybee (Apis mellifera) colony growth and failure. Ecol. Model.; 2013; 265, pp. 138-169.

17. Georgiev, S.G.; Vulkov, L.G. Parameter identification approach for a fractional dynamics model of honeybee population. Large-Scale Scientific Computing; Lirkov, I.; Margenov, S. Springer: Cham, Switzerland, 2022; Volume 13127, pp. 40-48.

18. Yıldız, T.A. A fractional dynamical model for honeybee colony population. Int. J. Biomath.; 1991; 11, 1850063.

19. 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, pp. 213-219. [DOI: https://dx.doi.org/10.1016/j.jtbi.2017.03.009]

20. Cascante, J.E.; Rojas, M.C.; Salinas, A.; Serna, M.; Torres, H.D.; Varga, E.; Cordovez, J.M.; Artega, J.R. Influence of brood deaths on Honey Bee population dynamics and the potential impact of insecticides. IBIO 4299 Int. Res.; 2017; Available online: https://www.academia.edu/download/55628215/ires-2017-honey-2.pdf (accessed on 5 December 2024).

21. Atanasov, A.Z.; Georgiev, S.G.; Vulkov, L.G. Parameter identification of Colony Collapse Disorder in honeybees as a contagion. Modelling and Development of Intelligent Systems; Simian, D.; Stoica, L.F. Springer: Dordrecht, The Netherland, 2021; Volume 1341, pp. 363-377.

22. Rihan, F.A. Delay Differential Equations and Applications to Biology. Forum for Interdisciplinary Mathematics; Springer: Singapore, 2021.

23. Chen, J.; DeGrandi-Hoffman, G.; Ratti, V.; Kang, Y. Review on mathematical modeling of honeybee population dynamics. Math. Biosci. Eng.; 2021; 18, pp. 9606-9650. [DOI: https://dx.doi.org/10.3934/mbe.2021471]

24. Winston, W.L. The Biology of the Honey Bee; Harvard University Press: Cambridge, MA, USA, 1991.

25. Myerscough, 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. The Role and Importance of Mathematics in Innovation: Proceedings of the Forum “Math-for-Industry” 2015; Springer: Singapore, 2017; pp. 35-50.

26. Hale, J. Theory of Functional Differential Equations; Springer: New York, NY, USA, 1977.

27. Kuang, Y. Delay Differential Equations with Applications. Population Dynamics; Academic Press: Boston, MA, USA, 1993.

28. Harris, J.L. A Population Model and Its Application to the Study of Honey Bee Colonies. Master’s Thesis; Department of Entomology, The University of Manitoba: Winnipeg, MB, Canada, 1980.

29. Muir, P.H. Optimal discrete and continuous mono-implicit Runge-Kutta schemes for BVODEs. Adv. Comput. Math.; 1999; 10, pp. 135-167. [DOI: https://dx.doi.org/10.1023/A:1018926631734]

30. Doha, E.H.; Rihan, E.A.; Hassan, M.I.; Kamel, N.M. Mono-implicit Runge-Kutta method for delay differential equations. J. Egypt. Math. Soc.; 2009; 17, pp. 213-232.

31. Butcher, J.C.; Chartier, P. A generalization of singly-implicit Runge- Kutta methods. Appl. Numer. Math.; 1997; 24, pp. 343-350. [DOI: https://dx.doi.org/10.1016/S0168-9274(97)00031-7]