A mathematical model for the spread of Varroa mites in honeybee populations: two simulation scenarios with seasonality.

We formulate and study a mathematical model for a honeybee colony infected with Varroa mites which describes the parallel phenomena of the spread of both the mites and the virus transmitted by them. We extend our previous model by including infected forager bees and considering model parameters as time-periodic functions. Firstly, we study the autonomous model and show the stability of equilibria. We present two simulation scenarios to study the impact of seasonality on the spread of Varroa mites and the disease they carry. Numerical studies are given to show how the parameter changes might lead to the colony's failure.

Introduction

The Western honeybee (Apis mellifera), first domesticated several thousands of years ago for honey production, is also the most important pollinator of food crops [1], [2], hence, it is not only ecologically [3] but also economically one of the most important insect species [4]. In recent years, honeybee colonies have been menaced by several factors including climate, pesticides and infections, resulting in high losses in honey bee populations [5]. Colony collapse disorder, consisting of a sudden colony death, and a lack of healthy adult bees in the hive has been observed in recent years in several countries of the world, however, the reason of the phenomenon could not be completely identified yet. Several factors may contribute to this disorder, including various infections, among them Varroa mites carrying acute bee paralysis virus, Kashmir bee virus, sac-brood virus, Israeli acute paralysis virus and deformed wing virus.

Varroa mite is an ectoparasite mite that infests honeybee colonies. The life cycle of Varroa mites is tightly adapted to the development of honeybees. Varroa mites can only reproduce in a honey bee colony. A female mite gets into the brood cell before it is capped, feeds on the juvenile bee and it also reproduces [6]. The adult female mite becomes attached to the bee and weakens it by sucking its fat bodies [7]. Varroa infestation, in high levels, may result in the collapse of the colony, which happens most often between the end of fall and the beginning of spring [8]. Mites can spread both vertically and horizontally in the bee population. Vertical transmission of the mites happens when the bee colonies swarm out for reproduction and the attached mites move to the new nest site with the swarming bees. Horizontal transmission between colonies occurs mainly via migrating worker bees to other colonies and also when stronger colonies rob honey stores of weaker colonies as well as when infested brood or bees are moved by beekeepers. A further way of horizontal transmission involves mites moving from one colony to another by attaching to foragers from a mite-free colony after having been removed by an infested bee onto a flower [9]. Varroa destructor mites are capable of rapidly infesting honey bees foraging at a feeder or on flowers of various species. Mites can rapidly ascend foraging honey bees, and although the bees try to repel the mites, in most occasions they succeed to leave the foraging site still attached to the bee. The transmission of Varroa mites from flower to bee can occur within just two seconds spent on a flower foraging [9]. Varroa mites are the parasites with the most important economic burden on beekeeping. It is also regarded as one of several stress factors [10] responsible for the increasing degree of bee loss in several countries of the world.

Several mathematical models have been established to describe the evolution of honeybee populations and to model the spread of infections affecting them. In [11], Sumpter et al. established a model to describe the impact of the parasites on the brood and on adult bees. Ratti et al. [12] studied virus transmission via Varroa mites, with the bee population divided into virus-free and infected individuals. The same bee compartments were considered with mite-induced mortality [13]. Kang et al. [14] also considered that transmission of the virus happens at different biological stages of honeybees and mites. In [15], the authors extended earlier models by considering uninfected forager bees and adding different types of mortality. Petric et al. [16] studied a model for a colony infected with Nosema ceranae. The model with time-periodic coefficients includes healthy and infected hive and forager bees, as well as disease potential in the hive. Recently, in [17], the authors formulated a model for a honey bee colony infected with the parasite Nosema ceranae. The model includes healthy and infected hive and forager bees, environmental spores of N. ceramae and also considered sugar and food storage. Numerical analysis of the model revealed that parasite extinction, co-existence between bees and the parasites, or colony failure are all possible. In [18], a compartmental model for the dynamics of honey bee colonies was developed to study the effects of pollen on colony failure. Seasonal changes affecting the food amount collected by foragers are also considered. Variations in the availability of pollen and nectar within a year were represented by trigonometric functions. In [19], interaction among brood, adult bees and mites was considered in which the time from brood to adult bee is also taken into account. The authors proposed a single-patch delay equation model with stage structure to investigate the effects of age structure and parasitism on colony dynamics under periodically changing circumstances. Seasonality is reflected in the model by the seasonally varying egg-laying rate of the queen bee. The time-dependent egg-laying rate of a fully mated queen is represented by a trigonometric function.

