Impacts of seasonality and parasitism on honey bee population dynamics ^††thanks:

If there is no mites, i.e., $v(0)=0$, the model (2) reduces to the following bee-only population model with seasonality:

with $\bar{r}=r_{0}(1+\epsilon\cos(\frac{2\pi(t-\psi)}{\gamma}))$ which satisfies a Lipschitz condition for all $u\geq 0$. Thus according to Theorem 3.1, the initial value problem with $u(0)\geq 0$ has a unique non-negative and bounded solution.

In order to study the effects of the strength of seasonality ($\epsilon$) and the length of seasonality ($\gamma$) on bee populations, we start with the dynamics of the Honeybee-only model (4) when $\bar{r}(t)=r_{0}$ is a constant. The honeybee-only system without seasonality (4) has two equilibria $u_{i}^{*},i=1,2$ shown as below provided $r_{0}>2\bar{d}_{h}\sqrt{\hat{K}}$:

The global dynamics of (4) when $\bar{r}(t)=r_{0}$ can be summaries as the following proposition:

In Figure 1, without seasonality $\epsilon=0$, the honey bee colony with $r_{0}=1$, $\bar{d}_{h}=0.5$, $\hat{K}=1/4$, $\psi=0$ and $\gamma=100$ can survive under its initial condition $u(0)=1$ (red curve in Figure 1(a)) while it collapses under its initial condition $u(0)=0.1$ (red curve in Figure 1(b)). When the intensity of seasonality is not too high, i.e., $\epsilon=0.2$ or $0.5$, the honey bee colony can survive under its initial condition $u(0)=0.1$ (black and green curves in Figure 1(b)). This is an example showing that seasonality can promote the survival of a honey bee colony. On the other hand, When the intensity of seasonality is high, i.e., $\epsilon=0.8$ (blue curve in Figure 1(a)), the honey bee colony collapses with the initial condition of $u(0)=1$ when the honey bee colony can survive without seasonality. This is an example showing that seasonality can make honey bee colony collapse under certain conditions.

In order to explore the impact of the intensity of seasonality $\epsilon$, we first define the minimum and maximum value of the egg-laying rate function: $r_{m}=\min{\bar{r}(t)}=r_{0}(1-\epsilon)$ and $r_{M}=\max{\bar{r}(t)}=r_{0}(1+\epsilon)$. Motivated by Proposition 1, the intensity of seasonality can be classified into the following three cases:

1. 1.

   The low egg-laying rate if $r_{M}=r_{0}(1+\epsilon)\leq 2\bar{d}_{h}\sqrt{\hat{K}}$. This case is equivalent to

   $$0\leq\epsilon\leq 1-\frac{2\bar{d}_{h}\sqrt{\hat{K}}}{r_{0}}$$

2. 2.

   The high egg-laying rate if $r_{m}=r_{0}(1-\epsilon)\geq 2\bar{d}_{h}\sqrt{\hat{K}}$. This case is equivalent to

   $$0\leq\epsilon\leq\frac{2\bar{d}_{h}\sqrt{\hat{K}}}{r_{0}}-1$$

3. 3.

   The intermediate egg-laying rate if $r_{m}=\min{\bar{r}(t)}=r_{0}(1-\epsilon)<2\bar{d}_{h}\sqrt{\hat{K}}\leq r_{M}=\max{\bar{r}(t)}=r_{0}(1+\epsilon)$. This is the case when

   $$\max\{1-\frac{2\bar{d}_{h}\sqrt{\hat{K}}}{r_{0}},\frac{2\bar{d}_{h}\sqrt{\hat{K}}}{r_{0}}-1\}\leq\epsilon\leq 1.$$

Now we have the following theorem:

By applying Proposition 3.1 and the method used in Ratti et al.(2015) ratti2015mathematical , we obtain the stability condition when Model 4 processes a periodic solution $u^{*}$ as the following theorem:

To further address the impacts of seasonality on honey bee population dynamics, we provide basins of attractions for Model (4) in Figure 2 and Figure 3 by setting $\bar{d}_{h}=0.5,\hat{K}=1/4,r_{0}=1$. We set $\psi=0$ in Figure 2. The x-axis is the initial honey bee population $u(0)$, and the y-axis is the intensity of seasonality measured by $\epsilon$. Those parameter values gives $\epsilon_{c}=0.5$ which is a white horizontal line in Figure 2 and Figure 3. The blue region in Figures is the value of the strength of seasonality ($\epsilon$) and the corresponding initial conditions that lead the colony to collapse, while the red region is the value of $\epsilon$ and $u(0)$ that lead to the colony survival.

Figure 2 and Figure 3 suggest that the strength of seasonality ($\epsilon$), the length of seasonality ($\gamma$), and the time of the maximum laying rate ($\psi$) impact the survival of honey bee colony in the synergistic ways:

1. 1.

   The length of seasonality ($\gamma$) is small, e.g., $\gamma=4$:

   • •

     If the time of the maximum laying rate ($\psi$) is less than the half period $\gamma$, the seasonality seems to promote the survival of the colony in the sense that the initial bee population that originally leads to collapsing but it leads to colony survival with seasonality.

   • •

     If the time of the maximum laying rate ($\psi$) is larger than the half period $\gamma$, the seasonality seems to suppress the survival of the colony in the sense that the initial bee population that originally leads to survival but it leads to colony collapsing with seasonality.

2. 2.

   When the length of seasonality ($\gamma$) is larger, e.g., $\gamma=40,400$, the large intensity of seasonality $\epsilon$ can lead to the collapsing of the colony while the impacts of the smaller intensity of seasonality $\epsilon$ depends on the timing of the maximum laying rate ($\psi$) as follows:

   • •

     If the time of the maximum laying rate ($\psi$) is less than the half period $\gamma$, the seasonality seems to promote the survival of the colony.

   • •

     If the time of the maximum laying rate ($\psi$) is larger than the half period $\gamma$, the seasonality seems to suppress the survival of the colony.

In this subsection, we focus on dynamics of the honeybee-parasite interaction model (3) in the absence of seasonality, i.e., $\bar{r}(t)=\bar{r}$. Thus, we have the following rescaled model (5):

that would allow us to obtain biological insights on how parasitism impacts the honey bee population by comparing the dynamics of $v(0)=0$ versus $v(0)>0$. In the case that $v(0)=0$, the model (3) reduces to the honey bee only model in the constant environment (4) whose dynamics are summarized in Proposition 1.

Let $(u^{*},v^{*})$ be an equilibrium of Model (3), then it satisfies the following equations:

Solving Eqt.7 gives $v^{*}=0$ or $u^{*}=\frac{\bar{d}_{m}}{\omega-\bar{d}_{m}}$. And if $v^{*}=0$, then Eqt.6 is

which gives the following two positive solutions provided that $\bar{r}>2\bar{d}_{h}\sqrt{\bar{K}}$,

or

In the case that $\omega>\bar{d}_{m}$, we have $u^{*}=\frac{\bar{d}_{m}}{\omega-\bar{d}_{m}}$ and $v^{*}=\frac{\left[\bar{r}u^{*}-\bar{d}_{h}\left((u^{*})^{2}+\hat{K}\right)\right](1+u^{*})}{\omega((u^{*})^{2}+\hat{K})}$ as the unique interior equilibrium of Model 3. The stability of the equilibrium point can be evaluated through the following Jacobean matrix of Model 3 is

Now we are the following on the dynamics of the Honeybee-Parasite system (3):

Item (3) needs further theoretical exploration regarding how may parasitism destabilize the colony dynamics. For example, the colony destabilizes to show fluctuating dynamics through supercritical Hopf-bifurcation; or to collapse supercritical Hopf-bifurcation.

By applying the results in wang2011predator , our system 3 undergoes a Hopf-bifurcation. To study further, we re-scaled the system 3 to the following model:

where $g(u)=\frac{\omega u}{1+u}$ and $f(u)=\frac{\bar{r}}{g(u)}\cdot\frac{u^{2}}{\hat{K}+u^{2}}-\frac{\bar{d}_{h}}{g(u)}\cdot u$. We can verify that our system 8 satisfies the following conditions:

(a1) $f\in C^{1}(\bar{\mathbb{R}}),f(a)=f(b)=0$, where $0<a<b$; $f(u)$ is positive for $a<u<b$, and $f(u)$ is negative otherwise; there exists $\bar{\lambda}\in(a,b)$ such that $f^{\prime}(u)>0$ on $[a,\bar{\lambda})$, $f^{\prime}(u)<0$ on $(\bar{\lambda},b]$;

(a2) $g\in C^{1}(\bar{\mathbb{R}}),g(0)=0$; $g(u)>0$ for $u>0$ and $g^{\prime}(u)>0$ for $u>0$, and there exists $\lambda>0$ such that $g(\lambda)=d$.

(a3) $f(u)$ and $g(u)$ are $C^{3}$ near $\lambda=\bar{\lambda}$ and $f^{\prime\prime}(\bar{\lambda})<0$.

Then according to Theorem 3.1 in Wei et al. (2011) wang2011predator , we can conclude that our system 8 exists the first Lyapunov coefficient

where

Thus, we have the following results on Hopf-bifurcations:

We choose a typical example of our honeybee-parasite interaction model (3) by setting

which has a bistability between the colony collapsing state $(0,0)$ and the survival equilibrium at the locally stable equilibrium point $(5,5.5769)$ whose basins of attractions are red area shown in Figure 6(a).

To further explore the impacts of the seasonality strength $\epsilon$ and the period of seasonality $\gamma$ on the colony survival and population dynamics, without loss of generosity, we set the queen laying her maximum number of eggs at time $\psi=0$, and we perform the following simulations (Figure 6, 7, 8) on basin’s attractions of our honeybee-parasite model (3).

1. 1.

   When the period of the seasonality $\gamma$ is small, e.g., $\gamma=4$, comparisons of areas of basin attractions for the colony survival among Figure 6(a) (no seasonality), 6(c) (the seasonality strength $\epsilon=0.2$), and 6(d) (the seasonality strength $\epsilon=0.8$), suggest that seasonality strength $\epsilon$ may not impact the basin attractions of the colony survival but it impacts the population dynamics as shown in Figure 6(b). Simulations suggest that the larger value of the strength of seasonality $\epsilon$, the larger amplitude of the population.
   

   

   

   

   

2. 2.

   When the period of the seasonality $\gamma$ is in the intermediate range, e.g., $\gamma=80$, the impacts from the strength seasonality $\epsilon$ can be very complicated. For example, Figure 7(d) shows that basins of attractions for the colony survival are splitted into two red areas, and Figure 7(b) shows larger $\epsilon$ gives larger population amplitude.
   

   

   

   

   

3. 3.

   When the period of seasonality $\gamma$ is large, e.g., $\gamma=100,250$, comparisons of the basin attractions for the colony survival suggest that the small strength seasonality $\epsilon$ may not impact the basin attractions of the colony survival while its large value may cause the colony collapsing (see Figure 8(c) & 8(d)). In some cases, the large strength of seasonality $\epsilon$ may have a positive influence on the colony survival by increasing the area of basin attractions of the colony survival (see the comparison of Figure 8(a) & 8(f)). From the population dynamics point of view, Figure 8(b) suggests that the population has a larger amplitude when $\epsilon$ is larger.
   

   

   

   

   

   

   

Next, we explore the impacts of the timing of the maximum egg-laying rate ($\psi$) on colony survival and population dynamics in Figure 9 by fixing

$\gamma=70:$ 1) $\epsilon=0.2$ (the order from large to small): $\psi=10$ is largest, then $\psi=0$, $\psi=60=30$ these two cases have same survival area, $\psi=35$, and $\psi=40$ is the smallest.

2) $\epsilon=0.35$ (the order from large to small): $\psi=60$ is largest, then $\psi=0$, $\psi=40$, $\psi=10$, $\psi=35$, and $\psi=30$ is the smallest.

