Impact of Infection on Honeybee Population Dynamics in a Seasonal Environment

.


Introduction
The bee colony develops according to a seasonal cycle subject to climatic variations and the influence of the environment in which the bees are located.In temperate zones, the activity of the colony is subject to the effect of the four seasons which punctuate its development.The biological cycle of the colony is regulated by the laying of the queen which starts more or less early at the end of winter and the beginning of spring depending on several parameters including the ecotypes of bees and the climate. .The peak population of bees in the colony is reached in May-June and their number gradually declines from July.At the end of the season, egg-laying depends heavily on late-season temperatures before arriving at overwintering where a cluster of winter bees protects the queen.These winter bees are biologically different from summer bees with more developed fat bodies, a higher level of proteins in the hemolymph and in the hypopharyngeal glands, a lower level of juvenile hormone.These differences are longevity traits.The colony builds up its winter bee reserve as soon as the queen's egg laying declines.This is a kind of internal hormonal regulation.Less brood, for example, allows nurses to store vitellogenin (longevity protein based on royal jelly) in their fatty bodies.Naturally, seasonal climatic variations have an impact on the development of the colony.Cold waves in spring when the queen's clutch is developing or warmer late seasons which lead to continued egg laying have consequences on the general development of the colony.Climatic peaks, characterized by brutal and extreme episodes (frost, rain, drought), are not without consequences on the life of the colony in the sense that they disrupt its development, influencing the entry of pollen and nectar, modifying the investment of the energy spent and facilitating the development of certain pathogens.The dynamics of the bee population in the colony is more or less directly subject to these climatic hazards.A honey bee colony is a population of related, interacting individuals.We are faced with a very complex society influenced by complex population dynamics.Individuals in the superorganism have roles or assignments that evolve within the group throughout their lives.Their lifespan is strongly influenced by their role.We know that the division of labor between the population of workers depends on the age and needs of the colony.According to a general principle, young workers are mainly used in the hive for colony maintenance and brood care tasks such as feeding.It is only the older workers who are responsible for supplying the colony with food and who have contact with the external environment via foraging.This principle responds to a development process based on utilitarian social behavior.This explains why this great principle can be called into question depending on the needs of the colony.If the foragers see their numbers drop drastically under the effect of environmental stress for example, the other workers in the colony will accelerate their behavioral development to enter into a dynamic of compensation: they will therefore forage early and undoubtedly die early.Conversely, if there is an overabundance of foragers and a lack of nurses, the mental behavioral development of some foragers may regress and they may reassume the role of nurses.We talk about social inhibition.All of these behavioral adaptations are governed by a now well-identified pheromonal mechanism.A clear interaction exists between the assignment of colony workers and their longevity.If workers start foraging early to compensate for the lack of foragers, their lifespan may be reduced and the time spent on brooding is also reduced, which can have a significant impact on the growth of the colony.Likewise, the various well-known stresses (diseases, varroa, etc.) can impact the growth of colonies and lead to the weakening of worker populations.The widespread collapse of honey bee colonies has been the subject of much discussion and research in recent years.Aside from their ecological importance, honey bee populations have a large economical impact on agriculture in North America, Europe, the Middle East, and Japan.The focus of research has been largely on environmental factors outside the hive, such as pesticides or insecticides, which may cause death or injury to foraging bees and jeopardize their return to the hive.The reduced number of foraging bees then leads to younger hive bees being recruited prematurely to perform foraging duties and this chain reaction ultimately leads to a disruption in the dynamics of the colony as a whole.A key element in this category of disruption to honey bee population dynamics is the untimely death of a certain proportion of foraging bees outside the hive and the consequences of this on the colony as a whole.
An important question here concerns the threshold in the death rate of foraging bees that would determine the survival or collapse of the bee colony.
In the present paper we consider a different category of disruption to the healthy dynamics of a bee colony in a seasonal environment, namely one in which the key hazard is an infection by a communicable disease acquired by foraging bees outside the hive.The key difference here is that foraging bees that have been infected would then transport the disease into the hive and go on to infect other members of the colony within the hive.Here too the affected bees will ultimately suffer an untimely death, but the effects on the dynamics of the colony are clearly more complex because the infection in this case may now involve all members of the colony.
The goal of this paper is to consider the influence of the seasonality on the spread of disease within a bee colony with the underlying demographic dynamics of the colony.The basic reproduction number, R 0 , was defined by using an integral linear operator.We perform the global analysis of the proposed system.It is deduced that the disease-free solution is globally asymptotically stable if R 0 < 1.However, for the case where R 0 > 1, we proved that the disease is persistent.The theoretical results were confirmed by several numerical tests.
The paper is organized as follows.In Section 2, we describe a generalised compartmental model for Honeybee population dynamics when it is influenced by the seasonality.We prove that the virus-free periodic solution is stable if R 0 < 1 however the disease will persist if R 0 > 1.We give in Section 3 several numerical tests confirming the theoretical results.We finish by giving some concluding remarks in section 4.