In [20], we took a novel approach by considering two types of infested bees: those infested by mites not carrying the virus and those infested by virus-carrying mites. Also, the bees were divided into the subpopulations of hive bees and forager bees. In that work, to enable a complete analytical characterization of the global dynamics of the model, we made the simplifying assumption that only hive bees could be infested. The aim of the present work is to extend that model by allowing forager bees to be infested by mites. Further, to make the model more realistic, we introduce time-periodic parameters following [16], [17], [18], [19]. First, we study the autonomous case without seasonality and then prove the existence and stability of the equilibria. We present two simulation scenarios to study the main colony dynamics and to show how seasonality affects the spread of Varroa mites in honey bee populations. In the first scenario, we assume piecewise constant function based on seasonal averages for each time-dependent parameter, while in the second scenario, we assume that the queen's egg-laying rate, as well as the transmission rates are continuous periodic functions represented by trigonometric functions. Finally, we use numerical simulations to investigate what kind of parameter changes might lead to colony failure.

Mathematical model

Our compartmental mathematical system is established based on the presence of a mite species, which is also a disease vector and can only be transmitted to a susceptible host through enough contact with an infested host. Varroa mites are considered vectors in this study, and the term “vector” refers to a Varroa mite throughout the text. By an infested honeybee, we mean a honeybee with (virus-free or virus-carrying) parasites, while we talk about infection in the case when the bee is infested by a mite carrying the virus, thus infecting the bee. Depending on the presence of the mites and the infections transmitted by them, the honeybee population, denoted by B, is split into six compartments as follows: Consequently, the total number of hive bees, total number of forager bees and total honeybee population, respectively, are given by

Healthy bees: those who can be infested by the vector. Healthy bees are classified as either healthy hive bees () or healthy foraging bees ().

Hive bees () and foraging bees () who are infested by vectors not carrying the virus.

Hive bees () and foraging bees () who are infested by vectors virus-carrying, and thus infected with the disease.

Model assumptions

We extend our earlier mathematical model of the disease in Dénes and Ibrahim [20] by including infected forager bees and taking into consideration model parameters as time-periodic functions. We follow the assumption given in [16], [17], [20] and add new items to describe the new terms. Hence, the following assumptions will be necessary for our model:

Hive bees () and foraging bees () infested by virus-carrying vectors can transmit the mites to healthy bees ( and ).

We assume that at least one Varroa mite not carrying the virus is attached to every single individual in and compartments.

Hive bees () infested by virus-carrying vectors can transmit the disease to healthy bees ().

Foraging bees () infested by virus-carrying vectors can transmit the disease to healthy bees ().

We assume that at least one virus-carrying Varroa mite is attached to every single individual in and compartments.

A bee infested by virus-carrying vectors can transmit the infection to a bee infested by virus-free vectors, i.e. an individual from the compartment can move to compartment following an adequate contact with a bee l from compartment ( or ). Similarly, a bee from the compartment can move to compartment after contacting a member of compartment ( or ).

We suppose that disinfestation rates from both compartments and to the susceptible compartment are equal and denoted by α.

The hive bees' death rates in the spring, summer, and fall months can be neglected in comparison with the rate at which hive bees are recruited to foraging work, as reported in [16], [17].

In accordance with [9], Varroa mites can move easily and fast from flowers to foraging honey bees, infecting them outside the hive. However, in our model, we only consider infection from bee to bee, with no infection from outside.

