Global dynamics of a mathematical model for a honeybee colony infested by virus-carrying Varroa mites

We establish a new four-dimensional system of differential equations for a honeybee colony to simultaneously model the spread of Varroa mites among the bees and the spread of a virus transmitted by the mites. The bee population is divided to forager and hive bees, while the latter are further divided into three compartments: susceptibles, those infested by non-infectious vectors and those infested by infectious vectors. The system has four potential equilibria. We identify three reproduction numbers that determine the global asymptotic stability of the four possible equilibria. By using Dulac’s criterion, Poincaré–Bendixson and persistence theory, we show that the solutions always converge to one of the equilibria, depending on those three reproduction numbers. Hence we completely describe the global dynamics of the system.

Colony collapse disorder (CCD) is the phenomenon of the majority of worker bees in a colony disappearing and leaving behind a queen, plenty of food and some nurse bees to care for the remaining immature bees. The phenomenon was renamed colony collapse disorder in 2006 when abnormally high die-offs (30-70% of hives) of common honeybee colonies have occurred in North America and at first no explanation could be given [6]. Several European countries have experienced the same phenomenon since 1998, and in the past few years countries in Africa and Asia also became affected by it. The reason seems to be a combination of factors, possibly including neonicotinoid pesticides or Israeli acute paralysis virus [7]. The collapse of honeybee colonies has become widespread in several regions of the world, and has been the subject of much discussion and research in recent years [8][9][10].
Varroa mite is an external parasitic mite that attacks the honey bees Apis cerana and Apis mellifera. It attaches to the body of the bee and weakens the bee by sucking fat bodies [11]. A significant mite infestation will lead to the death of a honey bee colony, usually in the late autumn through early spring. The Varroa mite is the parasite with the most pronounced economic impact on the beekeeping industry. Also, it is considered to be one of multiple stress factors [12] contributing to the higher levels of bee losses around the world. Varroa mites are carriers for many viruses that are damaging to bees, including Kashmir bee virus, sac-brood virus, acute bee paralysis virus, Israeli acuteparalysis virus, and deformed wing virus [13].
Several methods of treatment are currently applied to control these mites, which can be divided into chemical and mechanical controls. Usual chemical controls include "hard" synthetic chemicals such as amitraz, fluvalinate and coumaphos, while "soft" chemical controls (organic acids, essential oils) include thymol, sucrose octanoate esters oxalic acid and formic acid. Mechanical controls are usually based on disruption of some aspect of the mites' life cycle and they are generally intended not to eliminate all mites, but to keep the infestation at a tolerable level. Examples of mechanical controls include sacrifice of drone brood as Varroa mites most commonly attach to the drone brood, powdered sugar dusting which encourages cleaning behaviour and dislodges part of the mites, screened bottom boards which allows dislodged mites to fall through the bottom and away from the colony, brood interruption, application of heat to isolated brood combs or whole colonies and downsizing of the brood cell size. Another possibility in fighting the infestations is breeding more resistant colonies: several families of bees are able to coexist with Varroa mites (e.g. Africanized bees and Russian honey bees show a higher natural resistance against mites) [13][14][15].
A number of mathematical models have been established to model Varroa infestation and associated virus infections. Sumpter et al. [16] studied the disease in honey bee colonies and modeled the effects of Varroa mites on the brood and on the adult worker bees. The focus of the model was on the relationship between the mite population within a hive and its role in virus transmission within the hive. A study by Ratti et al. [17] examined the transmission of viruses via Varroa mites with the mites as vectors for transmission. In Ratti et al. [17], the population of honeybees was divided into honeybees that are virus free and that are infected with the virus, while the same population compartments of honeybees with adding mortality caused by mites in Ratti et al. [18].
The model presented in Kang et al. [19] follows Sumpter et al. [16] and Ratti et al. [18] with the same compartments for honeybees, but it also takes into account the fact that virus transmission occurs at different biological stages of Varroa mites and honeybees. The models of Ratti et al. [18] and Kang et al. [19] are extended in Ratti et al. [20] by introducing uninfected forager bees as a new honeybees compartment, also by adding different types of mortality.
The purpose of this work is to establish and analyse a novel mathematical model for the dynamics of a honeybee colony affected by Varroa mites by providing a new approach in comparison with earlier works: instead of considering separate compartments for the mites, we divide the honeybees into compartments depending on their infestation and infection status, based on earlier models for ectoparasite-borne diseases [21][22][23][24]. In the present paper we set up the simplest model that is complex enough to capture the fundamental features of the parallel transmission of the mite infestation and the infection. The model is based on the ectoparasite model given in [22], however, in the present model we consider different birth and death rates and additional disease-induced mortality for infected bees. Most importantly, we extend the model with a new compartment by differentiating susceptible hive bees and forager bees. Hence, our model consists of four nonlinear ODEs dividing the population of hive bees into susceptible (uninfected and uninfested) hive bees, hive bees infested by non-infectious vectors and hive bees infested by infectious vectors (infected hive bees) with mortality caused by mites. The fourth compartment stands for forager bees. We identify key variables that determine the collapse or survival of the bee colony, namely the severity of the disease and the rate of transmission, and examine different scenarios using different combinations of these variables.
The structure of the paper is the following. In Sect. 2, we establish our compartmental model describing the spread of the infestation and the disease. In Sect. 3, we study the three-dimensional subsystem and determine the reproduction numbers and prove the existence of the positive equilibria. In Sect.s 4 and 5, we investigate the local stability of equilibria and prove the persistence of the compartments, while, using the results of Sects. 4 and 5, we study the global stability of equilibria. In the last section, we use the results of the previous sections to study the two-dimensional subsystem of susceptible.