Mathematical Modeling for Honeybee Population Dynamics
In what follows we present a mathematical model that combines the normal demographic dynamics of a honey bee colony with the dynamics of an infection affecting foraging bees outside the hive at first and then spreading to the rest of the colony.This mathematical model generalise the one given in [20] to a seasonal environment.We assume that adult bee population is divided into a number of hive bees H, and a number of foraging bees F. In the model to be described below, we extend this division into four categories, namely susceptible hive bees H s , infected hive bees H i , susceptible foraging bees F s , and infected foraging bees F i .The proposed model is governed by a system of four ordinary differential equations [20]: with positive initial condition (H i (0), F i (0), H s (0), F s (0)) ∈ R 4  + and the food equation with the positive initial condition f (0)

Preliminary. Consider a T-periodic m × m continuous matrix function denoted by A(t) that
it is irreducible and cooperative and consider the following equation admitting a fundamental matrix with positive entries as a solution.Let r(β A (T)) to be the spectral radius of the matrix β A (T).According to the Perron-Frobenius theorem, r(β A (T)) is the principal eigenvalue of β A (T).By using [21], we obtain Lemma 2.1.[21].(2.3) admits a positive T-periodic function x(t) such that w(t) = x(t)e at with a = 1 T ln(r(β A (T))).
Proposition 2.1.The compact set in is a positively invariant and attractor of trajectories of dynamics (2.1) with Proof.Using the dynamics (2.1), we obtain if In section 2.2, we aim to define the basic reproduction number; R 0 , the disease-free and then its global stability for R 0 ≤ 1.Later, in section 2.3, we aim to prove that compartments H i (t) and F i (t) 2.2.Disease-free trajectory.In this section, we shall define the expression of the basic reproduction number; R 0 , according to the definition given by the theory in [19].
The dynamics (2.7) admits a disease-free periodic trajectory Ȳ(t) = (0, 0, Hs (t), Fs (t)).Let where and Y i are the i-th components of f (t, Y(t)) and Y, respectively.An easy calculus gives us Therefore, the trajectory Ȳ(t) is linearly asymptotically stable in Therefore, the condition (A6) in [19, Section 1] also holds.
Let us define A + (t) and A − (t) to be two matrices defined by and where V(t, Y), respectively.An easy calculus gives us the expressions of matrices A + (t) and A − (t) as the following: Consider Z(s 1 , s 2 ) to be the two by two matrix solution of the system for any s 1 ≥ s 2 , with Z(s 1 , s 1 ) = I 2 , i.e., the 2 × 2 identity matrix.Therefore, condition (A7) is also fulfilled.
Denote by C T the ordered Banach space of T-periodic functions that are defined on R → R 2 , with the maximum norm .∞ and the positive cone C + T = {ϕ ∈ C T : ϕ(s) ≥ 0, for any s ∈ R}.Consider the linear operator L : C T → C T given by Therefore, the basic reproduction number, R 0 , of system (2.1) is given by R 0 = r(L).Thus, the local stability of the disease-free periodic trajectory, E 0 (t) = (0, 0, Hs (t), Fs (t)), of the dynamics (2.1) with respect to R 0 is given hereafter.
We deduce that E 0 (t) is asymptotically stable if R 0 < 1 and it is unstable if R 0 > 1.Now, we show that if R 0 < 1 then the disease-free periodic solution E 0 (t) = (0, 0, Hs (t), Fs (t)) is globally asymptotically stable and then the disease dies out.
Therefore, it remains to satisfy the global attractivity of E 0 (t) once R 0 < 1.Using (2.5) in Proposition 2.1, for any for t > T 1 .Let M 2 (t) be the two by two matrix function given hereafter Then, we deduce that the disease-free periodic trajectory E 0 (t) is globally attractive.Now, we show that if R 0 > 1 then H i (t) and F i (t) are uniform persistence and then the disease persists.