Notice that the seasonality period is $\gamma=70$ and $\epsilon=0.35$. We choose the timing of the maximum egg-laying rate $\psi\in\{0,10,30,35,40$ and $60\}$ and observe that the red area of the basin attractions for the colony survival is largest when the timing of the maximum laying rate ($\psi$) is $\psi=60$ (Figure 9(l)), then the second largest in the case when $\psi=0$ (Figure 9(g)), the smallest one is $\psi=30$ (Figure 9(i)), and the second smallest in the case when $\psi=35$ (Figure 9(j)). These observations from Figure 9 regarding the impacts of the maximum laying rate ($\psi$) of our honeybee-parasite model 3 seem to show similar trends of our honey bee-only model 4 (see Figure 3): as the $\psi$ increases, the seasonality can suppress the survival of the colony; and after the minimum survival area, the $\psi$ can promote the survival of the colony. But the significant difference with the bee-only model is the smallest area is not $\psi=\frac{\gamma}{2}$. Figure 9(m) and 9(n) show how different timing of the maximum egg-laying rate $\psi$ can lead to different colony dynamics.

To further understand the impacts of the timing of the maximum laying rate ($\psi$) of our honeybee-parasite model 3, we set the strength of the seasonality being $\epsilon=0.2$, and choose the timing of the maximum egg-laying rate $\psi\in\{0,10,30,35,40$ and $60\}$, respectively. The basin attractions for the colony survival is largest when the timing of the maximum laying rate ($\psi$) is $\psi=10$(Figure 9(b)), the second largest in the case when $\psi=0$ (Figure 9(a)), the smallest one is $\psi=40$ (Figure 9(e)), and the second smallest in the case when $\psi=35$ (Figure 9(d)). These observations are different than the case of $\epsilon=0.35$ shown in Figure 9 and the case of the honey bee only model 4 (see Figure 3). The significant difference is that $\psi$ can promote the survival of the colony at the very beginning of $\psi$ growth ($\psi=10$ in our simulation). These comparisons and our further simulations suggest that the impacts of the timing of the maximum laying rate ($\psi$) on the honey bee colony survival in the presence of parasitism are very complicated. The area of the basin attractions for the colony survival may be increasing or decreasing with respect to the value of $\psi$ and $\epsilon$ without clear patterns.

By comparing the basins of attractions of the honeybee-mite system without seasonality in Figure 7(a) to the honeybee-mite system with seasonality in Figure 7(d), 9(e) & 9(g), we observe that seasonality can split the basins of attractions into disconnected regions. This may lead to two scenarios after adding seasonality: (1) the colony may survive from collapsing (see Point A in Figure 10(c) versus Figure 10(d)), and (2) the colony may be prone to collapsing (see Point B in Figure 10(c) versus Figure 10(d)). This suggests that seasonality may generate varied outcomes depending on initial conditions. For instance, while an initial rise in the parasite population is generally perceived as detrimental, it can enhance colony survival under specific circumstances, particularly when considering seasonal factors (compare points A and B in Figure 10(d)). This phenomenon has been observed in experimental data degrandi2020can . To illustrate those observations, we use Figure 7(d) as an example where we list four cases. Among them, the initial bees population of Colony 1 (Case 1, blue) and Colony 3 (Case 3, gray) are similar, and Colony 2 (Case 2, red) and Colony 4 (Case 4, black) are close. While, the initial mite population is increasing in the order of Case 1, Case 2, Case 3, and Case 4. Figure 10(a) shows the colony of Case 1 and Case 3 survived while Case 2 and Case 4 collapsed, especially Colony 3 has fewer bees and more mites than Colony 2, but survives.