Following [13], [16], [21], [22], we suppose that the queen does not die or is replaced with no negative impact on the performance of the colony. We further assume that the disease has no effect on the queen bee.

Without taking into account larval and pupal stages, following [11], [13], [16], [22], we assume bees emerge as adult worker bees. The eclosion rate for hive bees is determined by the product of the daily maximum potential eclosion rate (the average number of eggs laid per day by the queen) and a measure of the colony's brood-rearing capacity (a sigmoidal function of the total colony population), as in [22], [23]. As a result, low worker bee populations can hamper the eclosion rate.

We assume our parameters to be time-dependent, with a one-year periodicity, as in [13], because both the biology and population dynamics of honey bees, as well as the dynamics of Varroa mites, fluctuate with the seasons.

Model equations

Our model is based on the bees model presented in [20], with various modifications and extensions. In accordance with the transmission diagram in Fig. 1 and using the above assumptions, the corresponding system of differential equations takes the form where

Figure 1

Transmission diagram showing honey bee colony dynamics combined with the dynamics of infestation by virus-carrying or virus-free Varroa mites. The progression of infection is represented by black solid arrows. Disinfestation is shown by blue dashed arrows. The movement here between forager bee and hive bee compartments is represented by solid blue arrows. Birth and death functions are shown by solid green arrows.

The maximum emergence rate of healthy hive bees is given by where is the queen's egg laying rate. The term represents brood maintenance function, in [17], [23], with half-saturation constant K. To incorporate the seasonal changes all parameters are assumed to be time-dependent. More precisely, we assume that these parameters are periodic with a periodicity of one year.

In system (1), is the maximum rate at which hive bees switch to become foragers when no foragers are present in the colony. The parameter stands for a reduction in recruitment of foragers and possible reversion if the distribution of hive and forager bees is perturbed [16]. In accordance with [24], we assume that the recruitment of hive bees to foragers is proportional to the actual number of hive bees. Social inhibition is the phenomenon of a surplus of foragers causing foragers to change back to become hive bees. Our assumption is that both healthy and infected bees have the same rates of recruitment and social inhibition. The description of the model parameters is listed in Table 1.

Table 1

Description of parameters of model (1).

Parameters	Description
E(B,t)	maximum emergence rate
K	brood maintenance constant
β1(t),β2(t),β3(t)	transmission rate from hive bee-to-hive bee
β4(t),β5(t),β6(t)	transmission rates from forager bee-to-forager bee
β__FH(t)	transmission rate from forager bee-to-hive bee
β__HF(t)	transmission rate from hive bee-to-forager bee
α	disinfestation rate from (Hm and Hi) to Hs
σ2(t)	rate of reversion
σ1(t)	rate of recruitment
ds(t),dm(t),di(t)	hive bees death rates
μs(t),μm(t),μi(t)	forager bees death rates

To address the existence and uniqueness of the solutions of system (1), we introduce some notations Define

For any, every forward solution of(1)eventually enters Ω for all positive t and is bounded from above by where .

By the Lipschitz condition we can prove the existence and uniqueness of model (1) solutions. From Equation (1), we have Applying the comparison theorem [25], we complete the proof. □

Analysis in the autonomous case

In this section, we study the stability of the autonomous model solutions obtained from Equation (1) by setting the time-dependent parameters to be constant. Although this is an impractical simplification in quantitative terms, it gives a qualitative understanding of the interplay of several parameters of the model. These results would also lead to the formulation and explanation of the full time-dependent model numerical simulations.

From Equation (1), we obtain

The following proposition is analogous with [15, Proposition 3.1] and [16, Proposition 3.3]. This result indicates that a sufficient number of healthy worker (hive and forager) bees are necessary to preserve the hive in order for it to establish itself as a properly working colony.

Existence of equilibria

We defineThe trivial equilibriumof the model(2)always exists. Ifthen the trivial one is the only equilibrium. If the above inequality is reversed, then two positive equilibriaandexist. Furthermore,