Numerical Examples
In this section, we adopted the numerical simulations validating analytical findings.For all numerical simulations, the periodic functions are given by with w | describe the seasonal cycles frequencies, however, φ describes the phase shift.The numerical values of Parameter Three scenarios were consider here.The first one was allocated to the case of fixed environment.
However, the second was allocated to the case where only the contact rates are seasonal.Finally, the last case were allocated to the case where all parameters are periodic.The numerical resolution was done using explicit Runge-Kutta formulas of orders 4 and 5 under Matlab.
3.1.Case of autonomous system.Let us start by the simple case where there is no influence of the seasonality on the dynamics.Thus, we restrict our attention on the autonomous dynamics (3.2), i.e., all parameters are positive constants.
with an initial condition (H i (0), The trivial steady state is given by We apply the next-generation matrix method introduced by Diekmann [23,24] to calculate the basic reproduction number for our system (3.2).See [25][26][27][28][29][30] for other applications.Let us define the matrices F and V by  with the inverse matrix V −1 of V given by , and the next generation matrix is defined by .
The characteristic polynomial is given by . The discriminant of the previous quadratic equation is given by Therefore the characteristic polynomial admits two roots the basic reproduction number for model (3.2) that it is defined as the spectral radius of the next-generation matrix, FV −1 is then given by.
However, in Figure 2, the calculated trajectories of the dynamics (3.2) converge to the disease-free steady state E 0 , then confirming the global asymptotic stability of E 0 if R 0 ≤ 1.
In Figure 1, the calculated trajectories of the dynamics (3.2) converge asymptotically to the periodic solution corresponding to the disease persistence if R 0 > 1.In Figure 2, different initial conditions were considered and for each one of them, the solution converge to the same periodic solution.In Figures 3 and 4, the calculated trajectories of the dynamics (3.2) converge to the disease-free steady    All the rest of parameters are fixed.We obtain the following system.
with the positive initial condition (H i (0), F i (0), H s (0), F s (0)) ∈ R 4 + .We give the results of some numerical simulations confirming the stability of the steady states of system (3.3).The approximation of the basic reproduction number R 0 was performed using the time-averaged system as in [14,15].Other definitions of R 0 can be found in [17,18].
In Figure 5, the calculated trajectories of the dynamics (3.3) converge asymptotically to the periodic solution corresponding to the disease persistence if R 0 > 1.In Figure 6, different initial conditions were considered and for each one of them, the solution converge to the same periodic solution.
In Figures 7 and 8, the calculated trajectories of the dynamics (3.3) converge to the disease-free periodic solution E 0 (t) = (0, 0, Hs (t), Fs (t)) for the case where R 0 ≤ 1.    3.3.Case of full periodic system.In the third step, we performed numerical simulations for the system (2.1)where all parameters were set as T-periodic functions.Thus the model is given by  with the positive initial condition (H i (0), F i (0), H s (0), F s (0)) ∈ R 4 + .We give the results of some numerical simulations confirming the stability of the steady states of system (3.4).The basic reproduction number R 0 was approximated by using the time-averaged system as in [14,15].Other definitions of R 0 can be found in [17,18].
In Figure 9, the calculated trajectories of the dynamics (3.4) converge asymptotically to the periodic solution corresponding to the disease persistence if R 0 > 1.In Figure 10, different initial conditions were considered and for each one of them, the solution converge to the same periodic solution.

Conclusions
In this paper, we consider the Honeybee population dynamics in a seasonal environment observed in real life.We defined the basic reproduction number, R 0 by using an integral operator.
It is proved that once R 0 ≤ 1, all solution of the dynamics converge to the disease-free periodic trajectory and that the disease persists if R 0 > 1.

Conflicts of Interest:
The authors declare that there are no conflicts of interest regarding the publication of this paper.
∈ R + .H s (t), H i (t), F s (t), F i (t) and f (t) describe the susceptible hive bees, infected hive bees, susceptible foraging bees, infected foraging bees, and the amount of food available at time t, respectively.L in (t), β HH (t), β HF (t), β FF (t), m w (t), d H (t), d F (t), R(t), γ A (t), γ L (t), and c(t) are continuous, positive T-periodic functions reflecting the influence of seasonality of the environment on the Honeybee population dynamics.β HH Contact rate between hive bees β HF Contact rate between hive bees and foraging bees.β FF Contact rate between foraging bees m w Natural death rate of bees during the winter season d H Death rate of hive bees due to infection d F Death rate of foraging bees due to infection R Recruitment rate of maturing hive bees to foraging duties c Foraging rate (gm/day) γ A Consumption rate of food by foragers and hive bees (gm/day) γ L Consumption rate of food by larvae (gm/day) m w L in The queen's egg laying rate per day Let ρ(t) to be a continuous, positive T-periodic function.Let us denote by ρ u = max t∈[0,T) ρ(t) and .15) P applied to the dynamics (2.13) admits a fixed point ( Hα s , Fα s ) that it is globally attractive with Hα s (t) > Hs − ξ, and Fα s (t) > Fs (t) − ξ; then, ∃ T 2 > 0 such that H s (t) > Hs (t) − ξ and F s (t) > Fs (t) − ξ for t > T 2 .Then, for t > T 2 , we have Hi , Fi , Hs , Fs ) ∈ Ω 0 with Ỹ0 ∈ Int(R + ) 2 × R 2 + .Suppose that Hs = 0. From the first equation of the dynamics (2.1), Hs (t) satisfies which contradicts the boundedness of the solution.Therefore, (2.14) is satisfied and P is weakly uniformly persistent with respect to (Ω 0 , ∂Ω 0 ).By applying Proposition 2.1, P has a global attractor.We deduce that Y 1 is an isolated invariant set inside Ω and that W s (Y 1 ) ∩ Ω 0 = ∅.All trajectories inside M ∂ converges to Y 1 which is acyclic in M ∂ .Applying [22, Theorem 1.3.1 and Remark 1.3.1],we deduce that P is uniformly persistent with respect to (Ω 0 , ∂Ω 0 ).Moreover, by using [22, Theorem 1.3.6],P has a fixed point Ỹ0 = ( = [β HH (t)H i (t) + β HF (t)F i (t)]H s (t) − (m w (t) + d H (t) + R(t))H i (t), Ḟi (t) = R(t)H i (t) + [β HF (t)H i (t) + β FF (t)F i (t)]F s (t) − (m w (t) + d F (t))F i (t), Ḣs (t) = m w (t)L in (t) − (m w (t) + R(t))H s (t) − [β HH (t)H i (t) + β HF (t)F i (t)]H s (t), Ḟs (t) = R(t)H s (t) − m w (t)F s (t) − [β HF (t)H i (t) + β FF (t)F i (t)]F s (t).