The potential biological explanation for this phenomenon lies in the heart of seasonality impacts on the egg laying rate incorporated in the model, and mite population is impacted through bee population. For the mite population to grow, colonies must have enough bees. If the system doesn’t consider the impacts of seasonality (Model 5), then a higher parasite level ($v_{0}$: Point A $>$ Point B) leads to colony collapse (see red curves in Figure 11(b) & 11(d)), because parasitism reduces colony population growth by shortening the lifespan of adult workers, then the population of bees is reduced and so will Varroa population growth. However, the system with seasonality (Model 3) leads to a switch in the outcomes of these two colonies, which is survival colony goes to collapse because of seasonality, whereas the collapsing colony becomes survival. The point is that the egg-laying rate of bees is periodic due to seasonal effects, then the number of bees will increase at some time intervals (the green curve in Figure 11(a)). At a higher parasite level, fewer bees will bring the mite population down ($\frac{\omega u}{1+u}v$) to a manageable level, and the seasonality egg-laying rate helps the colony grow up periodically (seasonality in Point A). At a lower parasite level, seasonality also leads mites to grow up more than without seasonality effects (Figure 11(d)). Seasonality and high numbers of bees may lead to excessive mite growth beyond the colony’s sustainable threshold and colony collapse. This principle is similar to one method of controlling Varroa mites: removing the brood from the hive and interrupting the brood reproductive cycle. With no brood present, mites are compelled to feed on adult bees, which can limit the mites’ ability to reproduce, helping to control their populations jack2021integrated . Nevertheless, this method will be affected by seasonality. Removing lots of broods in the fall may have strong negative impacts on overwintering survival jack2020evaluating .

Now we explore the impacts of parasitism $\omega$ on honey bee population dynamics and its colony’s survival in Figure 12. Comparison of black areas (which is the basins of attractions of only honey bee survival) in Figure 12(d) & 12(a) suggest that small parasitism (e.g., $\omega=0.18$) with seasonality is more likely to lead to the colony survival than the case without seasonality. When parasitism is not small (e.g., $\omega=0.30$) (see Figure 12(b)), seasonality can destabilize the system and decrease the average population of the honey bee. When $\omega$ is large (e.g., $\omega=0.5$), parasitism has negative impacts on the honey bee colony that lead the colony to collapse (see all blue areas in Figure 12(c)). Figure 12(e) shows increasing parasitism, colonies may still survive but the average population of honey bees decreases (see black and green curves). These observations are in line with our theorem 3.4 for the case without seasonality.

We choose a stable limit cycle example of our honeybee-parasite interaction model 3 by setting

which has a stable collapsing state $(0,0)$ for the colony, and a stable limit cycle around the source interior equilibrium $(2.33,0.3875)$ whose basins of attractions are red area shown in Figure 13(e).

We explore the impacts of the seasonality strength $\epsilon$, the period of seasonality $\gamma$, the queen laying her maximum number of eggs at time $\psi$, and the parasitism effects $\omega$ on the colony survival and population dynamics. We perform the following simulation in Figure 13 on the basin’s attractions of our honeybee-parasite model 3. We set the queen laying her maximum number of eggs at time $\psi=0$ to observe the impacts of $\gamma$ and $\epsilon$.

1. 1.

   When the period of the seasonality is small, i.e., $\gamma=4$, comparisons of areas of basin attractions for the colony survival among Figure 13(e) (no seasonality), 13(a) (the seasonality strength $\epsilon=0.2$) suggest that small seasonality strength $\epsilon$ may not significantly impact the survival of the colony much but larger seasonality strength $\epsilon$ can generate larger population amplitude (see Figure 13(f)).

2. 2.

   When the period of the seasonality is in the intermediate range, e.g., $\gamma=40$, comparisons of areas of basin attractions for the colony survival among Figure 13(e) (no seasonality), 13(c) (the seasonality strength $\epsilon=0.2$) and 13(d) (the seasonality strength $\epsilon=0.5$) suggest that seasonality strength $\epsilon$ seems to suppress the survival of the colony.

3. 3.

   When the the seasonality strength $\epsilon$ is fixed, increasing the period of the seasonality $\gamma$ seems to suppress the survival of the colony (See Figures 13(a) & 13(c) and Figures 13(b) & 13(d)).

4. 4.

   The large $\gamma$ and $\epsilon$ would lead to the colony collapsing as we observe that the colony collapses when $\gamma>60$.