From Equation (2), it is easy to see that is always an equilibrium. In the absence of mites and viruses, we have From (4) we get By solving Equation (5) with respect to , we obtain We ignore the lower branch to keep the positivity of the disease free-solution and obtain Adding Equations (3) and (4), we get Substituting Equation (6) into Equation (7), we find that is a positive root a polynomial of degree n. With Descartes' rule of signs, we obtain that, depending on the parameters, there are either none or two such positive roots. We have and for sufficiently large H. Positive roots exist if the maximum of f is positive. The maximum is reached for Thus, we have If then there exist two positive disease-free equilibria and with . □

Stability of trivial equilibrium

The empty-hive equilibriumof model(2)is locally asymptotically stable.

By the tangent criterion, it can be seen that positive initial values yield non-negative solutions. We define .

From (2), we have Consider the following auxiliary equation: Applying the comparison theorem [25], we obtain with initial value .

For given , the set is positively invariant. Equation (8) has the asymptotically stable equilibrium for . Therefore, for small enough , we get Thus, for small enough all solutions starting in will tend to the equilibrium . □

Stability of disease-free equilibrium

Assume is a disease-free equilibrium of (2) in accordance with Proposition (3.1) . Then is asymptotically stable if the inequalities are satisfied.

Linearizing system (2) around the disease-free equilibrium , we obtain the Jacobian where The eigenvalues of the matrix J are negative if Then, starting with Equation (10) and substituting in Equation (11), we have Therefore, the inequalities in Equation (9) are satisfied, and hence the proof is complete. □

These results indicate that the disease-free equilibrium is stable if inequalities (9) are satisfied. If, however, the disease-free equilibrium is unstable, it is unclear whether the system will converge to the trivial equilibrium or an endemic equilibrium. We note that these findings are in accordance with the results concerning the stability of the disease-free equilibrium in [15], [16], [20].

Simulation results: non-autonomous model

In this section, we present some numerical experiments concerning model (1). Following [16], we generally assume that the Hill exponent n takes a value . Thus, if , emergence is close to maximal rate, while for a very low number of new bees emerge. In [23], the emergence rate's sigmoidal form had simple Holling type 2 saturation, implying that emergence is proportional to the number of worker bees at low bee populations. Since bees are eusocial insects with colonies that cannot survive in the long term with extremely low population sizes, we conceive that Holling type 3 form is ecologically more realistic, as in [16].

Using the Piecewise function of Wolfram Mathematica software, we create piecewise constant function based on seasonal averages for each time-dependent parameter, except where stated differently. Parameter values given in Table 2 and Table 3 represent average value of each season. In particular, we assume that when for spring, summer and fall, while in winter we set (all forager bees rapidly revert back into the hive) and we maintain , as in [16], [17], [24], [26]. We assume that hive bees die at the same rates equal to zero for spring, summer and fall, while in winter these values are assumed to be . For the foragers, following [16], [17], [24], [26], we assume that death rates of healthy bees, those infested by non-infectious vectors and infected bees are , respectively, for spring, summer and fall, while in winter the death rates are assumed to be as foragers return to the hive.

Table 2

Seasonal averages of some of the model (1) parameters, derived from the data in [11], [13], [16], [22], [24], [27], [28], [29] and used in scenario 1.

Parameters	Spring	Summer	Fall	Winter	References
E˜(t)	500	1500	500	0	[11], [16], [24]
K(t)	8000	12000	8000	6000	[13], [16], [24]
βi(t),i = 1,…,6	5 × 10−6	5 × 10−6	5 × 10−6	5 × 10−6	estimated
β__FH(t)=β__HF(t)	5 × 10−6	5 × 10−6	5 × 10−6	5 × 10−6	estimated
σ2(t)	1.5	1.5	1.5	1.5	[16], [24]
σ1(t)	0.25	0.25	0.25	0	[22], [24]
ds(t),dm(t),di(t)	0	0	0	0.00649	[16], [27]
μs(t),μm(t)	0.08511	0.08511	0.08511	0	[16], [28]
μi(t)	0.16936	0.16936	0.16936	0	[16], [28], [29]