Mathematical model
Our mathematical model is based on the presence of a mite species which is a vector for a disease as well and transmitted to a susceptible host only upon adequate contact with an infested host. The honeybee population is divided into four compartments depending on the presence of the vectors and the disease transmitted by them as follows.
(i) Susceptibles: those who can be infested by the vector. Following [20,25] the healthy bee population is divided into hive bees and foragers. We denote by H s the compartment of susceptible hive bees and F s the compartment of susceptible forager bees. (ii) H m denotes hive bees infested by non-infectious vectors. (iii) H i denotes hive bees infested by infectious vectors, and thus infected with the disease.
In our model we assume that the following assumptions hold. infection. So, we assume that the death rate of infected hive bees H i is equal to d + δ where δ is the death rate caused by mites. 9. In accordance with [26], we assume that the proportion of recruitment of susceptible hive bees H s to become forager bees and the rate of healthy forager bees F s that are reverting to hive duties following social inhibition has the same proportion R(H s , F s ). We formulate this process of recruitment and social inhibition as a Holling-type II functional response, see [20,25,27].
Using the above assumptions, our model takes the form: where the term R(H s , F s ) in the first two equations represents the effect of social inhibition on the recruitment rate and is formulated as where the parameter σ 1 is the maximum rate at which hive bees are recruited as foragers when there are no foragers present in the colony. The term σ 2 F s H s +F s represents social inhibition, that is, the process whereby a surplus of foragers causes the foragers to revert to being hive bees. We assume that social inhibition is directly proportional to the forager population present in the colony.
In the present paper, following [16] and for technical reasons, we will consider a constant birth rate and study the special case The transmission chart of the model is depicted in Fig. 1. Let us introduce the notation S(t) = H s (t) + F s (t). With this notation, we can transcribe system (3) to the three-dimensional system Let us note that this system has a similar structure as the model given by Dénes and Röst in [21,22] and the rodent subsystem given by Dénes and Röst in [23], however, in the present case, in contrary to those papers, birth and death rates are different and there is additional mortality for infected hive bees. Because of these differences in the models and to make the present paper self-contained, we give all the proofs of the results on the dynamics of model (3). We also note that the methods applied in the analysis are also different: in the earlier papers, Lyapunov functions were used in the proofs of global stability, however, in the present paper we apply Bendixson-Dulac criterion instead.