Let the period of the seasonality be $\gamma=40$ and the strength of the seasonality be $\epsilon=0.2$. We explore the impacts of the timing of the maximum egg-laying rate $\psi$ by varying $\psi$=0, $15(<\frac{\gamma}{2}=20)$, $35(>\frac{\gamma}{2}=20)$ in Figure 14. We observe that the basin attractions for the colony survival seem to have similar shapes: the largest area is $\psi=0$ (Figure 14(a)), the second largest being $\psi=35$ (Figure 14(c)), and the smallest one is $\psi=15$ (Figure 14(b)). The observation under this particular parameter set regarding the impacts of $\psi$ of our honeybee-parasite model 3 seems to show similar trends as our honey bee only model 4 (see Figure3). Figure 14(d) provides some visual insights on how may $\gamma$ and $\epsilon$ impact population dynamics.

Let the period of the seasonality be $\gamma=40$ and the strength of the seasonality be $\epsilon=0.2$. We explore the impacts of the parasitism by varying $\omega\in[0.1,0.35]$ in Figure 15. We observe follows:

1. 1.

   When the parasitism $\omega$ is small (e.g., $\omega=0.1$), the honey bee can survive while the parasite dies out (see Figure 15(a)).

2. 2.

   When parasitism is increased to $\omega=0.292$, colonies can survive with parasitism, but the area of basins of attractions for survival decreases as parasitism increases (see Figure 15(a), 15(b) & 15(c)). Thus parasitism has a negative influence on the colonies’ survival.

3. 3.

   When the parasitism is large (e.g., $\omega>0.3$ when $u_{0}=5,v_{0}=0.04$), colonies collapse.

We observe that (1) If the colony can survive, increasing the parasitism attack degree can decrease the average population of honey bees. (2) Large parasitism can lead to a colony collapsing. In general, Seasonality with parasitism can have negative impacts in terms of either decreasing the average honey bee population or the colony collapsing.

Without seasonality, the value of parasitism rate $\omega$ can lead to destability through hopf-bifurcation (see Theorem 3.5). To further explore how parasitism may impact the honeybee population dynamics with or without seasonality, we perform bifurcation on the impacts of parasitism $\omega=[0.2,0.33]$ with (see Figure 15(e)) or without (see Figure 15(f)) seasonality by setting $\bar{r}_{0}=1$, $\bar{d}_{h}=0.2$, $\bar{d}_{m}=0.21$, $\psi=0$, $\gamma=40$, $\epsilon=0.2$, $\hat{K}=4.49$, $u(0)=5$ and $v(0)=0.04$.

In the absence of seasonality (Model 3), the $\omega_{3}$ is the bifurcation value where the mite-free equilibrium ($(\bar{N}^{*}_{h},0)$) changes from being locally stable to unstable (see Theorem 3.4 item (2)), and the coexistence of bee and mite population emerges as the locally stable interior equilibrium, and the interior equilibrium become unstable (see Theorem 3.4 item (3)) through supercritical Hopf-bifurcation at $\omega_{4}$, where exists the stable limit cycle (see Theorem 3.5). After the value of $\omega_{5}$, the colony collapses. Therefore, the bifurcation diagram in Figure 15(f)) suggests that: (1) when the severity of parasitism ($\omega$) is small, the colony survives with non-parasites; (2) when the value of $\omega$ rise, bees and parasites coexist in the colony and gradually decreases the population of bees; (3) under the conditions of supercritical Hopf-bifurcation, bees and parasites coexist in a periodic state; (4) when $\omega$ is large enough, the parasites leads to the colony collapse.

In the seasonality model (Model 4 and see Figure 15(e)), before the value of $\omega_{1}$, the system is locally stable around mite-free solutions (see Theorem 4.1); after this bifurcation point, bees and parasites coexist as the periodic interior solutions. The $\omega_{2}$ is the critical value when colony collapses.

We observe that seasonality can delay the impact of parasitism in two bifurcation points: (1) $\omega_{1}>\omega_{3}$: parasite needs larger attacking rates to survive in the periodic environment. And (2) $\omega_{2}>\omega_{5}>\omega_{4}$: colony can still survive with the larger attacking rates from parasites in the periodic environment.