Table 3

Seasonal averages of some of the parameters in model (1), derived from the data in [13], [16], [22], [24], [27], [28], [29] and used in scenario 2.

Parameters	Spring	Summer	Fall	Winter	Reference
K(t)	8000	12000	8000	6000	[13], [16], [24]
β˜_i,i=1,…,6	10−6	10−6	10−6	10−6	estimated
β˜__FH=β˜__HF	10−6	10−6	10−6	10−6	estimated
σ2(t)	1.5	1.5	1.5	1.5	[16], [24]
σ1(t)	0.25	0.25	0.25	0	[22], [24]
ds(t),dm(t),di(t)	0	0	0	0.00649	[16], [27]
μs(t),μm(t)	0.08511	0.08511	0.08511	0	[16], [28]
μi(t)	0.16936	0.16936	0.16936	0	[16], [28], [29]

The results of numerical simulations of two different scenarios that illustrate the main dynamics of the bee colony are presented below. The base case with no disease in the system is shown in each scenario. In addition, we show three other situations in which the bee population has been infested by non-infectious vectors, infectious vectors, or both. Except where otherwise noted, our simulations began on the first day of spring in year 1, with initial conditions In all cases, the year is assumed to begin with spring and that the length of all four seasons is the same, 91.25 days. Introducing e, spring occurs from 0 to l, summer from l to 2l, fall from 2l to 3l and winter from 3l to 4l days. Simulations are performed for a period of 3,650 days (approximately 10 years), except in the disease-free case, which is obtained for a period of approximately 3 years.

Scenario 1: piecewise constant birth and transmission functions

In this scenario, we integrate the system (1) numerically with piecewise constant birth and transmission functions. We assume that the maximum emergence rate of healthy adult hive bees is set to be 1500 for summer, 500 for spring and fall, and 0 for winter. The transmission rates are assumed to be constant for each season. The seasonal average values of the parameters are given in Table 2. Fig. 2 depicts the curve of the queen's egg laying rate that will be used in this scenario.

Figure 2

The curve of the queen's egg laying rate with seasonal values is shown in Table 2.

Case 1: disease-free case

In the absence of the disease, we first show possible model simulation outcomes. The bee populations fluctuate with the seasons, as shown in Fig. 3, and the colony rapidly settles into a healthy periodic solution. The number of bees in a hive decreases during the winter, when no new eggs are produced and there are no new adults emerging. Hive bees are recruited for foraging in the spring, and the colony strengthens during the summer. It begins to drop again in the fall, when the no eggs are laid and foragers return to the hive.

Figure 3

Disease-free periodic solution of model (1) with parameter values given in Table 2 and α = 0.1237.

Case 2: the bee colony infested by virus-free vectors only

We assume that the colony is only infested by virus-free vectors. Also, we increase the level of to 10 with a small change in the disinfection rate. Fig. 4 shows that the hive and forager bees ( and ) infested by virus-free vectors are presented in the bee population and fluctuate with the season.

Figure 4

Solution of model (1) with parameter values given in Table 2 and α = 0.145.

Case 3: the bee colony infested by infectious vectors only

In this case, we assume that the colony is only infested by infected vectors. Thus, we increase the level of to 10 and increase to . We obtain that the bees infected by infectious vectors ( and ) periodically reappear in the colony, as shown in Fig. 5.

Figure 5

The model (1) solution with parameter values given in Table 2 except β2 = 7.8 × 10−6 and α = 0.15.

Case 4: endemic periodic solutions