Equilibria, reproduction numbers
To determine the equilibria of the full system (3), we start by calculating those of the subsystem (4). These are easily obtained by solving the algebraic system of equations We can determine four possible equilibria of system (4), one of which is disease-and infestation-free, one is disease-free with infestation, one is endemic where all vectors are infectious, and one is endemic where both infectious and non-infectious vectors are present: Due to the biological interpretation of the model, we are only interested in nonnegative equilibria. In the sequel, we say that a given equilibrium exists, if each of its three coordinates is non-negative.
We can determine four reproduction numbers by introducing a single infested (infectious or non-infectious) individual into a population in which neither infected nor non-infected parasites are present (E s ), or into a population where only noninfected parasites are present (E m ), or into a population where only infected parasites are present (E i ), and calculating the expected number of generated secondary cases.
If we introduce an individual infested by non-infectious mites into the infestationand disease-free equilibrium E s , we obtain the reproduction number by introducing an individual infested by infested mites into the same equilibrium we obtain the reproduction number Calculating the expected number of secondary infections caused by the introduction of an individual infested by infectious mites into a population in the equilibrium E m , we obtain the same reproduction number R 2 , as the transmission rate from H i individuals is the same for the S and H m compartment. Finally, let us introduce an individual infested by non-infectious mites into a population in the equilibrium E i . In this case, the expected sojourn time of an individual infected with the first strain in the H m compartment is (β 2 H * i + α + d) −1 , and the number of new infections generated by this individual is β 1 S * , where S * and H * i stand for the first, resp. Taking the product of these two expressions and substituting the values of S * and H * i at the equilibrium E i , we obtain the reproduction number In the next proposition we show how the reproduction numbers determine the existence of the four equilibria.
. All of the eigenvalues are negative if R 1 < 1 and R 2 < 1, while at least one of them is positive if R 1 > 1 and R 2 > 1.
(ii) If we linearize around the equilibrium E m , we find the eigenvalues λ H m 1 = λ S 1 , λ H m 2 = −λ S 2 and λ H m 3 = λ S 3 , thus we can argue similarly as in case (i). (iii) Linearization around the equilibrium E i yields the three eigenvalues . R 2 > 1 is needed for the existence of E i . If we add the terms in λ H i 1 , it is easy to see that the numerator of the fraction is the difference of the numerator and the denominator of the reproduction numberR 3 , which means that it is negative if and only if R 3 < 1. The absolute value of the term under the square root in the nominator of λ H i 2 resp. λ H i 3 is less than that of the first term which itself is negative as dα < Aβ 2 follows from R 2 > 1. Thus λ H i 2 and λ H i 3 always have negative real parts if R 2 > 1. (iv) Linearizing around E mi , we get the eigenvalues λ H mi 1 = −λ H i 1 , so it is negative if R 3 > 1 and λ H mi 2 = λ H i 2 and λ H mi 3 = λ H i 3 , from which the assertion follows.

Persistence
Before we can state our results on the persistence of the three compartments, we will need some notions and theorems from [28].
Definition 5.1 Let X be a nonempty set and ρ : System (4) generates a continuous flow on the feasible state space Proof To prove the persistence of S(t) we use the method of fluctuation (see for example Appendix A of [28]). Let S ∞ denote the limit inferior of S(t) as t → ∞.
From the fluctuation lemma we know that there exists a sequence t k → ∞ such that S(t k ) → S ∞ and S (t k ) → 0 as k → ∞. If we apply this to the equation for S(t), we obtain It is easy to see that for the total bee population we have lim t→∞ ( Using this and letting k → ∞, we obtain To show persistence of the infested compartments, we need some theory from [28]. We use the notation x = (S, H m , H i ) ∈ X for the state of the system and the usual notation ω(x) for the ω-limit set of a point x defined as ω(x) := {y ∈ X : ∃ t n such that t n → ∞ and Φ(t n , x) → y as n → ∞}.
which is positive if ε is sufficiently small as β 1 From R 2 R 3 > 1, we have A d β 1 > αδ+Aβ 2 d+δ , so we can write the following estimation for H m (t): This expression is positive for ε small enough as A d β 2 > d(d +α +δ) , which follows from R 2 > 1. This contradicts H m (t) → 0. Now, let us assume that M 2 is not weakly ρ-repelling. Then there exists a solution such that lim t→∞ ( and H m (t) > 0. Then for sufficiently large t, we have S(t) > d+α+δ For such t, we can give the following estimation for H m (t): which is positive if ε is sufficiently small as β 1 To prove the persistence of H i (t), we choose ρ(x) = H i (t). We have the equilibrium E s if R 1 ≤ 1 and the two equilibria E s and E m if R 1 > 1. We define the extinction space as Similarly as in the proof of the persistence of H m , is invariant, and M 1 and M 2 are isolated and acyclic.
Firstly, to show that M 1 is weakly ρ-repelling, let us suppose that M 1 is not weakly ρ-repelling. Then there exists a solution such that lim t→∞ (S(t) + H m (t) + H i (t)) = A d , 0, 0 and H i (t) > 0. Then for any ε > 0, for sufficiently large t, we have For such t, we can give the following estimation for H i (t): which is positive if ε is sufficiently small as β 2 and H i (t) > 0. Then for any ε > 0, for sufficiently large t, we have S(t) > d+α For such t, we can give the following estimation for H i (t): This expression is positive if ε is sufficiently small as follows from R 2 > 1, which contradicts H i (t) → 0.
We have shown uniform weak persistence in all cases, to show uniform (strong) persistence, we apply Theorem 4.5 from [28]. Our flow is clearly continuous, the subspaces X m , X i , X \X m and X \X i are invariant. The existence of a compact attractor is also clear, as all solutions enter a compact region after some time. This means that all conditions of [28,Theorem 4.5] hold and thus we obtain uniform strong persistence.

Global stability
In this section we extend the statements about local stability in the previous section to global asymptotic stability by means of the Bendixson-Dulac criterion and the Poincaré-Bendixson theorem, where we also apply the persistence results of the previous section. (2) For R 2 > 1 the following statements hold.
(i) If R 1 ≤ 1 and R 3 ≤ 1, then E i is globally asymptotically stable on X\X i and E s is globally asymptotically stable on X i . (ii) If R 1 > 1 and R 3 ≤ 1, then E i is globally asymptotically stable on X\X i and E m is globally asymptotically stable on X i . (iii) If R 3 > 1, then E mi is globally asymptotically stable on X\(X m ∪ X i ), E m is globally asymptotically stable on X i and E i is globally asymptotically stable on X m .
Proof Let us introduce the notation G(t) := S(t) + H m (t). With this notation, we can transcribe system (4) to the two-dimensional system This system has the two positive equilibria A d , 0 and To show that the limit of each solution of this system is one of these two equilibria, according to the Poincaré-Bendixson theorem, all we have to prove is that system (6) does not have any periodic solutions. To show this, we apply the Bendixson-Dulac criterion using the Dulac function D(G, In the first case, when R 2 ≤ 1, only the first equilibrium exists. Thus, in this case H i (t) → 0 and G(t) → A d as t → ∞, and therefore, the second Eq. of (6) takes the following form on the limit set: Clearly, for R 1 ≤ 1 (which is equivalent to γ ≤ 0), the solutions tend to zero on the limit set, therefore, for all solutions, H m (t) → 0 as t → ∞ and the limit set of any solution is the equilibrium E s .
It is easy to see that γ > 0 if and only if R 1 > 1. Thus, for R 1 > 1 we obtain on the limit set and by using the persistence of H m (t) we obtain that for all solutions, H m (t) → A d − α+d β 1 as t −→ ∞. From this follows the assertion of the first case of the theorem.
Secondly, for R 2 > 1 also the second equilibrium exists. However, we know from the previous section that H i (t) is persistent for R 2 > 1, thus for each solution of (6) started in X \X i we have lim t→∞ (S(t) + H m (t)) = G * . Thus, on the limit set the equation for H m (t) takes the form Similarly to the previous case, the nontrivial solutions of this logistic equation have the form It is easy to see that η > 0 if and only if R 3 > 1. Thus, for R 3 ≤ 1, lim t→∞ H m (t) = 0 and the limit of solutions started in X \X i is E i . In the case R 3 > 1 we have thus we obtain that the limit of solutions started in X \(X m ∪ X i ) is E mi . Solutions started in X m tend to E i . The limit set of solutions of Eq. (6) started in X m is the equilibrium A d , 0 . Thus, in this case the equation for H m (t) on the limit set has the form (7). Similarly to the previous cases, the nontrivial solutions of this equation have the form We have γ > 0 if and only if R 1 > 1. Thus, for R 1 ≤ 1, H m (t) → 0 (as t → ∞) and the limit of solutions started in X i is E s , while for R 1 > 1 we obtain the solutions started in X i tend to E mi . To complete the proof of the theorem, we notice that R 2 > 1 and R 3 > 1 imply R 1 > 1: and for R 1 > 1 we already established the global asymptotic stability of E m on X i . The proof of (iii) is complete.
The results of Theorem 6.1 are summarized in Table 1. In Figs. 2, 3, 4, 5 and 6, we denoted E s by black , E m by red , E i by blue and E mi by green.

The two-dimensional subsystem of susceptible bees
In this section we consider the two-dimensional subsystem of (3) consisting of the first two equations and make a discussion about the equilibrium points, local and global stability. Also, we assume that the subsystem (4) is in a steady state, and substitute any of the equilibria of the subsystem (4) into these to obtain the two-dimensional system where and H * m and H * i are the second, resp. third coordinates in any of the four equilibria (5).
We calculate the equilibria of the subsystem (9) by solving the algebraic system of equations This system has two equilibrium points: denoted by E 1 and E 2 respectively, with J = L 2 + 2σ 1 − 3σ 2 , and L 1 and L 2 are defined in (10). It can easily be seen that the first equilibrium E 1 is always positive and the second equilibrium E 2 is positive if J is negative.

Proposition 7.1
The local stability of the equilibria of (9) is determined as follows: Proof The statements of the proposition can be shown in a straightforward way by calculating the eigenvalues of the Jacobian of the linearized system at the equilibria. The Jacobian has the form (i) The eigenvalues of the Jacobian of the linearized equation around the equilibrium E 1 are λ 11 = −L 2 and λ 12 = − J 2 + 8σ 2 (J − σ 1 + 2σ 2 ). Both eigenvalues have negative real part, which means that the first equilibrium E 1 is locally asymptotically stable. (ii) If we linearize around the equilibrium E 2 , we find the eigenvalues λ 21 = −L 2 and λ 22 = J 2 + 8σ 2 (J − σ 1 + 2σ 2 ). Thus, the second equilibrium E 2 is unstable.
In the sequel, we want to find the limit set of each solution of H s and F s . So, let us consider F s = S * − H s and substituting this into (9) we get where Also, it is easy to prove that r 1 − r 2 > 0. The solutions of (11) take the form for c ∈ R + . Thus, we have Let us first consider the case H * i = 0, (i.e. when the subsystem (4) tends to the equilibrium E s or E m , which is equivalent to R 2 ≤ 1). In this case, if R 1 ≤ 1 the solutions of subsystem (4) tend to the equilibrium E s with S * = A d , L 1 = A, L 2 = d and So, we obtain If R 1 > 1 the solutions of subsystem (4) tends to the equilibrium E m with S * = d+α In this case, we get Thus, we obtain that all solutions of (9) tend to the equilibrium E 1 .
Similarly, in the case R 2 > 1, i.e. when H * i > 0 , the subsystem (4) tends to the equilibrium E i or E mi . If R 3 ≤ 1 the solutions tends to the equilibrium E i and if R 3 > 1 the solutions tends to the equilibrium E mi . Also, in this case we have that all solutions of (9) tend to the equilibrium E 1 .
Consequently, for all of the different cases all solutions of (9) tend to the equilibrium E 1 globally asymptotically stable.

Discussion
In this paper we established a novel compartmental model consisting of a fourdimensional system of differential equations for a honeybee colony to simultaneously describe the spread Varroa mites among the bees and the spread of the disease transmitted by the mites. Our work provides a different approach from existing models for Varroa infestation of honeybees as instead of separate compartments for the mites, we consider bee compartments, differentiated by the presence of infectious, resp. noninfectious parasites, hence, in contrast to earlier works, our model allows to keep track of infestation by non-infectious, resp. infectious mites separately and hence to give suggestions concerning the level control measures necessary to eliminate one or the other.
We calculated four potential equilibria of the system and three reproduction numbers. These threshold numbers determine which of the four possible equilibria of the system is globally asymptotically stable, depending on the parameters. By using persistence theory, the Bendixson-Dulac criterion and the Poincaré-Bendixson theorem, we showed that the solutions always converge to one of the equilibria, depending on the reproduction numbers. Consequently, we gave a complete description of the global dynamics of the system for all possible cases. Our model allows us to study the interplay between forager and hive bees, also to determine whether the recruitment and social inhibition rate will effect the stability of our system or not.
Our results suggest that to eradicate the disease, we have to decrease R 2 to be less than 1, which is possible by reducing the transmission rate β 2 or increasing the disinfestation rate α. If we also want to eliminate the infestation, then besides decreasing R 2 , we also have to decrease R 1 (possible by reducing the transmission rate β 1 or increasing α). The reduction of transmission rates is possible by breeding more resistant bees, while there are several methods to increase the disinfestation rate α, including the use of synthetic, resp. organic chemicals, as well as mechanical control measures e.g. heating, powdered sugar dusting or drone brood sacrifice. In case of the presence of the infection, decreasing only R 1 is not enough for the elimination of the infestation. The reproduction number R 3 is a threshold parameter which, in the case when the disease persists, shows whether all parasites become infectious.
We emphasize that the model given in this paper is the minimal model that is complex enough to include the essential features of the parallel transmission mechanism of the infestation and the infection. In this simplest model, which is a starting point for further research, we applied a simple birth function for technical reasons. It is a topic of a future work to make the model more realistic by including further compartments, different types of birth functions, as well as periodic parameters to account for seasonal variability.