To obtain the full periodic solution where the populations fluctuate with the seasons, we increase the level of and to 10 and increase to with a slight decrease in the infestation rate from both and . Fig. 6 shows the bee colony is infested by both virus-free and virus-carrying vectors and the disease is endemic in the whole bee population. Fig. 5 shows that the bee population is infested by virus-carrying mites only, while Fig. 6 shows that the bees are infested by both virus-free and virus-carrying mites, and the disease is endemic in the bee population (hive and forager). The hive does not die out in Fig. 5, but it does not attain an endemic periodic solution. Fig. 6 illustrates the disease being controlled to the point that the hive does not die out but tends to an endemic periodic solution. There are only a small number of infected bees, but it is still sufficient for the disease to reemerge each year.

Figure 6

The model (1) solution with parameter values given in Table 2 except β2 = 6.56 × 10−6 and α = 0.128.

Scenario 2: continuous birth and transmission functions

In this scenario, we assume that the queen's egg laying rate () as well as the transmission rates ( and ) are periodic continuous functions.

Based on earlier works where trigonometric functions were applied for seasonal parameters (see e.g. [19], [30], [31]), these time-dependent parameters take the form , , and where for . Here, denotes the day of the year with the maximum egg-laying rate, is the baseline egg-laying rate, are the baseline transmission rates, b is the seasonality parameter, and t is time measured in days. The seasonal values of the parameters used in this scenario are given in Table 3. The curves of the queen's egg laying rate and the transmission rate that will be applied in this scenario are shown in Fig. 7.

Figure 7

The curve of in (a) the queen's egg laying rate where E1 = 2000 and E2 = 320, and in (b) transmission rates, for example, where and b = 120.

In this scenario, we show also four different cases similar to Section 4.1:

Case 1: disease-free case

Similarly, the bee populations fluctuate with the seasons when no infected bees are introduced, and the colony rapidly reaches an infestation-free periodic solution, as can be seen in Fig. 8.

Figure 8

The model (1) solution with parameter values given in Table 3 and α = 0.07.

Case 2: the bee colony infested by non-infectious vectors only

In this case, we assume that the colony is only infested by non-infectious vectors and increase the level of to 10. Fig. 9 shows that the hive and forager bees ( and ) are introduced in the bee population and fluctuate with the season.

Figure 9

The model (1) solution with parameter values given in Table 3 except , α = 0.13, and .

Case 3: the bee colony infested by infectious vectors only

Here, we assume that the colony is only infested by infected vectors. Hence, we increase the level of to 10 and increase and to . The infected by infectious vectors ( and ) are presented in the colony, as shown in Fig. 10.

Figure 10

The model (1) solution with parameter values given in Table 3 except , α = 0.075, and .

Case 4: endemic periodic solutions

To show the endemic periodic solution, we increase the level of and to 10, and rise to and to . Fig. 11 shows the bee colony is periodically reappearance and the disease is endemic in the bee populations.

Figure 11

The model (1) solution with parameter values given in Table 3 except , , and α = 0.134.

Seasonal average and annual average values for the piecewise constant functions and the trigonometric functions applied in scenario 1 and scenario 2 are calculated and compared in Table 4. The seasonal average observed in the parameter (the queen's egg laying rate) by interpolating the piecewise constant seasonal averages and the trigonometric functions are shown in Fig. 12. The seasonal average for by using the trigonometric functions is higher than the seasonal average of the piecewise constant function for all seasons, while the seasonal average for (the transmission rate) using the trigonometric functions is smaller than the seasonal average of the piecewise constant function for all seasons as shown in Fig. 12. The annual average for by using the trigonometric functions is higher than the annual average from the piecewise constant function, while the annual average for using the trigonometric functions is smaller than the annual average from the piecewise constant function as given in Table 4.

Table 4

Comparison of the average values of the parameters computed in Scenario 1 and Scenario 2 and shown in Fig. 12.

Parameter	Function	Spring	Summer	Fall	Winter	Annual average
E˜(t)	Piecewise	500	1500	500	0	625
	Trigonometric	581.19	1809.38	600.56	8.85	750
β(t)	Piecewise	5 × 10−6	5 × 10−6	5 × 10−6	5 × 10−6	5 × 10−6
	Trigonometric	4.38 × 10−9	5.54 × 10−14	9.16 × 10−8	6.08 × 10−7	1.76 × 10−7

Figure 12

The curve of in (a) the queen's egg laying rate used in scenario 1 and scenario 2.

Colony failure

The colony fails if, after a given number of years, depending on initial conditions and parameters, the colony is not strong enough to recoup at the end of the winter season and goes extinct owing to wintering losses, i.e. . In this section, we show numerical simulation to study what kind of parameter changes might lead to the colony failure. We study the impact of the transmission of the disease between individuals by contact, for example, increasing the transmission rates between hive bees and hive-forager bees. We assume the base case as the disease free periodic solution is eventually obtained and . Moreover, we assume that emergence is proportional to the population size of worker bees at low bee populations and considering the emergence rate's sigmoidal form had simple Holling type 2 saturation.

We first examine both scenarios by increasing the direct transmission from infected hive bees to healthy hive bees. By increasing the value of in scenario 1 to , we observe that after two years the number of bees decreased but still not extinct, as shown in Fig. 13(b). For the number of bee populations tend to zero and goes extinct after three years, see Fig. 13(c). Using the same parameters and values in scenario 2 with similar changes in the value of , as given in Fig. 14, we obtained the extinction of the bee population after four years and the colony fails if . Fig. 13 and Fig. 14 show that increasing the transmission rate from infected hive bees to healthy hive bees due to direct contact between individuals might lead to the colony failure after three or four years.

Figure 13

The curve of the total bee population size when β2(t)=5 × 10−6, β2(t)=3.5 × 10−5 and β2(t)=5 × 10−5 with parameter values given in Table 2.

Figure 14

The curve of the total bee population size when , and with parameter values given in Table 3.

To investigate the effect of the disease transmission between hive and forger bees on colony strength, we repeat the simulation for selected values of transmission rate from hive to forager bees () and transmission rate from forager to hive bees (), the rest of parameters values the same as used in the base case for both scenarios (see Table 2 and Table 3). For the total number of the bee population decreased to after two years. As increases further, e.g. at , the colony fails after three years that is shown in Fig. 15. For selected values of we plot in Fig. 16 the total bee population, we observed that the colony losses it's strength and fails at after three years. Fig. 15 and Fig. 16 indicate that if the disease is endemic in the hive due to a significant number of diseased hive bees, who subsequently spread the disease to forager bees, the colony will collapse after three years in both cases.

Figure 15

The curve of the total bee population size when , and with parameter values given in Table 2.

Figure 16

The curve of the total bee population size when , and with parameter values given in Table 3.

Varroa mites can jump from flowers to foraging honey bees, infecting them outside the hive [9]. As a result, if a high number of infected foraging bees return to the hive, the population of infected hive bees will grow. Similarly, we run the simulation again with different values of to assess the effect of disease transmission rate from forager to hive bees. For chosen the colony fails as shown in Fig. 17(c), while in Fig. 18(c) the colony fails at . The colony fails after two to three years, as can be seen in Fig. 17 and Fig. 18, due to a large number of infected foraging bees returning to the hive and/or forager losses.

Figure 17

The curve of the total bee population size when , and with parameter values given in Table 2.

Figure 18

The curve of the total bee population size when , and with parameter values given in Table 3.

As can be seen in Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, increasing one of the transmission rates () up to a certain level causes the colony to lose strength and become too weak to rebound at the end of the winter, potentially leading to colony failure after two to five years.

Discussion

Honeybees serve a significant role in the production of many agricultural products and in the preservation of plant variety in undisturbed ecosystems. For more than a decade, honeybee colonies have seen a significant decline in bee populations. This collapse may not have a single cause, rather there are a variety of factors at work [32], [33]. Colony failure has been strongly linked to the presence of Varroa mites [14], [34]. In this article, we developed a novel non-autonomous mathematical model to investigate the spread of Varroa mites and the viruses they carry in a honeybee colony. The main novelty of the model is that it describes the parallel spread of the virus and that of the mites carrying them and transmitting them to honeybees by introducing compartments for honeybees infested by virus-free and virus-carrying mites. Seasonality is also integrated into our model because the temporal dynamics of honeybee colonies are affected by seasonality (temperature, photoperiod, etc.). This model extends our recent mathematical model [20] of the disease by including infected forager bees, and model parameters are simulated as time-periodic functions as well.

We initially investigated the autonomous case without seasonality and demonstrated the existence and stability of the equilibria. Since the disease, its viability, and its level of infection undergo a drastic change throughout the year, we consider a non-autonomous system where coefficients are assumed to be periodic. The resulting model is algebraically rather complex, hence, there is no possibility to perform rigorous mathematical analysis via Floquet theory and further usual tools from the theory of ordinary differential equations. For this reason, the dynamical behaviour of the model is investigated numerically in computer simulations.

Besides piecewise constant parameters, often used to describe the effect of seasonal variation of the parameters in the spread of mites and the viruses they carry [12], [13], [15], [16], [17], [23], [24], we also used periodic functions. This way, two simulation scenarios are presented to investigate the main dynamics of the bee colony and to show how seasonality affects the spread of Varroa mites. Four different possibilities are provided to the reader in both scenarios. In each scenario, the base case with no disease in the system is shown, see Fig. 3 and Fig. 8. Moreover, the system converges to the disease-free periodic solution. In the last three cases, we investigated the impact of Varroa mites on bee populations in which virus-free vectors, virus-carrying vectors, or both have infested the bee population, see Figure 4, Figure 5, Figure 6 and Figure 9, Figure 10, Figure 11. The system converges to the endemic periodic solution only in one case as shown in Figure 6, Figure 11.

Depending on initial conditions and parameters, the colony fails if it is too feeble to rebound at the end of winter and goes extinct after a certain number of years, which concurs with the findings in [22], [24], [26], [35]. We performed numerical simulation to demonstrate what kinds of parameter modifications could result in colony failure. Our results from both scenarios indicate that if the transmission rate from hive bees infested by virus-carrying vectors to healthy hive bees reaches a specific threshold, the colony will collapse within three to four years. The interplay between hive and forager bees has also contributed to colony collapse. Furthermore, if one of the transmission rates from hive to forager bees or forager to hive bees exceeds a certain level, the colony loses strength and collapses within two to three years. Mite control is critical for the colony's survival. As a result, our recommendation to beekeepers is to lower transmission rates while also increasing disinfestation rates, both of which have a substantial impact. To control these mites, a variety of treatments are now being used, which may be classified into chemical and mechanical controls [36], [37], [38].

Comparing the results of the two simulation scenarios, one may find that depending on the parameters, the risk of colony collapse is might be shown smaller when using piecewise constant functions (see Fig. 13(c)), and might be shown smaller when using periodic functions 17(c). At the same time, it is important to note that without a secure knowledge on the exact shape of the functions describing the time-dependent parameters, it is a very difficult question to judge which functions reflect reality the most. However, we may conclude that the choice of non-autonomous parameters affects the time to colony failure predicted by the model.

Our model has several possibilities for further development. We extended the model to account for increased forager losses by including the effects of foragers being exposed to pesticides while foraging on crops treated with such chemicals. This offers us a better understanding of how these effects interact to affect a colony's survival. A further possibility is to study how Varroa mite infection affects food storage dynamics and how that affects colony survival as well as yearly honey supply.

Declarations

Author contribution statement

M.A. Ibrahim, A. Dénes: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This research was supported by the ÚNKP-20-3-New National Excellence Program of the Ministry for Innovation and Technology from the source of the . M.A.I. was supported by the Stipendiary Hungarian scholarship and a scholarship from the of the Arab Republic of Egypt. A.D. was supported by the Hungarian grants NAIF PD 128363 and NAIF KKP 129877.

Data availability statement

No data was used for the research described in the article.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.