Mathematical Analysis for Honeybee Dynamics Under the Influence of Seasonality

Abstract

In this paper, we studied a mathematical model for honeybee population diseases under the influence of seasonal environments on the long-term dynamics of the disease. The model describes the dynamics of two different beehives sharing a common space. We computed the basic reproduction number of the system as the spectral radius of either the next generation matrix for the autonomous system or as the spectral radius of a linear integral operator for the non-autonomous system, and we deduced that if the reproduction number is less than unity, then the disease dies out in the honeybee population. However, if the basic reproduction number is greater than unity, then the disease persists. Finally, we provide several numerical tests that confirm the theoretical findings.

1. Introduction

The term “nosemosis” appeared in 1914, and comes from the translation of the German term “nosema-seuche” invented at the time by Enoch Zander. It is one of the first bee diseases described, with its etiological agent Nosema apis having been identified in 1909. It was then widely studied by European and North American scientists, and quickly considered “one of the most interesting diseases of adult bees” (Toumanoff, 1930). Today, nosemosis due to N. apis is an integral part of universal beekeeping medical culture, in the same way as foulbrood (American and European) and varroosis. However, the latest epidemiological data available to us seem to indicate a gradual disappearance of this disease in France. This biological phenomenon, a priori fortunate, hides complex host–parasite relationships and evolutionary mechanisms. Indeed, the recent discovery of another species of Nosema (N. ceranae) parasites on the European honeybee (Apis mellifera) has highlighted a competition that was previously ignored between the native species and this new emerging species for the same ecological niche, the bee intestine. This competition now seems to largely benefit the new parasite, thus is responsible for the near disappearance of the native parasite and the disease associated with it. Nosemosis is a parasitic disease that affects the three castes of adult bees (queen, drone, and worker bee), due to the proliferation in the epithelium of the midgut of bees of a fungus of the genus Nosema (phylum Microsporidia, kingdom Fungi), Nosema apis. Its parasitic cycle includes a vegetative phase (in the ventricle of bees) and a phase of resistance and dissemination (in the external environment, by a spore). The complete cycle lasts between two and five days. The spore, with an average size of 5–6 µm in length and 2–3 µm in width, has a thick wall (except on the anchoring disk) and an evaginable polar tube of 100 to 400 µm in length. This wall effectively and durably protects the spore. Bees become contaminated by ingesting these spores during inter-individual exchanges (trophallaxis) and during their cleaning (grooming). Once ingested, the spore finds, in the bee’s ventricle, the conditions necessary for its germination (favorable chemical environment), which consists of the evagination of its polar tube, through which the sporoplasm (cytoplasm, nucleus and organelles of the spore) will be transmitted to the host epithelial cell by perforation. The multiple intracellular duplications of the parasite will ultimately lead to the destruction of the host cell, thus releasing numerous spores into the intestinal lumen, which will themselves be able to renew this cycle (auto-infection). The clinical consequences of this infection are, among others and for the most obvious, diarrhea, trailing bees, and major depopulations. Targeting the adult population of the colony, nosemosis is a protozoan (microorganism) disease that attacks the lining of the intestinal tract, weakening the hive by reducing the lifespan of sick adult bees. This disease can be mild, reducing a colony’s production by nearly half, or more severe, killing an entire colony. Since bees depend on the resources available in nature, the colony varies according to the seasons. In winter, with no flowers to forage while the cold is severe, the swarm lives grouped together in the hive. In spring, when pollen and nectar are abundant, the queen lays eggs; the colony grows and produces honey in quantity. There are many reasons why we may need to insert two colonies together in the same natural space. It may be that, during the season, our colony numbers have increased due to taking in swarms or performing artificial splits/swarms for swarm prevention or control.

To understand how to control and eradicate infectious diseases in the honeybee population is one of the main goals of mathematical epidemiology. Several sophisticated mathematical models predicting the honeybee population dynamics have been proposed [1,2,3,4]. Seasonal changes can considerably affect the food amount for both colonies, which is why taking account of the seasonality in a mathematical model predicting the behavior of the honeybees dynamics become a necessity. Several previous mathematical models describing other infectious diseases considering the seasonality [5,6,7,8,9,10,11,12] are studied. In particular, in [13], the authors studied a dynamical system with model parameters as time-periodic functions describing the dynamics of a honeybee colony that is infected with Varroa mites. In [14], the authors studied a non-autonomous nonlinear honeybee–parasite dynamical model considering the seasonality’s relation to the queen’s egg-laying rate. In [15], the authors considered a non-autonomous dynamics describing the spread of a disease within a bee colony in a seasonal environment by considering time-dependent parameters reflecting the periodicity of weather.

In this paper, we assume that we are bringing two colonies of bees from two different hives together in the same natural space, and we assume that both colonies are target of the nosemosis. We proposed a mathematical model predicting the behavior of both colonies asymptotically. In Section 2, we present the mathematical model studied in this paper. In Section 3, we studied the case of fixed environment by presenting the methods and main results of the study. Precisely, we will analyze the dynamical behavior of the model. In Section 4, we studied the case of seasonal environment where we discussed the methods and main findings. In particular, we will present the periodic dynamical behavior of the model. In Section 5, the analytical findings in both Section 3 and Section 4 will be supported by some numerical examples simulated under MATLAB R2024a software. Finally, in Section 6, we provide a brief conclusion.

2. Mathematical Model

Let us consider two beehives (A and B) sharing the same natural space, where each colony is subdivided into four classes as follows:

• Healthy bees, which are classified as either healthy bees from the first beehive ($A_s$) or healthy bees from the second beehive ($B_s$).

• The first beehive ($A_l$) and the second beehive ($B_l$), which are latently infected (latency stage).

• The first beehive ($A_i$) and the second beehive ($B_i$), which are already infected with the Nosema disease.

• The first beehive ($A_r$) and the second beehive ($B_r$), who recovered from the Nosema disease.

• We assume that the infected bees from the first beehive ($A_i$) and the infected bees from the second beehive ($B_i$) do not contribute equally in the transmission of the Nosema disease. The Nosema disease is caused by intracellular microsporidia parasites in the genus Nosema [16,17]. Therefore, we added a compartment reflecting the contribution of each infected bees’ category in the Nosema disease transmission in the environment, denoted by $I_c$.

Let us denote the honeybee population at time t in the first and second beehives by $A\left(t\right)$ and $B\left(t\right)$, such that $C\left(t\right)=A\left(t\right)+B\left(t\right)$ is the total number of honeybees in the considered space. We subdivided each population of bees into compartments (Figure 1) as follows: susceptible bees from the first beehive $A_s\left(t\right)$, susceptible foraging bees from the second beehive $B_s\left(t\right)$, latently infected bees from the first beehive $A_l\left(t\right)$, latently infected bees from the second beehive $B_l\left(t\right)$, infectious bees from the first beehive $A_i\left(t\right)$, and infectious bees from the second beehive $B_i\left(t\right)$. We denote the compartment reflecting the impact of each infected bees’ category in the Nosema disease transmission in the colony by $I_c\left(t\right)$.

We aim to give an explanation, as simply as possible, of the honeybees dynamics in a colony. The mathematical model that we proposed here is governed by the following ninth-dimensional dynamical system of ordinary differential equations. In reality, it consists of merging two “SEIR” epidemic models with a common vector virus-carrying [18].

(1)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{A}}_s\left(t\right)=m_a\left(t\right){\mathsf{\Lambda }}_a\left(t\right)-{\beta }_a\left(t\right)A_s\left(t\right)I_c\left(t\right)-m_a\left(t\right)A_s\left(t\right),\hfill \\ {\dot{A}}_l\left(t\right)={\beta }_a\left(t\right)A_s\left(t\right)I_c\left(t\right)-({\zeta }_a\left(t\right)+m_a\left(t\right))A_l\left(t\right),\hfill \\ {\dot{A}}_i\left(t\right)={\zeta }_a\left(t\right)A_l\left(t\right)-({\delta }_a\left(t\right)+m_a\left(t\right))A_i\left(t\right),\hfill \\ {\dot{A}}_r\left(t\right)={\delta }_a\left(t\right)A_i\left(t\right)-m_a\left(t\right)A_r\left(t\right),\hfill \\ {\dot{B}}_s\left(t\right)=m_b\left(t\right){\mathsf{\Lambda }}_b\left(t\right)-{\beta }_b\left(t\right)B_s\left(t\right)I_c\left(t\right)-m_b\left(t\right)B_s\left(t\right),\hfill \\ {\dot{B}}_l\left(t\right)={\beta }_b\left(t\right)B_s\left(t\right)I_c\left(t\right)-({\zeta }_b\left(t\right)+m_b\left(t\right))B_l\left(t\right),\hfill \\ {\dot{B}}_i\left(t\right)={\zeta }_b\left(t\right)B_l\left(t\right)-({\delta }_b\left(t\right)+m_b\left(t\right))B_i\left(t\right),\hfill \\ {\dot{B}}_r\left(t\right)={\delta }_b\left(t\right)B_i\left(t\right)-m_b\left(t\right)B_r\left(t\right),\hfill \\ \dot{I_c}\left(t\right)={\eta }_a\left(t\right)A_i\left(t\right)+{\eta }_b\left(t\right)B_i\left(t\right)-m\left(t\right)I_c\left(t\right),\hfill \end{array}\right.\end{array}$

with initial condition $(A_s^0,A_l^0,A_i^0,A_r^0,B_s^0,B_l^0,B_i^0,B_r^0,I_c^0)\in {\mathbb{R}}_{+}^9$.

All model parameters are assumed to be continuous, bounded, and T-periodic positive functions. ${\mathsf{\Lambda }}_a$ (respectively, ${\mathsf{\Lambda }}_b$) denotes the recruitment rate of bees through birth and movement of susceptible bees for the first beehive (respectively, the second beehive), $m_a$ (respectively, $m_b$) is the natural mortality rate of the bees from the first beehive (respectively, the second beehive). ${\displaystyle \frac{1}{m_a}}$ and ${\displaystyle \frac{1}{m_b}}$ are the average lifespans of bees in the first and second beehives, respectively. ${\beta }_a$ (respectively, ${\beta }_b$) describes the indirect Nosema disease transmission rate from the first beehive (respectively, the second beehive), ${\zeta }_a$ (respectively, ${\zeta }_b$) is the incubation rate of the bees from the first beehive (respectively, the second beehive), and then ${\displaystyle \frac{1}{{\zeta }_a}}$ and ${\displaystyle \frac{1}{{\zeta }_b}}$ are the incubation periods, ${\delta }_a$ (respectively, ${\delta }_b$) is the recovery rate of the bees from the first beehive (respectively, the second beehive). ${\displaystyle \frac{1}{{\delta }_a}}$ (respectively, ${\displaystyle \frac{1}{{\delta }_b}}$) is the average infectious period of the bees from the first beehive (respectively, the second beehive). Infected bees contribute to the Nosema disease transmission in the environment at rates ${\eta }_a$ and ${\eta }_b$, and have a natural death rate m. A summary of the model parameters’ signification is given in Table 1.

Since the fourth and eighth equations of the dynamics of (1) do not affect the other equations, we can restrict the dynamics to ${\mathbb{R}}_{+}^7$ instead of ${\mathbb{R}}_{+}^9$ as follows:

(2)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{A}}_s\left(t\right)=m_a\left(t\right){\mathsf{\Lambda }}_a\left(t\right)-{\beta }_a\left(t\right)A_s\left(t\right)I_c\left(t\right)-m_a\left(t\right)A_s\left(t\right),\hfill \\ {\dot{A}}_l\left(t\right)={\beta }_a\left(t\right)A_s\left(t\right)I_c\left(t\right)-({\zeta }_a\left(t\right)+m_a\left(t\right))A_l\left(t\right),\hfill \\ {\dot{A}}_i\left(t\right)={\zeta }_a\left(t\right)A_l\left(t\right)-({\delta }_a\left(t\right)+m_a\left(t\right))A_i\left(t\right),\hfill \\ {\dot{B}}_s\left(t\right)=m_b\left(t\right){\mathsf{\Lambda }}_b\left(t\right)-{\beta }_b\left(t\right)B_s\left(t\right)I_c\left(t\right)-m_b\left(t\right)B_s\left(t\right),\hfill \\ {\dot{B}}_l\left(t\right)={\beta }_b\left(t\right)B_s\left(t\right)I_c\left(t\right)-({\zeta }_b\left(t\right)+m_b\left(t\right))B_l\left(t\right),\hfill \\ {\dot{B}}_i\left(t\right)={\zeta }_b\left(t\right)B_l\left(t\right)-({\delta }_b\left(t\right)+m_b\left(t\right))B_i\left(t\right),\hfill \\ \dot{I_c}\left(t\right)={\eta }_a\left(t\right)A_i\left(t\right)+{\eta }_b\left(t\right)B_i\left(t\right)-m\left(t\right)I_c\left(t\right),\hfill \end{array}\right.\end{array}$

with initial condition $(A_s^0,A_l^0,A_i^0,B_s^0,B_l^0,B_i^0,I_c^0)\in {\mathbb{R}}_{+}^7$.

3. The Case of a Fixed Environment

We start by studying the autonomous dynamics where the model parameters are constants.

(3)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{A}}_s=m_a{\mathsf{\Lambda }}_a-{\beta }_a{A}_s{I}_c-m_a{A}_s,\hfill \\ {\dot{A}}_l={\beta }_a{A}_s{I}_c-({\zeta }_a+m_a)A_l,\hfill \\ {\dot{A}}_i={\zeta }_a{A}_l-({\delta }_a+m_a)A_i,\hfill \\ {\dot{B}}_s=m_b{\mathsf{\Lambda }}_b-{\beta }_b{B}_s{I}_c-m_b{B}_s,\hfill \\ {\dot{B}}_l={\beta }_b{B}_s{I}_c-({\zeta }_b+m_b)B_l,\hfill \\ {\dot{B}}_i={\zeta }_b{B}_l-({\delta }_b+m_b)B_i,\hfill \\ \dot{I_c}={\eta }_a{A}_i+{\eta }_b{B}_i-m{I}_c,\hfill \end{array}\right.\end{array}$

with initial condition $(A_s^0,A_l^0,A_i^0,B_s^0,B_l^0,B_i^0,I_c^0)\in {\mathbb{R}}_{+}^7$. Since all parameters are non-negative, we should prove that the trajectory of the dynamics of (3) is non-negative. The biologically feasible region of the dynamics of (3), denoted by $\mathsf{\Sigma }$, is defined as follows:

(4)$\begin{array}{c}\hfill \begin{array}{c}\mathsf{\Sigma }=\{(A_s,A_l,A_i,B_s,B_l,B_i,I_c)\in {\mathbb{R}}_{+}^7:A_s+A_l+A_i\le {\mathsf{\Lambda }}_a,\hfill \\ \hfill B_s+B_l+B_i\le {\mathsf{\Lambda }}_b,I_c\le {\displaystyle \frac{{\eta }_a{\mathsf{\Lambda }}_a+{\eta }_b{\mathsf{\Lambda }}_b}{m}}\}.\end{array}\hfill \end{array}$

Note that the set $\mathsf{\Sigma }$ is positively attracting of the solutions of (3) according to the following Lemma.

Lemma 1. 

Σ is a positively invariant, attracting, and compact set for the dynamics of (3).

Proof. 

Since we have ${\dot{A}}_s{\mid }_{A_s=0}=m_a{\mathsf{\Lambda }}_a>0$, ${\dot{A}}_l{\mid }_{A_l=0}={\beta }_a{A}_s{I}_c>0$, ${\dot{A}}_i{\mid }_{A_i=0}={\zeta }_a{A}_l>0$, ${\dot{B}}_s{\mid }_{B_s=0}=m_b{\mathsf{\Lambda }}_b>0$, ${\dot{B}}_l{\mid }_{B_l=0}={\beta }_b{B}_s{I}_c>0$, ${\dot{B}}_i{\mid }_{B_i=0}={\zeta }_b{B}_l>0$, $\dot{I_c}{\mid }_{I_c=0}={\eta }_a{A}_i+{\eta }_b{B}_i>0$. Therefore, ${\mathbb{R}}_{+}^7$ is invariant by the dynamics of (3). Let us denote by $A:=A_s+A_l+A_i$ and $B:=B_s+B_l+B_i$ the sizes of the total bees’ compartments from the first and second beehives, respectively. From the dynamics of (3), we have $\dot{A}=m_a{\mathsf{\Lambda }}_a-m_a{A}_s-m_a{A}_l-({\delta }_a+m_a)A_i\le m_a({\mathsf{\Lambda }}_a-A)$. Hence, $A\le {\mathsf{\Lambda }}_a$ if $A\left(0\right)\le {\mathsf{\Lambda }}_a$. Similarly, $\dot{B}=m_b{\mathsf{\Lambda }}_b-m_b{B}_s-m_b{B}_l-({\delta }_b+m_b)B_i\le m_b({\mathsf{\Lambda }}_b-B)$. Hence, $B\le {\mathsf{\Lambda }}_b$ if $B\left(0\right)\le {\mathsf{\Lambda }}_b$. Similarly, ${\dot{I}}_c={\eta }_a{A}_i+{\eta }_b{B}_i-m{I}_c\le {\eta }_a{\mathsf{\Lambda }}_a+{\eta }_b{\mathsf{\Lambda }}_b-m{I}_c$. Hence, $I_c\le {\displaystyle \frac{{\eta }_a{\mathsf{\Lambda }}_a+{\eta }_b{\mathsf{\Lambda }}_b}{m}}$ if $I_c\left(0\right)\le {\displaystyle \frac{{\eta }_a{\mathsf{\Lambda }}_a+{\eta }_b{\mathsf{\Lambda }}_b}{m}}$. □

Now, we need to compute the basic reproduction number of (3), denoted by ${\mathcal{R}}_0$, by using the next generation matrix [19]. We need first to determine the Nosema disease-free steady state of (3) having the following form:

(5)$\begin{array}{c}\hfill \begin{array}{c}{\mathcal{E}}_0=\left({\mathsf{\Lambda }}_a,0,0,{\mathsf{\Lambda }}_b,0,0,0\right).\hfill \end{array}\end{array}$

In our case, we have five equations that describe the production of new infections and changes in state among infected individuals (Equations (2), (3), (5), (6), and (7) of System (3)). In order to construct the next generation matrix, we define the matrix F as representing the rate of appearance of new infections in these five equations, and V as representing the rate of transfer of individuals into and out of these compartments by all other means. The non-negative matrix F and the non-singular matrix V are then given by

$\begin{array}{c}\hfill \begin{array}{c}F=\left[\begin{array}{ccccc}0& 0& 0& 0& {\beta }_a{\mathsf{\Lambda }}_a\\ 0& 0& 0& 0& {\beta }_b{\mathsf{\Lambda }}_b\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right] \mathrm{and} V=\left[\begin{array}{ccccc}{\zeta }_a+m_a& 0& 0& 0& 0\\ 0& {\zeta }_b+m_b& 0& 0& 0\\ -{\zeta }_a& 0& {\delta }_a+m_a& 0& 0\\ 0& -{\zeta }_b& 0& {\delta }_b+m_b& 0\\ 0& 0& -{\eta }_a& -{\eta }_b& m\end{array}\right].\hfill \end{array}\end{array}$

The spectral radius of (3) is the spectral radius of the next generation matrix $F{V}^{-1}$, expressed as follows:

(6)$\begin{array}{c}\hfill \begin{array}{c}{\mathcal{R}}_0={\displaystyle \frac{{\beta }_a{\zeta }_a{\eta }_a{\mathsf{\Lambda }}_a}{m({\zeta }_a+m_a)({\delta }_a+m_a)}}+{\displaystyle \frac{{\beta }_b{\zeta }_b{\eta }_b{\mathsf{\Lambda }}_b}{m({\zeta }_b+m_b)({\delta }_b+m_b)}}={\mathcal{R}}_1+{\mathcal{R}}_2.\hfill \end{array}\end{array}$

Theorem 1.

If ${\mathcal{R}}_0>1$, then (3) admits a unique endemic steady state.

Proof. 

Let the endemic steady state of (3) be denoted by ${\mathcal{E}}^{\ast }=(A_s^{\ast },B_s^{\ast },A_l^{\ast },B_l^{\ast },A_i^{\ast },B_i^{\ast },I_c^{\ast })$. This steady state is determined by solving the following system of equations:

(7)$\begin{array}{c}\hfill \left\{\begin{array}{c}{m}_a{\mathsf{\Lambda }}_a-{\beta }_a{A}_s^{\ast }{I}_c^{\ast }-m_a{A}_s^{\ast }=0,\hfill \\ {\beta }_a{A}_s^{\ast }{I}_c^{\ast }-({\zeta }_a+m_a)A_l^{\ast }=0,\hfill \\ {\zeta }_a{A}_l^{\ast }-({\delta }_a+m_a)A_i^{\ast }=0.\hfill \end{array}\right.\end{array}$

(8)$\begin{array}{c}\hfill \left\{\begin{array}{c}{m}_b{\mathsf{\Lambda }}_b-{\beta }_b{B}_s^{\ast }{I}_c^{\ast }-m_b{B}_s^{\ast }=0,\hfill \\ {\beta }_b{B}_s^{\ast }{I}_c^{\ast }-({\zeta }_b+m_b)B_l^{\ast }=0,\hfill \\ {\zeta }_b{B}_l^{\ast }-({\delta }_b+m_b)B_i^{\ast }=0.\hfill \end{array}\right.\end{array}$

and

(9)$\begin{array}{c}\hfill \begin{array}{c}{\eta }_a{A}_i^{\ast }+{\eta }_b{B}_i^{\ast }-m{I}_c^{\ast }=0.\hfill \end{array}\end{array}$

We obtain

(10)$\begin{array}{c}\hfill \begin{array}{c}{A}_l^{\ast }={\displaystyle \frac{({\delta }_a+m_a)}{{\zeta }_a}}{A}_i^{\ast },B_l^{\ast }={\displaystyle \frac{({\delta }_b+m_b)}{{\zeta }_b}}{B}_i^{\ast },I_c^{\ast }={\displaystyle \frac{{\eta }_a}{m}}{A}_i^{\ast }+{\displaystyle \frac{{\eta }_b}{m}}{B}_i^{\ast }.\hfill \end{array}\end{array}$

Substituting (10) into the second equations of (7) and (8) leads to

(11)$\begin{array}{c}\hfill \begin{array}{c}\left({\displaystyle \frac{{\beta }_a{\eta }_a{A}_s^{\ast }}{m}}-{\displaystyle \frac{({\delta }_a+m_a)({\zeta }_a+m_a)}{{\zeta }_a}}\right)A_i^{\ast }+\left({\displaystyle \frac{{\beta }_b{\eta }_b{B}_s^{\ast }}{m}}-{\displaystyle \frac{({\delta }_b+m_b)({\zeta }_b+m_b)}{{\zeta }_b}}\right)B_i^{\ast }=0.\hfill \end{array}\end{array}$

We can write

(12)$$\begin{gathered} \begin{array}{}\hfill \begin{array}{cc}& {\displaystyle \frac{({\delta }_a+m_a)({\zeta }_a+m_a)}{{\zeta }_a}}\left({\displaystyle \frac{{\beta }_a{\eta }_a{A}_s^{\ast }}{m({\delta }_a+m_a)({\zeta }_a+m_a)}}-1\right)A_i^{\ast }\hfill \end{array} \\ +{\displaystyle \frac{({\delta }_b+m_b)({\zeta }_b+m_b)}{{\zeta }_b}}\left({\displaystyle \frac{{\beta }_b{\eta }_b{B}_s^{\ast }}{m({\delta }_b+m_b)({\zeta }_b+m_b)}}-1\right)B_i^{\ast }=0.\end{array} \end{gathered}$$

Equation (12) implies that either $A_i^{\ast }=B_i^{\ast }=0$, or that

$\begin{array}{c}\hfill \begin{array}{c}{A}_s^{\ast }={\displaystyle \frac{m({\delta }_a+m_a)({\zeta }_a+m_a)}{{\beta }_a{\eta }_a{\zeta }_a}}={\displaystyle \frac{{\mathsf{\Lambda }}_a}{{\mathcal{R}}_1}}\mathrm{and}{B}_s^{\ast }={\displaystyle \frac{m({\delta }_b+m_b)({\zeta }_b+m_b)}{{\beta }_b{\eta }_b{\zeta }_b}}={\displaystyle \frac{{\mathsf{\Lambda }}_b}{{\mathcal{R}}_2}}.\hfill \end{array}\end{array}$

From the first equations of (7) and (8), we have

$\begin{array}{c}\hfill \begin{array}{c}\left(m_a{\mathsf{\Lambda }}_a-{\beta }_a{\displaystyle \frac{{\mathsf{\Lambda }}_a}{{\mathcal{R}}_1}}{\displaystyle \frac{{\eta }_a}{m}}{A}_i^{\ast }-{\displaystyle \frac{m_a{\mathsf{\Lambda }}_a}{{\mathcal{R}}_1}}\right)+\left(m_b{\mathsf{\Lambda }}_b-{\beta }_b{\displaystyle \frac{{\mathsf{\Lambda }}_b}{{\mathcal{R}}_2}}{\displaystyle \frac{{\eta }_b}{m}}{B}_i^{\ast }-{\displaystyle \frac{m_b{\mathsf{\Lambda }}_b}{{\mathcal{R}}_2}}\right)=0,\hfill \end{array}\end{array}$

and this gives $A_i^{\ast }={\displaystyle \frac{m{m}_a{\mathcal{R}}_1}{{\beta }_a{\eta }_a}}\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_1}}\right)$ and $B_i^{\ast }={\displaystyle \frac{m{m}_b{\mathcal{R}}_2}{{\beta }_b{\eta }_b}}\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_2}}\right)$. From (10), we have

$\begin{array}{c}\hfill \begin{array}{c}{A}_l^{\ast }={\displaystyle \frac{m{m}_a({\delta }_a+m_a){\mathcal{R}}_1}{{\beta }_a{\zeta }_a{\eta }_a}}\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_1}}\right),B_l^{\ast }={\displaystyle \frac{m{m}_b({\delta }_b+m_b){\mathcal{R}}_2}{{\beta }_b{\zeta }_b{\eta }_b}}\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_2}}\right),\hfill \end{array}\end{array}$

and $I_c^{\ast }={\displaystyle \frac{m_a{\mathcal{R}}_1}{{\beta }_a}}\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_1}}\right)+{\displaystyle \frac{m_b{\mathcal{R}}_2}{{\beta }_b}}\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_2}}\right)$. As we can note, ${\mathcal{E}}^{\ast }$ exists whenever $A_s^{\ast }>0$, $B_s^{\ast }>0$, $A_l^{\ast }>0$, $B_l^{\ast }>0$, $A_i^{\ast }>0$, $B_i^{\ast }>0$ and $I_c^{\ast }>0$, and this only feasible if ${\mathcal{R}}_0>1$. Therefore, we deduce that the endemic equilibrium point ${\mathcal{E}}^{\ast }$ exists, and is unique only if ${\mathcal{R}}_0>1$. □

Next, we determine the global stability of the Nosema disease-free equilibrium point and the endemic steady state. The global stability of these equilibria will be carried out by using the Lyapunov theory.

Theorem 2. 

Whenever ${\mathcal{R}}_0\le 1$, then ${\mathcal{E}}_0$ is globally asymptotically stable inside Σ.

Proof. 

Let us consider the candidate Lyapunov function,

$\begin{array}{c}\hfill \begin{array}{ccc}{\mathcal{L}}_0\hfill & =\hfill & {\zeta }_a{\eta }_a(m_b+{\zeta }_b)(m_b+{\delta }_b)A_l+{\eta }_a(m_a+{\zeta }_a)(m_b+{\zeta }_b)(m_b+{\delta }_b)A_i\hfill \\ & & +{\zeta }_b{\eta }_b(m_a+{\zeta }_a)(m_a+{\delta }_a)B_l+{\eta }_b(m_b+{\zeta }_b)(m_a+{\zeta }_a)(m_a+{\delta }_a)B_i\hfill \\ & & +(m_a+{\zeta }_a)(m_b+{\zeta }_b)(m_a+{\delta }_a)(m_b+{\delta }_b)I_c.\hfill \end{array}\end{array}$

By differentiating ${\mathcal{L}}_0$ along the trajectories of the dynamics of (3), we obtain

$\begin{array}{c}\hfill \begin{array}{ccc}\dot{{\mathcal{L}}_0}\hfill & =\hfill & {\zeta }_a{\eta }_a(m_b+{\zeta }_b)(m_b+{\delta }_b)\dot{A_l}+{\eta }_a(m_a+{\zeta }_a)(m_b+{\zeta }_b)(m_b+{\delta }_b)\dot{A_i}\hfill \\ & & +{\zeta }_b{\eta }_b(m_a+{\zeta }_a)(m_a+{\delta }_a)\dot{B_l}+{\eta }_b(m_b+{\zeta }_b)(m_a+{\zeta }_a)(m_a+{\delta }_a)\dot{B_i}\hfill \\ & & +(m_a+{\zeta }_a)(m_b+{\zeta }_b)(m_a+{\delta }_a)(m_b+{\delta }_b)\dot{I_c}.\hfill \end{array}\end{array}$

Since $A_s\le {\mathsf{\Lambda }}_a$ and $B_s\le {\mathsf{\Lambda }}_b$, it follows that

$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \dot{{\mathcal{L}}_0}\hfill & =\hfill & m(m_a+{\zeta }_a)(m_b+{\zeta }_b)(m_a+{\delta }_a)(m_b+{\delta }_b)[{\displaystyle \frac{{\beta }_a{\zeta }_a{\eta }_a{A}_s}{m(m_a+{\zeta }_a)(m_a+{\delta }_a)}}\hfill \\ & & \hfill +{\displaystyle \frac{{\beta }_b{\zeta }_b{\eta }_b{B}_s}{m(m_b+{\zeta }_b)(m_b+{\delta }_b)}}-1]I_c\hfill \\ & \le \hfill & m(m_a+{\zeta }_a)(m_b+{\zeta }_b)(m_a+{\delta }_a)(m_b+{\delta }_b)[{\displaystyle \frac{{\beta }_a{\zeta }_a{\eta }_a{\mathsf{\Lambda }}_a}{m(m_a+{\zeta }_a)(m_a+{\delta }_a)}}\hfill \\ & & \hfill +{\displaystyle \frac{{\beta }_b{\zeta }_b{\eta }_b{\mathsf{\Lambda }}_b}{m(m_b+{\zeta }_b)(m_b+{\delta }_b)}}-1]I_c\hfill \\ & =\hfill & m(m_a+{\zeta }_a)(m_b+{\zeta }_b)(m_a+{\delta }_a)(m_b+{\delta }_b)\left({\mathcal{R}}_0-1\right)I_c\hfill \\ & \le \hfill & 0\hfill \end{array}\hfill \end{array}$

where ${\mathcal{R}}_0\le 1$. When ${\mathcal{R}}_0<1$, $\dot{{\mathcal{L}}_0}=0$ yields $I_c=0$. Therefore, it can be deduced from System (3) that, as $t\to \infty$, $A_l\to 0$, $B_l\to 0$, $A_i\to 0$, $B_i\to 0$, $A_s\to {\mathsf{\Lambda }}_a$, and $B_s\to {\mathsf{\Lambda }}_b$. Thus, the invariant set for which $\dot{{\mathcal{L}}_0}=0$ is the singleton is ${\mathcal{E}}_0=\left({\mathsf{\Lambda }}_a,0,0,{\mathsf{\Lambda }}_b,0,0,0\right)$. By using Lasalle’s Invariance Principle [20], we deduce that any solution of the dynamics of (3) with initial conditions in $\mathsf{\Sigma }$, converges to ${\mathcal{E}}_0$ as $t\to \infty$. Thus, the equilibrium point ${\mathcal{E}}_0$ is globally asymptotically stable once ${\mathcal{R}}_0<1$. When ${\mathcal{R}}_0=1$, $\dot{{\mathcal{L}}_0}=0$ implies either $I_c=0$, or

$\begin{array}{c}\hfill \begin{array}{ccc}1\hfill & =\hfill & {\displaystyle \frac{{\beta }_a{\zeta }_a{\eta }_a{A}_s}{m(m_a+{\zeta }_a)(m_a+{\delta }_a)}}+{\displaystyle \frac{{\beta }_b{\zeta }_b{\eta }_b{B}_s}{m(m_b+{\zeta }_b)(m_b+{\delta }_b)}}\hfill \\ & \le \hfill & {\displaystyle \frac{{\beta }_a{\zeta }_a{\eta }_a{\mathsf{\Lambda }}_a}{m(m_a+{\zeta }_a)(m_a+{\delta }_a)}}+{\displaystyle \frac{{\beta }_b{\zeta }_b{\eta }_b{\mathsf{\Lambda }}_b}{m(m_b+{\zeta }_b)(m_b+{\delta }_b)}}\hfill \\ & =\hfill & {\mathcal{R}}_0.\hfill \end{array}\end{array}$

The latter case yields $A_s={\mathsf{\Lambda }}_a$, $B_s={\mathsf{\Lambda }}_b$ and, consequently, $A_l=A_i=B_l=B_i=I_c=0$. Therefore, the invariant set for which $\dot{{\mathcal{L}}_0}=0$ is the singleton ${\mathcal{E}}_0=({\mathsf{\Lambda }}_a,0,0,{\mathsf{\Lambda }}_b,0,0,0)$. Therefore, the steady state ${\mathcal{E}}_0$ is globally asymptotically stable if ${\mathcal{R}}_0=1$. □

In contrast, if ${\mathcal{R}}_0>1$, then by continuity, $\dot{{\mathcal{L}}_0}\left(t\right)>0$ in a neighborhood of ${\mathcal{E}}_0$ in $\mathsf{\Sigma }$. Trajectories in $\mathsf{\Sigma }$ sufficiently close to ${\mathcal{E}}_0$ move away from ${\mathcal{E}}_0$, implying that the DFE is unstable. By using the persistence theory [21] and an argument used in the proof of Proposition 3.3 of [22], it can be shown that, if ${\mathcal{R}}_0>1$, the instability of ${\mathcal{E}}_0$ implies uniform persistence of the dynamics of (3).

Next, we study the global stability of the endemic equilibrium point ${\mathcal{E}}^{\ast }$. By using the approach considered in [23], we impose the following assumptions for the dynamics of (3):

• (A1). For all $A_s$, $B_s$, $A_l$, $B_l$, $I_c>0$, $\left({\displaystyle \frac{I_c}{I_c^{\ast }}}-1\right)\left(1-{\displaystyle \frac{I_c^{\ast }{A}_l}{I_c{A}_l^{\ast }}}\right)\le 0$ and $\left({\displaystyle \frac{I_c}{I_c^{\ast }}}-1\right)\left(1-{\displaystyle \frac{I_c^{\ast }{B}_l}{I_c{B}_l^{\ast }}}\right)\le 0$.

• (A2). For all $A_l>0$, $B_l>0$, $A_i>0$, $B_i>0$, $\left({\displaystyle \frac{A_l}{A_l^{\ast }}}-1\right)\left(1-{\displaystyle \frac{A_l^{\ast }{A}_i}{A_l{A}_i^{\ast }}}\right)\le 0$, and $\left({\displaystyle \frac{B_l}{B_l^{\ast }}}-1\right)\left(1-{\displaystyle \frac{B_l^{\ast }{B}_i}{B_l{B}_i^{\ast }}}\right)\le 0$.

• (A3). For all $A_i>0$, $B_i>0$, $I_c>0$, $\left({\displaystyle \frac{A_i}{A_i^{\ast }}}-1\right)\left(1-{\displaystyle \frac{A_i^{\ast }{I}_c}{A_i{I}_c^{\ast }}}\right)\le 0$, and $\left({\displaystyle \frac{B_i}{B_i^{\ast }}}-1\right)\left(1-{\displaystyle \frac{B_i^{\ast }{I}_c}{B_i{I}_c^{\ast }}}\right)\le 0$.

These hypotheses are motivated by the proof of Theorem 6.1 of [23]. We now prove that if ${\mathcal{R}}_0>1$, then the endemic steady states is globally asymptotically stable.

Theorem 3. 

If ${\mathcal{R}}_0>1$, and under Assumptions (A1)–(A3), the endemic steady state of the dynamics of (3), ${\mathcal{E}}^{\ast }$, exists, and is globally asymptotically stable.

Proof. 

Let us consider the second candidate Lyapunov function,

$\begin{array}{c}\hfill \begin{array}{ccc}{\mathcal{L}}^{\ast }\left(t\right)\hfill & =\hfill & {\displaystyle {\int }_{A_s^{\ast }}^{A_s}{\displaystyle \frac{\sigma -A_s^{\ast }}{\sigma }}d\sigma +[A_l-A_l^{\ast }-A_l^{\ast }ln{\displaystyle \frac{A_l}{A_l^{\ast }}}]+[A_i-A_i^{\ast }-A_i^{\ast }ln{\displaystyle \frac{A_i}{A_i^{\ast }}}]}\hfill \\ & & {\displaystyle +{\int }_{B_s^{\ast }}^{B_s}{\displaystyle \frac{\sigma -B_s^{\ast }}{\sigma }}d\sigma +[B_l-B_l^{\ast }-B_l^{\ast }ln{\displaystyle \frac{B_l}{B_l^{\ast }}}]+[B_i-B_i^{\ast }-B_i^{\ast }ln{\displaystyle \frac{B_i}{B_i^{\ast }}}]}\hfill \\ & & {\displaystyle +[I_c-I_c^{\ast }-I_c^{\ast }ln{\displaystyle \frac{I_c}{I_c^{\ast }}}].}\hfill \end{array}\end{array}$

At the endemic steady state, we have the following identities:

(13)$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\mathsf{\Lambda }}_a={\beta }_a{A}_s^{\ast }{I}_c^{\ast }+m_a{A}_s^{\ast },{\mathsf{\Lambda }}_a={\beta }_b{B}_s^{\ast }{I}_c^{\ast }+m_b{B}_s^{\ast },({\zeta }_a+m_a)A_l^{\ast }={\beta }_a{A}_s^{\ast }{I}_c^{\ast },}\hfill \\ ({\zeta }_b+m_b)B_l^{\ast }={\beta }_b{B}_s^{\ast }{I}_c^{\ast },({\delta }_a+m_a)A_i^{\ast }={\zeta }_a{A}_l^{\ast },({\delta }_b+m_b)B_i^{\ast }={\zeta }_b{B}_l^{\ast },\hfill \\ \mathrm{and} m{I}_c^{\ast }={\eta }_a{A}_i^{\ast }+{\eta }_b{B}_i^{\ast }.\hfill \end{array}\end{array}$

By differentiating ${\mathcal{L}}^{\ast }$ along the solutions of the dynamics of (3) and using the equilibrium points Equation (13), we have

$\begin{array}{c}\hfill \begin{array}{ccc}{\displaystyle \frac{d{\mathcal{L}}^{\ast }}{dt}}\hfill & =\hfill & {\displaystyle [1-{\displaystyle \frac{A_s^{\ast }}{A_s}}]{\dot{A}}_s+[1-{\displaystyle \frac{B_s^{\ast }}{B_s}}]{\dot{B}}_s+[1-{\displaystyle \frac{A_l^{\ast }}{A_l}}]{\dot{A}}_l+[1-{\displaystyle \frac{B_l^{\ast }}{B_l}}]{\dot{B}}_l+[1-{\displaystyle \frac{A_i^{\ast }}{A_i}}]{\dot{A}}_i+}\hfill \\ & & [1-{\displaystyle \frac{B_i^{\ast }}{B_i}}]{\dot{B}}_i+[1-{\displaystyle \frac{I_c^{\ast }}{I_c}}]{\dot{I}}_c\hfill \\ & =\hfill & {\displaystyle [1-{\displaystyle \frac{A_s^{\ast }}{A_s}}][{\beta }_a{A}_s^{\ast }{I}_c^{\ast }+m_a{A}_s^{\ast }-{\beta }_a{A}_s{I}_c-m_a{A}_s]}\hfill \\ & & +[1-{\displaystyle \frac{B_s^{\ast }}{B_s}}][{\beta }_b{B}_s^{\ast }{I}_c^{\ast }+m_b{B}_s^{\ast }-{\beta }_b{B}_s{I}_c-m_b{B}_s]+{\eta }_a{A}_i-{\eta }_a{A}_i{\displaystyle \frac{I_c^{\ast }}{I_c}}-{\eta }_a{A}_i^{\ast }{\displaystyle \frac{I_c}{I_c^{\ast }}}\hfill \\ & & +{\eta }_a{A}_i^{\ast }+{\eta }_b{B}_i-{\eta }_b{B}_i{\displaystyle \frac{I_c^{\ast }}{I_c}}-{\eta }_b{B}_i^{\ast }{\displaystyle \frac{I_c}{I_c^{\ast }}}+{\eta }_b{B}_i^{\ast }+{\beta }_a{A}_s{I}_c-{\beta }_a{A}_s{I}_c{\displaystyle \frac{A_l^{\ast }}{A_l}}-{\beta }_a{A}_s^{\ast }{I}_c^{\ast }{\displaystyle \frac{A_l}{A_l^{\ast }}}\hfill \\ & & +{\beta }_a{A}_s^{\ast }{I}_c^{\ast }+{\beta }_b{B}_s{I}_c-{\beta }_b{B}_s{I}_c{\displaystyle \frac{B_l^{\ast }}{B_l}}-{\beta }_b{B}_s^{\ast }{I}_c^{\ast }{\displaystyle \frac{B_l}{B_l^{\ast }}}+{\beta }_b{B}_s^{\ast }{I}_c^{\ast }+{\zeta }_a{A}_l-{\zeta }_a{A}_l{\displaystyle \frac{A_i^{\ast }}{A_i}}\hfill \\ & & {\displaystyle -{\zeta }_a{A}_l^{\ast }{\displaystyle \frac{A_i}{A_i^{\ast }}}+{\zeta }_a{A}_l^{\ast }+{\zeta }_b{B}_l-{\zeta }_b{B}_l{\displaystyle \frac{B_i^{\ast }}{B_i}}-{\zeta }_b{B}_l^{\ast }{\displaystyle \frac{B_i}{B_i^{\ast }}}+{\zeta }_b{B}_l^{\ast }.}\hfill \end{array}\end{array}$

After some algebraic manipulations, we have

$\begin{array}{c}\hfill \begin{array}{ccc}{\displaystyle \frac{d{\mathcal{L}}^{\ast }}{dt}}\hfill & =\hfill & -m_a{\displaystyle \frac{(A_s-A_s^{\ast }{)}^2}{A_s}}-m_b{\displaystyle \frac{(B_s-B_s^{\ast }{)}^2}{B_s}}+{\beta }_a{A}_s^{\ast }{I}_c^{\ast }[{\displaystyle \frac{I_c}{I_c^{\ast }}}-1][1-{\displaystyle \frac{I_c^{\ast }{A}_l}{I_c{A}_l^{\ast }}}]\hfill \\ & & +{\beta }_b{B}_s^{\ast }{I}_c^{\ast }[{\displaystyle \frac{I_c}{I_c^{\ast }}}-1][1-{\displaystyle \frac{I_c^{\ast }{B}_l}{I_c{B}_l^{\ast }}}]+{\beta }_a{A}_s^{\ast }{I}_c^{\ast }[3-{\displaystyle \frac{A_s^{\ast }}{A_s}}-{\displaystyle \frac{A_s{I}_c{A}_l^{\ast }}{A_s^{\ast }{I}_c^{\ast }{A}_l}}-{\displaystyle \frac{I_c^{\ast }{A}_l}{I_c{A}_l^{\ast }}}]\hfill \\ & & +{\beta }_b{B}_s^{\ast }{I}_c^{\ast }[3-{\displaystyle \frac{B_s^{\ast }}{B_s}}-{\displaystyle \frac{B_s{I}_c{B}_l^{\ast }}{B_s^{\ast }{I}_c^{\ast }{B}_l}}-{\displaystyle \frac{I_c^{\ast }{B}_l}{I_c{B}_l^{\ast }}}]+{\zeta }_a{A}_l^{\ast }[{\displaystyle \frac{A_l}{A_l^{\ast }}}-1][1-{\displaystyle \frac{A_l^{\ast }{A}_i}{A_l{A}_i^{\ast }}}]\hfill \\ & & +{\zeta }_b{B}_l^{\ast }[{\displaystyle \frac{B_l}{B_l^{\ast }}}-1][1-{\displaystyle \frac{B_l^{\ast }{B}_i}{B_l{B}_i^{\ast }}}]+{\zeta }_a{A}_l^{\ast }[2-{\displaystyle \frac{A_l{A}_i^{\ast }}{A_l^{\ast }{A}_i}}-{\displaystyle \frac{A_l^{\ast }{A}_i}{A_l{A}_i^{\ast }}}]\hfill \\ & & +{\zeta }_b{B}_l^{\ast }[2-{\displaystyle \frac{B_l{B}_i^{\ast }}{B_l^{\ast }{B}_i}}-{\displaystyle \frac{B_l^{\ast }{B}_i}{B_l{B}_i^{\ast }}}]+{\eta }_a{A}_i^{\ast }[{\displaystyle \frac{A_i}{A_i^{\ast }}}-1][1-{\displaystyle \frac{A_i^{\ast }{I}_c}{A_i{I}_c^{\ast }}}]\hfill \\ & & +{\eta }_b{B}_i^{\ast }[{\displaystyle \frac{B_i}{B_i^{\ast }}}-1][1-{\displaystyle \frac{B_i^{\ast }{I}_c}{B_i{I}_c^{\ast }}}]+{\eta }_a{A}_i^{\ast }[2-{\displaystyle \frac{A_i{I}_c^{\ast }}{A_i^{\ast }{I}_c}}-{\displaystyle \frac{A_i^{\ast }{I}_c}{A_i{I}_c^{\ast }}}]\hfill \\ & & +{\eta }_b{B}_i^{\ast }[2-{\displaystyle \frac{B_i{I}_c^{\ast }}{B_i^{\ast }{I}_c}}-{\displaystyle \frac{B_i^{\ast }{I}_c}{B_i{I}_c^{\ast }}}].\hfill \end{array}\end{array}$

It follows from Assumptions (A1)–(A3) that

$\begin{array}{c}\hfill \begin{array}{ccc}{\displaystyle \frac{d{\mathcal{L}}^{\ast }}{dt}}\hfill & \le \hfill & {\beta }_a{A}_s^{\ast }{I}_c^{\ast }[3-{\displaystyle \frac{A_s^{\ast }}{A_s}}-{\displaystyle \frac{A_s{I}_c{A}_l^{\ast }}{A_s^{\ast }{I}_c^{\ast }{A}_l}}-{\displaystyle \frac{I_c^{\ast }{A}_l}{I_c{A}_l^{\ast }}}]+{\beta }_b{B}_s^{\ast }{I}_c^{\ast }[3-{\displaystyle \frac{B_s^{\ast }}{B_s}}-{\displaystyle \frac{B_s{I}_c{B}_l^{\ast }}{B_s^{\ast }{I}_c^{\ast }{B}_l}}-{\displaystyle \frac{I_c^{\ast }{B}_l}{I_c{B}_l^{\ast }}}]\hfill \\ & & +{\zeta }_a{A}_l^{\ast }[2-{\displaystyle \frac{A_l{A}_i^{\ast }}{A_l^{\ast }{A}_i}}-{\displaystyle \frac{A_l^{\ast }{A}_i}{A_l{A}_i^{\ast }}}]+{\zeta }_b{B}_l^{\ast }[2-{\displaystyle \frac{B_l{B}_i^{\ast }}{B_l^{\ast }{B}_i}}-{\displaystyle \frac{B_l^{\ast }{B}_i}{B_l{B}_i^{\ast }}}]\hfill \\ & & +{\eta }_a{A}_i^{\ast }[2-{\displaystyle \frac{A_i{I}_c^{\ast }}{A_i^{\ast }{I}_c}}-{\displaystyle \frac{A_i^{\ast }{I}_c}{A_i{I}_c^{\ast }}}]+{\eta }_b{B}_i^{\ast }[2-{\displaystyle \frac{B_i{I}_c^{\ast }}{B_i^{\ast }{I}_c}}-{\displaystyle \frac{B_i^{\ast }{I}_c}{B_i{I}_c^{\ast }}}]\hfill \\ & \le \hfill & 0.\hfill \end{array}\end{array}$

The six terms are non-positive, since the geometric mean is smaller than the arithmetic one, as follows:

$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle \frac{A_s^{\ast }}{A_s}}+{\displaystyle \frac{A_s{I}_c{A}_l^{\ast }}{A_s^{\ast }{I}_c^{\ast }{A}_l}}+{\displaystyle \frac{I_c^{\ast }{A}_l}{I_c{A}_l^{\ast }}}\ge 3\sqrt[3]{{\displaystyle \frac{A_s^{\ast }}{A_s}}\times {\displaystyle \frac{A_s{I}_c{A}_l^{\ast }}{A_s^{\ast }{I}_c^{\ast }{A}_l}}\times {\displaystyle \frac{I_c^{\ast }{A}_l}{I_c{A}_l^{\ast }}}}=3,\hfill \\ & {\displaystyle \frac{B_s^{\ast }}{B_s}}+{\displaystyle \frac{B_s{I}_c{B}_l^{\ast }}{B_s^{\ast }{I}_c^{\ast }{B}_l}}+{\displaystyle \frac{I_c^{\ast }{B}_l}{I_c{B}_l^{\ast }}}\ge 3\sqrt[3]{{\displaystyle \frac{B_s^{\ast }}{B_s}}\times {\displaystyle \frac{B_s{I}_c{B}_l^{\ast }}{B_s^{\ast }{I}_c^{\ast }{B}_l}}\times {\displaystyle \frac{I_c^{\ast }{B}_l}{I_c{B}_l^{\ast }}}}=3,\hfill \\ & {\displaystyle \frac{A_l{A}_i^{\ast }}{A_l^{\ast }{A}_i}}+{\displaystyle \frac{A_l^{\ast }{A}_i}{A_l{A}_i^{\ast }}}\ge 2\sqrt{{\displaystyle \frac{A_l{A}_i^{\ast }}{A_l^{\ast }{A}_i}}\times {\displaystyle \frac{A_l^{\ast }{A}_i}{A_l{A}_i^{\ast }}}}=2,{\displaystyle \frac{B_l{B}_i^{\ast }}{B_l^{\ast }{B}_i}}+{\displaystyle \frac{B_l^{\ast }{B}_i}{B_l{B}_i^{\ast }}}\ge 2\sqrt{{\displaystyle \frac{B_l{B}_i^{\ast }}{B_l^{\ast }{B}_i}}\times {\displaystyle \frac{B_l^{\ast }{B}_i}{B_l{B}_i^{\ast }}}}=2,\hfill \\ & {\displaystyle \frac{A_i{I}_c^{\ast }}{A_i^{\ast }{I}_c}}+{\displaystyle \frac{A_i^{\ast }{I}_c}{A_i{I}_c^{\ast }}}\ge 2\sqrt{{\displaystyle \frac{A_i{I}_c^{\ast }}{A_i^{\ast }{I}_c}}\times {\displaystyle \frac{A_i^{\ast }{I}_c}{A_i{I}_c^{\ast }}}}=2 \mathrm{and} {\displaystyle \frac{B_i{I}_c^{\ast }}{B_i^{\ast }{I}_c}}+{\displaystyle \frac{B_i^{\ast }{I}_c}{B_i{I}_c^{\ast }}}\ge 2\sqrt{{\displaystyle \frac{B_i{I}_c^{\ast }}{B_i^{\ast }{I}_c}}\times {\displaystyle \frac{B_i^{\ast }{I}_c}{B_i{I}_c^{\ast }}}}=2.\hfill \end{array}\end{array}$

It is easy to see that the largest invariant set for which $\dot{{\mathcal{L}}^{\ast }}=0$ is the singleton is {${\mathcal{E}}^{\ast }$}. Then, by using LaSalle’s invariance principle, we deduce the global asymptotic stability of ${\mathcal{E}}^{\ast }$. □

4. Influence of Seasonality on the Dynamics

Let us consider time-varying parameters reflecting the periodicity influence, and let us rewrite the dynamics of (3) as follows:

(14)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{A}}_l\left(t\right)={\beta }_a\left(t\right)A_s\left(t\right)I_c\left(t\right)-({\zeta }_a\left(t\right)+m_a\left(t\right))A_l\left(t\right),\hfill \\ {\dot{B}}_l\left(t\right)={\beta }_b\left(t\right)B_s\left(t\right)I_c\left(t\right)-({\zeta }_b\left(t\right)+m_b\left(t\right))B_l\left(t\right),\hfill \\ {\dot{A}}_i\left(t\right)={\zeta }_a\left(t\right)A_l\left(t\right)-({\delta }_a\left(t\right)+m_a\left(t\right))A_i\left(t\right),\hfill \\ {\dot{B}}_i\left(t\right)={\zeta }_b\left(t\right)B_l\left(t\right)-({\delta }_b\left(t\right)+m_b\left(t\right))B_i\left(t\right),\hfill \\ \dot{I_c}\left(t\right)={\eta }_a\left(t\right)A_i\left(t\right)+{\eta }_b\left(t\right)B_i\left(t\right)-m\left(t\right)I_c\left(t\right),\hfill \\ {\dot{A}}_s\left(t\right)=m_a\left(t\right){\mathsf{\Lambda }}_a\left(t\right)-{\beta }_a\left(t\right)A_s\left(t\right)I_c\left(t\right)-m_a\left(t\right)A_s\left(t\right),\hfill \\ {\dot{B}}_s\left(t\right)=m_b\left(t\right){\mathsf{\Lambda }}_b\left(t\right)-{\beta }_b\left(t\right)B_s\left(t\right)I_c\left(t\right)-m_b\left(t\right)B_s\left(t\right),\hfill \end{array}\right.\end{array}$

with the initial condition $(A_l^0,B_l^0,A_i^0,B_i^0,I_c^0,A_s^0,B_s^0)\in {\mathbb{R}}_{+}^7$.

Let us consider a positive T-periodic continuous function, $\upsilon$ and we denote by ${\displaystyle {\upsilon }^u=\underset{t\in [0,T)}{max}\upsilon \left(t\right)}$ and ${\displaystyle {\upsilon }^l=\underset{t\in [0,T)}{min}\upsilon \left(t\right)}$ the upper and lower bounds, respectively.

4.1. Preliminary

Let us consider an $m\times m$ continuous matrix function denoted by $Z\left(t\right)$, which is T-periodic, irreducible, and cooperative. Consider $X_Z\left(t\right)$; the solution of the equation

(15)$\begin{array}{c}\dot{X}\left(t\right)=Z\left(t\right)X\left(t\right),\hfill \end{array}$

and $r\left(X_Z\left(T\right)\right)$, the spectral radius of $X_Z\left(T\right)$, have positive elements $\forall t>0$. By applying the famous theory of Perron–Frobenius [24,25], one can deduce that $X_Z\left(T\right)$ has the principal eigenvalue $r\left(X_Z\left(T\right)\right)$. Therefore, we need to use the following lemma several times.

Lemma 2 

([26]). The ordinary differential Equation (15) admits the solution $X\left(t\right)=x\left(t\right)e^{\ell t}$, where $\ell ={\displaystyle \frac{1}{T}}ln\left(r\left(X_Z\left(T\right)\right)\right)$, and the function $x\left(t\right)$ is positive and T-periodic.

Let us consider the two-dimensional dynamics,

(16)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{A}}_s\left(t\right)=m_a\left(t\right)({\mathsf{\Lambda }}_a\left(t\right)-A_s\left(t\right)),\hfill \\ {\dot{B}}_s\left(t\right)=m_b\left(t\right)({\mathsf{\Lambda }}_b\left(t\right)-B_s\left(t\right)).\hfill \end{array}\right.\end{array}$

such that $(A_s^0,B_s^0)\in {\mathbb{R}}_{+}^2$. The dynamics of (16) have a unique T-periodic trajectory $({\overline{A}}_s\left(t\right),{\overline{B}}_s\left(t\right))$ that it is globally attractive in ${\mathbb{R}}_{+}^2$, such that ${\overline{A}}_s\left(t\right)>0$ and ${\overline{B}}_s\left(t\right)>0$. Thus, the dynamics of (14) have a unique Nosema disease-free periodic trajectory $\overline{U}\left(t\right)=(0,0,0,0,0,{\overline{A}}_s\left(t\right),{\overline{B}}_s\left(t\right))$.

Proposition 1. 

The positive compact set

$\begin{array}{c}\hfill \begin{array}{c}\hfill {\displaystyle {\mathsf{\Sigma }}^u=\{(A_l,B_l,A_i,B_i,I_c,A_s,B_s)\in {\mathbb{R}}_{+}^7/A_s+A_l+A_i\le {\mathsf{\Lambda }}_a^u,B_s+B_l+B_i\le {\mathsf{\Lambda }}_b^u,}\\ \hfill I_c\le {\displaystyle \frac{{\eta }_a^u{\mathsf{\Lambda }}_a^u+{\eta }_b^u{\mathsf{\Lambda }}_b^u}{m^l}}\}\end{array}\end{array}$

is an attractor of all trajectories of System (14) with

(17) $\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \underset{t\to \infty }{lim}{A}_s\left(t\right)+A_l\left(t\right)+A_i\left(t\right)-{\overline{A}}_s\left(t\right)=0,}\hfill \\ {\displaystyle \underset{t\to \infty }{lim}{B}_s\left(t\right)+B_l\left(t\right)+B_i\left(t\right)-{\overline{B}}_s\left(t\right)=0.}\hfill \end{array}\end{array}$

Proof. 

It is easy to see that if $A_s\left(t\right)+A_l\left(t\right)+A_i\left(t\right)\ge {\mathsf{\Lambda }}_a^u$, then

$\begin{array}{c}\hfill \begin{array}{c}{\dot{A}}_s\left(t\right)+{\dot{A}}_l\left(t\right)+{\dot{A}}_i\left(t\right)\le m_a\left(t\right)[{\mathsf{\Lambda }}_a\left(t\right)-A_s\left(t\right)-A_l\left(t\right)-A_i\left(t\right)]\le 0,\hfill \end{array}\end{array}$

and if $B_s\left(t\right)+B_l\left(t\right)+B_i\left(t\right)\ge {\mathsf{\Lambda }}_b^u$, then

$\begin{array}{c}\hfill \begin{array}{c}{\dot{B}}_s\left(t\right)+{\dot{B}}_l\left(t\right)+{\dot{B}}_i\left(t\right)\le m_b\left(t\right)[{\mathsf{\Lambda }}_b\left(t\right)-B_s\left(t\right)-B_l\left(t\right)-B_i\left(t\right)]\le 0.\hfill \end{array}\end{array}$

Similarly, if $I_c\ge {\displaystyle \frac{{\eta }_a^u{\mathsf{\Lambda }}_a^u+{\eta }_b^u{\mathsf{\Lambda }}_b^u}{m^l}}$, then we obtain

${\dot{I}}_c\left(t\right)=m\left(t\right)\left[{\displaystyle \frac{{\eta }_a\left(t\right)A_i\left(t\right)+{\eta }_b\left(t\right)B_i\left(t\right)}{m\left(t\right)}}-I_c\left(t\right)\right]\le 0.$

Let us define $T_a\left(t\right)=A_s\left(t\right)+A_l\left(t\right)+A_i\left(t\right)$ and $T_b\left(t\right)=B_s\left(t\right)+B_l\left(t\right)+B_i\left(t\right)$. Consider $x_a\left(t\right)=T_a\left(t\right)-{\overline{A}}_s\left(t\right),x_b\left(t\right)=T_b\left(t\right)-{\overline{B}}_s\left(t\right),t\ge 0$; then, ${\dot{x}}_a\left(t\right)\le -m_a\left(t\right)x_a\left(t\right)$ and ${\dot{x}}_b\left(t\right)\le -m_b\left(t\right)x_b\left(t\right)$ and, therefore, ${\displaystyle \underset{t\to \infty }{lim}{x}_a\left(t\right)=\underset{t\to \infty }{lim}(T_a\left(t\right)-{\overline{A}}_s\left(t\right))=0}$ and ${\displaystyle \underset{t\to \infty }{lim}{x}_b\left(t\right)=\underset{t\to \infty }{lim}(T_b\left(t\right)-{\overline{B}}_s\left(t\right))=0}$. □

4.2. Nosema Disease-Free Solution

Let us start by defining the basic reproduction number, ${\mathcal{R}}_0$, as given in [27]. Let $U=(A_l,B_l,A_i,B_i,I_c,A_s,B_s)$, $\mathcal{F}(t,U)=\left(\begin{array}{c}{\beta }_a\left(t\right)A_s\left(t\right)I_c\left(t\right)\\ {\beta }_b\left(t\right)B_s\left(t\right)I_c\left(t\right)\\ {\zeta }_a\left(t\right)A_l\left(t\right)\\ {\zeta }_b\left(t\right)B_l\left(t\right)\\ {\eta }_a\left(t\right)A_i\left(t\right)+{\eta }_b\left(t\right)B_i\left(t\right)\\ 0\\ 0\end{array}\right)$, ${\mathcal{V}}^{+}(t,U)=$$\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ m_a\left(t\right){\mathsf{\Lambda }}_a\left(t\right)\\ m_b\left(t\right){\mathsf{\Lambda }}_b\left(t\right)\end{array}\right)$, ${\mathcal{V}}^{-}(t,U)=\left(\begin{array}{c}({\zeta }_a\left(t\right)+m_a\left(t\right))A_l\left(t\right)\\ ({\zeta }_b\left(t\right)+m_b\left(t\right))B_l\left(t\right)\\ ({\delta }_a\left(t\right)+m_a\left(t\right))A_i\left(t\right)\\ ({\delta }_b\left(t\right)+m_b\left(t\right))B_i\left(t\right)\\ m\left(t\right)I_c\left(t\right)\\ m_a\left(t\right)A_s\left(t\right)\\ m_b\left(t\right)B_s\left(t\right)\end{array}\right)$, and $\mathcal{V}(t,U)={\mathcal{V}}^{-}(t,U)-{\mathcal{V}}^{+}$$(t,U)$. Therefore, the dynamics of (14) are expressed as follows:

(18)$\begin{array}{c}\hfill \begin{array}{c}\dot{U}=f(t,U\left(t\right))=\mathcal{F}(t,U)-\mathcal{V}(t,U).\hfill \end{array}\end{array}$

As we can see here, Assumptions (A1)–(A5) of Section 1 of [27] were satisfied. System (18) has a Nosema disease-free periodic solution $\overline{U}\left(t\right)=(0,0,0,0,0,{\overline{A}}_s\left(t\right),{\overline{B}}_s\left(t\right))$. Let us define $M\left(t\right)={\left({\displaystyle \frac{\partial f_i(t,\overline{U}\left(t\right))}{\partial U_j}}\right)}_{6\le i,j\le 7}$, where $f_i(t,U\left(t\right))$ and $U_i$ are the i-th components of $f(t,U(t\left)\right)$ and U, respectively. We can easily verify that $M\left(t\right)=\left(\begin{array}{cc}-m_a\left(t\right)& 0\\ 0& -m_b\left(t\right)\end{array}\right)$. Then, the spectral radius satisfies $r\left(X_M\left(T\right)\right)<1$, and the trajectory $\overline{U}\left(t\right)$ is linearly asymptotically stable inside ${\mathsf{\Gamma }}_s=\left\{(0,0,0,0,0,A_s\left(t\right),B_s\left(t\right))\in R_{+}^7\right\}$. Thus, we deduce that Assumption (A6) of Section 1 of [27] is also verified.

Note that the first five equations of System (14) describe the production of new infections and changes in the state among infected individuals. Therefore, let us consider the two matrices, ${\mathbf{A}}^{+}\left(t\right)$ and ${\mathbf{A}}^{-}\left(t\right)$, representing the rate of the appearance of new infections in these five equations, and the rate of transfer of individuals into and out of these compartments by all other means, respectively. ${\mathbf{A}}^{+}\left(t\right)$ and ${\mathbf{A}}^{-}\left(t\right)$ are given by ${\mathbf{A}}^{+}\left(t\right)={\left({\displaystyle \frac{\partial {\mathcal{F}}_i(t,\overline{U}\left(t\right))}{\partial U_j}}\right)}_{1\le i,j\le 5}$ and ${\mathbf{A}}^{-}\left(t\right)={\left({\displaystyle \frac{\partial {\mathcal{V}}_i(t,\overline{U}\left(t\right))}{\partial U_j}}\right)}_{1\le i,j\le 5}$. An easy calculus gives us

$\begin{array}{c}\hfill {\mathbf{A}}^{+}\left(t\right)=\left(\begin{array}{ccccc}0& 0& 0& 0& {\beta }_a\left(t\right){\mathsf{\Lambda }}_a\left(t\right)\\ 0& 0& 0& 0& {\beta }_b\left(t\right){\mathsf{\Lambda }}_b\left(t\right)\\ {\zeta }_a\left(t\right)& 0& 0& 0& 0\\ 0& {\zeta }_b\left(t\right)& 0& 0& 0\\ 0& 0& {\eta }_a\left(t\right)& {\eta }_b\left(t\right)& 0\end{array}\right),\\ \hfill {\mathbf{A}}^{-}\left(t\right)=\left(\begin{array}{ccccc}({\zeta }_a\left(t\right)+m_a\left(t\right))& 0& 0& 0& 0\\ 0& ({\zeta }_b\left(t\right)+m_b\left(t\right))& 0& 0& 0\\ 0& 0& ({\delta }_a\left(t\right)+m_a\left(t\right))& 0& 0\\ 0& 0& 0& ({\delta }_b\left(t\right)+m_b\left(t\right))& 0\\ 0& 0& 0& 0& m\left(t\right)\end{array}\right).\end{array}$

Consider $Z(A_s,B_s)$, the trajectory of the dynamics ${\displaystyle \frac{d}{dt}}Z(A_s,B_s)=-{\mathbf{A}}^{-}\left(A_s\right)Z(A_s,B_s)$ for any $A_s\ge B_s$, with $Z(A_s,A_s)=I_5$. We deduce, then, that Assumption (A7) of Section 1 of [27] is also satisfied.

Let us now define the following integral operator:

(19)$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \left(I\varphi \right)\left(\chi \right)={\int }_0^{\infty }Z(\chi ,\chi -w){\mathbf{A}}^{+}(\chi -w)\varphi (\chi -w)dw,\forall \chi \in \mathbb{R},\varphi \in {\mathsf{\Delta }}_T}\end{array}\end{array}$

where ${\mathsf{\Delta }}_T$ is the ordered Banach space of T-periodic functions defined on $\mathbb{R}$ to ${\mathbb{R}}^2$. Then, the basic reproduction number, ${\mathcal{R}}_0$, of Model (14) is the spectral radius of the operator I,

${\mathcal{R}}_0=r\left(I\right).$

Theorem 4. 

According to the theory in Theorem 2.2 of [27], the statements hereafter are valid:

• ${\mathcal{R}}_0<1\iff r\left(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}}\left(T\right)\right)<1$.

• ${\mathcal{R}}_0=1\iff r\left(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}}\left(T\right)\right)=1$.

• ${\mathcal{R}}_0>1\iff r\left(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}}\left(T\right)\right)>1$.

Therefore, the T-periodic trajectory $\overline{U}\left(t\right)$ is only asymptotically stable whenever ${\mathcal{R}}_0<1$. However, $\overline{U}\left(t\right)$ is unstable whenever ${\mathcal{R}}_0>1$. Furthermore, we aim to prove that whenever ${\mathcal{R}}_0<1$, then $\overline{U}\left(t\right)$ is globally asymptotically stable, which means that the Nosema disease goes to extinction.

Theorem 5. 

$\overline{U}\left(t\right)$ is globally asymptotically stable whenever ${\mathcal{R}}_0<1$, and it is unstable whenever ${\mathcal{R}}_0>1$.

Proof. 

The proof is similar to the one given in [28]. According to Theorem 4, $\overline{U}\left(t\right)$ is locally stable whenever ${\mathcal{R}}_0<1$, and it remains to prove the global attractivity whenever ${\mathcal{R}}_0<1$. According to the limit results of (17) given in Proposition 1, we deduce that, for any ${\sigma }_1>0$, there exists $T_1>0$, such that $A_l\left(t\right)+A_i\left(t\right)+A_s\left(t\right)\le {\overline{A}}_s\left(t\right)+{\sigma }_1$ and $B_l\left(t\right)+B_i\left(t\right)+B_s\left(t\right)\le {\overline{B}}_s\left(t\right)+{\sigma }_1$ for $t>T_1$. Therefore, $A_s\left(t\right)\le {\overline{A}}_s\left(t\right)+{\sigma }_1$, $B_s\left(t\right)\le {\overline{B}}_s\left(t\right)+{\sigma }_1$, and

(20)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{A}}_l\left(t\right)\le {\beta }_a\left(t\right)({\overline{A}}_s\left(t\right)+{\sigma }_1)I_c\left(t\right)-({\zeta }_a\left(t\right)+m_a\left(t\right))A_l\left(t\right),\hfill \\ {\dot{B}}_l\left(t\right)\le {\beta }_b\left(t\right)({\overline{B}}_s\left(t\right)+{\sigma }_1)I_c\left(t\right)-({\zeta }_b\left(t\right)+m_b\left(t\right))B_l\left(t\right),\hfill \\ {\dot{A}}_i\left(t\right)={\zeta }_a\left(t\right)A_l\left(t\right)-({\delta }_a\left(t\right)+m_a\left(t\right))A_i\left(t\right),\hfill \\ {\dot{B}}_i\left(t\right)={\zeta }_b\left(t\right)B_l\left(t\right)-({\delta }_b\left(t\right)+m_b\left(t\right))B_i\left(t\right),\hfill \\ \dot{I_c}\left(t\right)={\eta }_a\left(t\right)A_i\left(t\right)+{\eta }_b\left(t\right)B_i\left(t\right)-m\left(t\right)I_c\left(t\right),\hfill \end{array}\right.\end{array}$

for $t>T_1$. Let us define the matrix $D\left(t\right)$ as follows:

(21)$\begin{array}{c}\hfill D\left(t\right)=\left(\begin{array}{ccccc}0& 0& 0& 0& {\beta }_a\left(t\right)\\ 0& 0& 0& 0& {\beta }_b\left(t\right)\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right).\end{array}$

Since $r\left(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}}\left(T\right)\right)<1$, then ${\sigma }_1>0$ can be chosen to satisfy $r\left(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}+{\sigma }_{1}D}\left(T\right)\right)<1$. Let us define the two-dimensional system given hereafter:

(22)$\begin{array}{c}\hfill \left\{\begin{array}{ccc}{\dot{\overline{A}}}_l\left(t\right)\hfill & =\hfill & {\beta }_a\left(t\right)({\overline{A}}_s\left(t\right)+{\sigma }_1){\overline{I}}_c\left(t\right)-({\zeta }_a\left(t\right)+m_a\left(t\right)){\overline{A}}_l\left(t\right),\hfill \\ {\dot{\overline{B}}}_l\left(t\right)\hfill & =\hfill & {\beta }_b\left(t\right)({\overline{B}}_s\left(t\right)+{\sigma }_1){\overline{I}}_c\left(t\right)-({\zeta }_b\left(t\right)+m_b\left(t\right)){\overline{B}}_l\left(t\right),\hfill \\ {\dot{\overline{A}}}_i\left(t\right)\hfill & =\hfill & {\zeta }_a\left(t\right){\overline{A}}_l\left(t\right)-({\delta }_a\left(t\right)+m_a\left(t\right)){\overline{A}}_i\left(t\right),\hfill \\ {\dot{\overline{B}}}_i\left(t\right)\hfill & =\hfill & {\zeta }_b\left(t\right){\overline{B}}_l\left(t\right)-({\delta }_b\left(t\right)+m_b\left(t\right)){\overline{B}}_i\left(t\right),\hfill \\ {\dot{\overline{I}}}_c\left(t\right)\hfill & =\hfill & {\eta }_a\left(t\right){\overline{A}}_i\left(t\right)+{\eta }_b\left(t\right){\overline{B}}_i\left(t\right)-m\left(t\right){\overline{I}}_c\left(t\right).\hfill \end{array}\right.\end{array}$

From the theory in [26], there exists a positive T-periodic function $x_1\left(t\right)$, such that $w\left(t\right)\le x_1\left(t\right)e^{a_{1}t}$, where

$w\left(t\right)={\left(A_l\left(t\right),B_l\left(t\right),A_i\left(t\right),B_i\left(t\right),I_c\left(t\right)\right)}^T\mathrm{and}{a}_1={\displaystyle \frac{1}{T}}ln\left(r(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}+{\sigma }_{1}D}\left(T\right)\right)<0.$

Thus, ${\displaystyle \underset{t\to \infty }{lim}{A}_l\left(t\right)=\underset{t\to \infty }{lim}{B}_l\left(t\right)=\underset{t\to \infty }{lim}{A}_i\left(t\right)=\underset{t\to \infty }{lim}{B}_i\left(t\right)=0}$ and ${\displaystyle \underset{t\to \infty }{lim}{I}_c\left(t\right)=0}$. Furthermore, we have that ${\displaystyle \underset{t\to \infty }{lim}{A}_s\left(t\right)-{\overline{A}}_s\left(t\right)=\underset{t\to \infty }{lim}{T}_a\left(t\right)-{\overline{A}}_s\left(t\right)-{\overline{A}}_l\left(t\right)-{\overline{A}}_i\left(t\right)=0}$ and ${\displaystyle \underset{t\to \infty }{lim}{B}_s\left(t\right)-{\overline{B}}_s\left(t\right)=\underset{t\to \infty }{lim}{T}_b\left(t\right)-{\overline{B}}_s\left(t\right)-{\overline{B}}_l\left(t\right)-{\overline{B}}_i\left(t\right)=0}$. Therefore, $\overline{U}\left(t\right)$ satisfies the globally attractivity for ${\mathcal{R}}_0<1$. □

Now, we show that if ${\mathcal{R}}_0>1$, then $A_l\left(t\right)$, $A_i\left(t\right)$$B_l\left(t\right)$, and $B_i\left(t\right)$ are uniformly persistent, which means that the Nosema disease will persist.

4.3. Disease Persistence

Let us define $U_0=(A_l^0,B_l^0,A_i^0,B_i^0,I_c^0,A_s^0,B_s^0)$ and $U_1=(0,0,0,0,0,A_s^0,B_s^0)$, and consider the Poincaré map $Q:{\mathbb{R}}_{+}^7\to {\mathbb{R}}_{+}^7$ associated with Model (14), such that $U_0\mapsto u(T,U^0)$ is the trajectory of the dynamics of (14) with the initial value $u(0,U^0)=U^0\in {\mathbb{R}}_{+}^7$. Q is point dissipative. Consider the sets

$\mathsf{\Gamma }=\left\{(A_l,B_l,A_i,B_i,I_c,A_s,B_s)\in {\mathbb{R}}_{+}^7\right\},{\mathsf{\Gamma }}_0=Int\left({\mathbb{R}}_{+}^7\right),\partial {\mathsf{\Gamma }}_0=\mathsf{\Gamma }\setminus {\mathsf{\Gamma }}_0 \mathrm{and}$

${\mathsf{\Gamma }}_{\partial }=\left\{U_0\in \partial {\mathsf{\Gamma }}_0:Q^p\left(U_0\right)\in \partial {\mathsf{\Gamma }}_0,\forall p\ge 0\right\}.$

$\mathsf{\Gamma }$ and ${\mathsf{\Gamma }}_0$ are invariant and, according to the theory in [26,29], we deduce that

(23)$\begin{array}{c}\hfill \begin{array}{c}{\mathsf{\Gamma }}_{\partial }=\left\{(0,0,0,0,0,A_s,B_s),A_s\ge 0,B_s\ge 0\right\}\hfill \end{array}\end{array}$

with ${\mathsf{\Gamma }}_{\partial }\supseteq \left\{(0,0,0,0,0,A_s,B_s),A_s\ge 0,B_s\ge 0\right\}$.

It remains to prove that ${\mathsf{\Gamma }}_{\partial }\setminus \left\{(0,0,0,0,0,A_s,B_s),A_s\ge 0,B_s\ge 0\right\}=\varnothing$. Consider $U_0\in {\mathsf{\Gamma }}_{\partial }\setminus \left\{(0,0,0,0,0,A_s,B_s),A_s\ge 0,B_s\ge 0\right\}$. If $A_l\left(0\right)=0$ and $0<B_l\left(0\right)$, then ${\dot{A}}_l{\left(t\right)}_{|t=0}={\beta }_a\left(0\right)A_s\left(0\right)I_c\left(0\right)>0$. If $A_l\left(0\right)>0$ and $B_l\left(0\right)=0$, then $A_l\left(t\right),B_s\left(t\right)>0$, $\forall t>0$. Therefore, we have, for any $t>0$,

$\begin{array}{c}\hfill \begin{array}{c}\hfill {\displaystyle B_l\left(t\right)=\left[B_l\left(0\right)+{\int }_0^t\left({\beta }_b\left(\omega \right)B_s\left(\omega \right)I_c\left(\omega \right)\right)e^{{\displaystyle {\int }_0^{\omega }({\zeta }_b\left(z\right)+m_b\left(z\right))dz}}d\omega \right]\times }\hfill \\ \hfill e^{{\displaystyle -{\int }_0^t({\zeta }_b\left(z\right)+m_b\left(z\right))dz}}>0\hfill \end{array}\end{array}$

which implies that $U\left(t\right)\notin \partial {\mathsf{\Gamma }}_0$ if $0<t\ll 1$ and ${\mathsf{\Gamma }}_0$ is positively invariant, and (23) is then then satisfied. Then, Q has a unique fixed point $U_1$ inside ${\mathsf{\Gamma }}_{\partial }$. We provide the main result of this section.

Theorem 6. 

Once ${\mathcal{R}}_0>1$, Model (14) admits at least one periodic solution, such that there exists $\eta >0$, verifying $\forall U_0\in Int{\left({\mathbb{R}}_{+}\right)}^5\times {\mathbb{R}}_{+}^2$, and we have

${\displaystyle \underset{t\to \infty }{lim\; inf}{A}_l\left(t\right)\ge \eta >0,\underset{t\to \infty }{lim\; inf}{B}_l\left(t\right)\ge \eta >0,\underset{t\to \infty }{lim\; inf}{A}_i\left(t\right)\ge \eta >0,\underset{t\to \infty }{lim\; inf}{B}_i\left(t\right)\ge \eta >0,}$

${\displaystyle \mathit{and}\underset{t\to \infty }{lim\; inf}{I}_c\left(t\right)\ge \eta >0.}$

Proof. 

We aim to demonstrate that the solutions of (14) are uniformly persistent when associated with $({\mathsf{\Gamma }}_0,\partial {\mathsf{\Gamma }}_0)$ through the Poincaré map properties given in Theorem 3.1.1 of [29]. We have $r\left(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}}\left(T\right)\right)>1$; therefore, $\exists {\sigma }_2>0$, such that $r\left(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}-{\sigma }_{2}D}\left(T\right)\right)>1$. Let consider the following two-dimensional dynamics:

(24)$\begin{array}{c}\hfill \left\{\begin{array}{ccc}{\dot{A}}_s^{\epsilon }\left(t\right)\hfill & =\hfill & m_a\left(t\right){\mathsf{\Lambda }}_a\left(t\right)-{\beta }_a\left(t\right)A_s^{\epsilon }\left(t\right)\epsilon -m_a\left(t\right)A_s^{\epsilon }\left(t\right),\hfill \\ {\dot{B}}_s^{\epsilon }\left(t\right)\hfill & =\hfill & m_b\left(t\right){\mathsf{\Lambda }}_b\left(t\right)-{\beta }_b\left(t\right)B_s^{\epsilon }\left(t\right)\epsilon -m_b\left(t\right)B_s^{\epsilon }\left(t\right).\hfill \end{array}\right.\end{array}$

The Poincaré map Q, associated with the dynamics of (24), admits a unique fixed point $({\overline{A}}_s^{\epsilon },{\overline{B}}_s^{\epsilon })$ that it is globally attractive. Let us apply the implicit function theorem; we deduce that $\epsilon \mapsto ({\overline{A}}_s^{\epsilon },{\overline{B}}_s^{\epsilon })$ is continuous. Thus, by choosing $\epsilon >0$ small enough, we can obtain ${\overline{A}}_s^{\epsilon }\left(t\right)>{\overline{A}}_s\left(t\right)-{\sigma }_2,$ and ${\overline{B}}_s^{\epsilon }\left(t\right)>{\overline{B}}_s\left(t\right)-{\sigma }_2,$$\forall t>0$. Now, as the trajectory is continuous respecting $U_0$, then $\exists {\epsilon }^{\ast }>0$, such that $\parallel U_0-u(t,U_1)\parallel \le {\epsilon }^{\ast }$; thus,

$\parallel u(t,U_0)-u(t,U_1)\parallel <\epsilon \mathrm{for} \mathrm{all} t\in [0,T].$

We aim to prove that

(25)$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \underset{p\to \infty }{lim\; sup}d(Q^p\left(U_0\right),U_1)\ge {\epsilon }^{\ast }\forall U_0\in {\mathsf{\Gamma }}_0}\hfill \end{array}\end{array}$

by contradiction. Assume that ${\displaystyle \underset{p\to \infty }{lim\; sup}d(Q^p\left(U_0\right),U_1)<{\epsilon }^{\ast }}$ for some $U_0\in {\mathsf{\Gamma }}_0$. We can assume that $d(Q^p\left(U_0\right),U_1)<{\epsilon }^{\ast },\forall p>0$. Therefore, we deduce that

$\parallel u(t,Q^p\left(U_0\right))-u(t,U_1)\parallel <\epsilon ,\forall p>0,t\in [0,T].$

For any $t\ge 0$, it can be expressed as follows: $t=pT+t_1$, where $t_1\in [0,T)$ and $p=\lfloor {\displaystyle \frac{t}{T}}\rfloor$. Then,

$\parallel u(t,U_0)-u(t,U_1)\parallel =\parallel u(t_1,Q^p\left(U_0\right))-u(t_1,U_1)\parallel <\epsilon ,\forall t\ge 0.$

Set $(A_l\left(t\right),B_l\left(t\right),A_i\left(t\right),B_i\left(t\right),I_c\left(t\right),A_s\left(t\right),B_s\left(t\right))=u(t,U_0)$.

Therefore, $0\le A_l\left(t\right),B_l\left(t\right),A_i\left(t\right),B_i\left(t\right),I_c\left(t\right)\le \epsilon ,t\ge 0$, and

(26)$\begin{array}{c}\hfill \left\{\begin{array}{ccc}{\dot{A}}_s\left(t\right)\hfill & \ge \hfill & m_a\left(t\right){\mathsf{\Lambda }}_a\left(t\right)-{\beta }_a\left(t\right)A_s\left(t\right)\epsilon -m_a\left(t\right)A_s\left(t\right),\hfill \\ {\dot{B}}_s\left(t\right)\hfill & \ge \hfill & m_b\left(t\right){\mathsf{\Lambda }}_b\left(t\right)-{\beta }_b\left(t\right)B_s\left(t\right)\epsilon -m_b\left(t\right)B_s\left(t\right).\hfill \end{array}\right.\end{array}$

Note that Q, associated with the dynamics of (24), admits a fixed point $({\overline{A}}_s^{\epsilon },{\overline{B}}_s^{\epsilon })$, which is globally attractive where ${\overline{A}}_s^{\epsilon }\left(t\right)>{\overline{A}}_s-{\sigma }_2,$ and ${\overline{B}}_s^{\epsilon }\left(t\right)>{\overline{B}}_s\left(t\right)-{\sigma }_2$; then, $\exists T_2>0$, satisfying $A_s\left(t\right)>{\overline{A}}_s\left(t\right)-{\sigma }_2$ and $B_s\left(t\right)>{\overline{B}}_s\left(t\right)-{\sigma }_2$ for $t>T_2$. Then, for $t>T_2$, we have

(27)$\begin{array}{c}\hfill \left\{\begin{array}{ccc}{\dot{A}}_l\left(t\right)\hfill & \ge \hfill & {\beta }_a\left(t\right)({\overline{A}}_s\left(t\right)-{\sigma }_2)I_c\left(t\right)-({\zeta }_a\left(t\right)+m_a\left(t\right))A_l\left(t\right),\hfill \\ {\dot{B}}_l\left(t\right)\hfill & \ge \hfill & {\beta }_b\left(t\right)({\overline{B}}_s\left(t\right)-{\sigma }_2)I_c\left(t\right)-({\zeta }_b\left(t\right)+m_b\left(t\right))B_l\left(t\right),\hfill \\ {\dot{A}}_i\left(t\right)\hfill & =\hfill & {\zeta }_a\left(t\right)A_l\left(t\right)-({\delta }_a\left(t\right)+m_a\left(t\right))A_i\left(t\right),\hfill \\ {\dot{B}}_i\left(t\right)\hfill & =\hfill & {\zeta }_b\left(t\right)B_l\left(t\right)-({\delta }_b\left(t\right)+m_b\left(t\right))B_i\left(t\right),\hfill \\ \dot{I_c}\left(t\right)\hfill & =\hfill & {\eta }_a\left(t\right)A_i\left(t\right)+{\eta }_b\left(t\right)B_i\left(t\right)-m\left(t\right)I_c\left(t\right).\hfill \end{array}\right.\end{array}$

Since $r\left(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}-{\sigma }_{2}D}\left(T\right)\right)>1$, we deduce the existence of a T-periodic positive function $x\left(t\right)$[26], such that $J\left(t\right)\ge e^{at}x\left(t\right)$, where $a={\displaystyle \frac{1}{T}}lnr\left(X_{{\mathbf{A}}^{+}-{\mathbf{A}}^{-}-{\sigma }_{2}D}\left(T\right)\right)>0$; thus, ${\displaystyle \underset{t\to \infty }{lim}{U}_i\left(t\right)=\infty }$, which is impossible, since the trajectory is bounded. Therefore, (25) is verified, and the map Q is weakly uniformly persistent with respect to $({\mathsf{\Gamma }}_0,\partial {\mathsf{\Gamma }}_0)$. By referring to Proposition 1, we deduce that the map Q has a global attractor. Thus, $U_1$ is an isolated invariant set of $\mathsf{\Gamma }$, and $W^s\left(U_1\right)\cap {\mathsf{\Gamma }}_0=\varnothing$. Therefore, every trajectory in ${\mathsf{\Gamma }}_{\partial }$ must converge to $U_1$, which is an acyclic in ${\mathsf{\Gamma }}_{\partial }$. By using Theorem 1.3.1 and Remark 1.3.1 from [29], we obtain the uniform persistence of the map Q with respect to $({\mathsf{\Gamma }}_0,\partial {\mathsf{\Gamma }}_0)$. Moreover, by referring to Theorem 1.3.6 of [29], Q admits a fixed point ${\tilde{U}}_0=({\tilde{A}}_l,{\tilde{B}}_l,{\tilde{A}}_i,{\tilde{B}}_i,{\tilde{I}}_c,{\tilde{A}}_s,{\tilde{B}}_s)\in {\mathsf{\Gamma }}_0$ with ${\tilde{U}}_0\in Int{\left(R_{+}\right)}^5\times R_{+}^2$.

Assume that ${\tilde{A}}_s=0$, and by injecting this into the dynamics of (14), ${\tilde{A}}_s\left(t\right)$ verifies

(28)$\begin{array}{c}\hfill \begin{array}{ccc}{\dot{\tilde{A}}}_s\left(t\right)\hfill & =\hfill & m_a\left(t\right){\mathsf{\Lambda }}_a\left(t\right)-{\beta }_a\left(t\right){\tilde{A}}_s\left(t\right){\tilde{I}}_c\left(t\right)-m_a\left(t\right){\tilde{A}}_s\left(t\right),\hfill \end{array}\end{array}$

with ${\tilde{A}}_s={\tilde{A}}_s\left(nT\right)=0,n=1,2,3,\cdots$. By using Proposition 1, $\forall {\sigma }_3>0$ and $\exists T_3>0$, satisfying

${\tilde{I}}_c\left(t\right)\le {\displaystyle \frac{{\eta }_a^u{\mathsf{\Lambda }}_a^u+{\eta }_b^u{\mathsf{\Lambda }}_b^u}{m^l}}+{\sigma }_3\mathrm{for}t>T_3.$

Therefore, we obtain

(29)$\begin{array}{c}\hfill \begin{array}{c}{\dot{\tilde{A}}}_s\left(t\right)\ge m_a\left(t\right){\mathsf{\Lambda }}_a\left(t\right)-{\beta }_a\left(t\right){\tilde{A}}_s\left(t\right)({\displaystyle \frac{{\eta }_a^u{\mathsf{\Lambda }}_a^u+{\eta }_b^u{\mathsf{\Lambda }}_b^u}{m^l}}+{\sigma }_3)-m_a\left(t\right){\tilde{A}}_s\left(t\right), \mathrm{for} t\ge T_3.\hfill \end{array}\end{array}$

$\exists \overline{n}>0$ such that $nT>T_3,\forall n>\overline{n}$. Therefore, we deduce that

$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\displaystyle {\tilde{A}}_s\left(nT\right)}\hfill & \ge \hfill & {\displaystyle \left[{\tilde{A}}_s\left(0\right)+{\int }_0^{nT}{m}_a\left(z\right){\mathsf{\Lambda }}_a\left(z\right)e^{{\displaystyle {\int }_0^z\left({\beta }_a\left(t\right)({\displaystyle \frac{{\eta }_a^u{\mathsf{\Lambda }}_a^u+{\eta }_b^u{\mathsf{\Lambda }}_b^u}{m^l}}+{\sigma }_3)+m_a\left(t\right)\right)dt}}dz\right]}\hfill \\ & & \hfill \times e^{{\displaystyle -{\int }_0^{nT}\left({\beta }_a\left(t\right)({\displaystyle \frac{{\eta }_a^u{\mathsf{\Lambda }}_a^u+{\eta }_b^u{\mathsf{\Lambda }}_b^u}{m^l}}+{\sigma }_3)+m_a\left(t\right)\right)dt}}\hfill \end{array}\end{array}$

for all $n>\overline{n}$, and this is impossible because ${\tilde{A}}_s\left(nT\right)=0$. Then, ${\tilde{A}}_s\left(0\right)>0$ and ${\tilde{U}}_0$ is a T-periodic trajectory of the dynamics of (14). □

5. Numerical Investigation

Our goal in this section is to perform the theoretical findings concerning the dynamics of (3) through some numerical examples. The model parameters will be periodic functions of the form $p\left(t\right)=p_0(1+p_{1}cos\left(q\pi (t+\alpha )\right))$, where $q\in \mathbb{N}$, $0\le p_0$, $0\le p_1\le 1$, and $0\le \alpha \le 1$ is the phase angle. Therefore, we define the model parameters as follows:

$\left\{\begin{array}{cccccc}\hfill {\mathsf{\Lambda }}_a\left(t\right)& =\hfill & {\mathsf{\Lambda }}_{a0}(1+{\mathsf{\Lambda }}_{a1}cos\left(q\pi (t+\alpha )\right)),\hfill & m_a\left(t\right)\hfill & \hfill =& m_{a0}(1+m_{a1}cos\left(q\pi (t+\alpha )\right)),\hfill \\ \hfill {\mathsf{\Lambda }}_b\left(t\right)& =\hfill & {\mathsf{\Lambda }}_{b0}(1+{\mathsf{\Lambda }}_{b1}cos\left(q\pi (t+\alpha )\right)),\hfill & m_b\left(t\right)\hfill & \hfill =& m_{b0}(1+m_{b1}cos\left(q\pi (t+\alpha )\right)),\hfill \\ \hfill {\beta }_a\left(t\right)& =\hfill & {\beta }_{a0}(1+{\beta }_{a1}cos\left(q\pi (t+\alpha )\right)),\hfill & {\zeta }_a\left(t\right)\hfill & \hfill =& {\zeta }_{a0}(1+{\zeta }_{a1}cos\left(q\pi (t+\alpha )\right)),\hfill \\ \hfill {\beta }_b\left(t\right)& =\hfill & {\beta }_{b0}(1+{\beta }_{b1}cos\left(q\pi (t+\alpha )\right)),\hfill & {\zeta }_b\left(t\right)\hfill & \hfill =& {\zeta }_{b0}(1+{\zeta }_{b1}cos\left(q\pi (t+\alpha )\right)),\hfill \\ \hfill {\eta }_a\left(t\right)& =\hfill & {\eta }_{a0}(1+{\eta }_{a1}cos\left(q\pi (t+\alpha )\right)),\hfill & {\delta }_a\left(t\right)\hfill & \hfill =& {\delta }_{a0}(1+{\delta }_{a1}cos\left(q\pi (t+\alpha )\right)),\hfill \\ \hfill {\eta }_b\left(t\right)& =\hfill & {\eta }_{b0}(1+{\eta }_{b1}cos\left(q\pi (t+\alpha )\right)),\hfill & {\delta }_b\left(t\right)\hfill & \hfill =& {\delta }_{b0}(1+{\delta }_{b1}cos\left(q\pi (t+\alpha )\right)),\hfill \\ \hfill m\left(t\right)& =\hfill & m_0(1+m_{1}cos\left(q\pi (t+\alpha )\right)).\hfill \end{array}\right.$

The seasonal cycle frequencies ${\mathsf{\Lambda }}_{a1}$, ${\mathsf{\Lambda }}_{b1}$, ${\beta }_{a1}$, ${\beta }_{b1}$, ${\eta }_{a1}$, ${\eta }_{b1}$, $m_1$, $m_{a1}$, $m_{b1}$, ${\zeta }_{a1}$, ${\zeta }_{b1}$, ${\delta }_{a1}$, and ${\delta }_{b1}$ are the amplitudes satisfying $|{\mathsf{\Lambda }}_{a1}|<1$, $|{\mathsf{\Lambda }}_{b1}|<1$, $|{\beta }_{a1}|<1$, $|{\beta }_{b1}|<1$, $|{\eta }_{a1}|<1$, $|{\eta }_{b1}|<1$, $|m_1|<1$, $|m_{a1}|<1$, $|m_{b1}|<1$, $|{\zeta }_{a1}|<1$, $|{\zeta }_{b1}|<1$, $|{\delta }_{a1}|<1$, and $|{\delta }_{b1}|<1$. The constants ${\beta }_{a0}$, ${\beta }_{b0}$, ${\eta }_{a0}$, ${\eta }_{b0}$, $m_0$, $m_{a0}$, $m_{b0}$, ${\zeta }_{a0}$, ${\zeta }_{b0}$, ${\delta }_{a0}$, ${\delta }_{b0}$, ${\mathsf{\Lambda }}_{a1}$, ${\mathsf{\Lambda }}_{b1}$, ${\beta }_{a1}$, ${\beta }_{b1}$, ${\eta }_{a1}$, ${\eta }_{b1}$, $m_1$, $m_{a1}$, $m_{b1}$, ${\zeta }_{a1}$, ${\zeta }_{b1}$, ${\delta }_{a1}$, ${\delta }_{b1}$, $\alpha$, and q are given in Table 2. Unfortunately, we have no biological data to use for our simulations. The parameter values considered for the numerical simulations are chosen arbitrarily.

The numerical method used here to resolve our system was the Runge–Kutta–Fehlberg method (RKF45), which is a numerical method of order 4 with an error estimator of order 5, and was applied under MATLAB R2024a Software (ode45). The first set of tests concerns the case of constant parameters. The second set of examples concerns the case of only T-periodic variable recruitment ${\mathsf{\Lambda }}_a\left(t\right)$ and ${\mathsf{\Lambda }}_b\left(t\right)$. The third set of examples concerns the case where all of the model parameters are periodic functions. In the first step, we consider the case where all parameters are assumed to be constant. We provide several examples validating the theoretical findings concerning the behavior of the trajectories of Model (3).

In Figure 2 (left), we provide the trajectories where ${\mathsf{\Lambda }}_{a0}=3$ and ${\mathsf{\Lambda }}_{b0}=2$, such that ${\mathcal{R}}_0\approx 4.8664>1$. The solution converges to the endemic steady state. However, in Figure 2 (right), we consider a set of parameters ${\mathsf{\Lambda }}_{a0}=0.4$ and ${\mathsf{\Lambda }}_{b0}=0.3$, such that ${\mathcal{R}}_0\approx 0.6802<1$. The trajectory converges to the HIV-free steady state ${\mathcal{E}}_0=\left({\mathsf{\Lambda }}_{a0},0,0,{\mathsf{\Lambda }}_{b0},0,0,0\right)$. In Figure 3, we provide the trajectories for several initial values, which converge to the same endemic steady state for the case where ${\mathsf{\Lambda }}_{a0}=3$ and ${\mathsf{\Lambda }}_{b0}=2$, and thus ${\mathcal{R}}_0\approx 4.8664>1$. In Figure 4, we provide the trajectories for several initial values which converge to the same Nosema disease-free steady state for the case where ${\mathsf{\Lambda }}_{a0}=0.4$ and ${\mathsf{\Lambda }}_{b0}=0.3$, and thus ${\mathcal{R}}_0\approx 0.6802<1$.

In the second step, we perform numerical simulations on the dynamics of (3) where only the seasonally forced T-periodic functions ${\mathsf{\Lambda }}_a\left(t\right)$ and ${\mathsf{\Lambda }}_b\left(t\right)$ depend on time. The approximation the basic reproduction number, ${\mathcal{R}}_0$, for the periodic non-autonomous epidemic models can be computed from the corresponding time-averaged autonomous epidemic models in several cases [30,31,32,33]. Several other approximations of ${\mathcal{R}}_0$ are used in [25,34]. We provide several examples validating the theoretical findings concerning the behavior of the trajectories of Model (3).

In Figure 5 (left), we provide the trajectories where ${\mathsf{\Lambda }}_{a0}=3$ and ${\mathsf{\Lambda }}_{b0}=2$, such that ${\mathcal{R}}_0\approx 4.8664>1$. The solution converges to the periodic trajectory expressing the persistence of the Nosema disease. However, in Figure 5 (right), we consider a set of parameters ${\mathsf{\Lambda }}_{a0}=0.4$ and ${\mathsf{\Lambda }}_{b0}=0.3$, such that ${\mathcal{R}}_0\approx 0.6802<1$. The trajectory converges to the Nosema disease-free steady state $\left({\overline{A}}_s\left(t\right),0,0,{\overline{B}}_s\left(t\right),0,0,0\right)$. In Figure 6, we provide the trajectories for several initial values, which converge to the same periodic trajectory for the case where ${\mathsf{\Lambda }}_{a0}=3$ and ${\mathsf{\Lambda }}_{b0}=2$, and thus ${\mathcal{R}}_0\approx 4.8664>1$. In Figure 7, we provide the trajectories for several initial values which converge to the same Nosema disease-free steady state for the case where ${\mathsf{\Lambda }}_{a0}=0.4$ and ${\mathsf{\Lambda }}_{b0}=0.3$, and thus ${\mathcal{R}}_0\approx 0.6802<1$.

In the last case, we assume that all the model parameters are periodic functions reflecting a full periodic environment. Again, the approximation of ${\mathcal{R}}_0$ is performed using the time-averaged dynamics. We provide several examples validating the theoretical findings concerning the behavior of the trajectories of Model (3).

In Figure 8 (left), we provide the trajectories where ${\mathsf{\Lambda }}_{a0}=3$ and ${\mathsf{\Lambda }}_{b0}=2$, such that ${\mathcal{R}}_0\approx 4.8664>1$. The solution converges to the periodic trajectory expressing the persistence of the Nosema disease. However, in Figure 8 (right), we consider a set of parameters ${\mathsf{\Lambda }}_{a0}=0.4$ and ${\mathsf{\Lambda }}_{b0}=0.3$, such that ${\mathcal{R}}_0\approx 0.6802<1$. The trajectory converges to the Nosema disease-free steady state $\left({\overline{A}}_s\left(t\right),0,0,{\overline{B}}_s\left(t\right),0,0,0\right)$. In Figure 9, we provide the trajectories for several initial values, which converge to the same periodic trajectory for the case where ${\mathsf{\Lambda }}_{a0}=3$ and ${\mathsf{\Lambda }}_{b0}=2$, and thus ${\mathcal{R}}_0\approx 4.8664>1$. In Figure 10, we provide the trajectories for several initial values which converge to the same Nosema disease-free steady state for the case where ${\mathsf{\Lambda }}_{a0}=0.4$ and ${\mathsf{\Lambda }}_{b0}=0.3$, and thus ${\mathcal{R}}_0\approx 0.6802<1$.

6. Conclusions

In this article, a differential equation-based model for two different beehives sharing the same space and infected by the Nosema disease, is proposed and analyzed. The proposed model is concerned with different host populations at the two colonies, but with a common space. It is deduced that, while each colony contribute differently to the transmission dynamics, the overall reproduction number, ${\mathcal{R}}_0$, is determined by the sum of the two individual reproduction numbers, and provide a threshold for the Nosema disease dynamics: if ${\mathcal{R}}_0\le 1$, the Nosema disease-free steady state is globally asymptotically stable, indicating that the infection would die out; however, if ${\mathcal{R}}_0>1$, the trajectories converge to a limit cycle reflecting the disease persistence.

Author Contributions

Conceptualization, M.E.H. and F.A.S.A.; methodology, M.E.H., F.A.S.A. and M.H.A.; software, M.E.H., F.A.S.A. and M.H.A.; investigation, M.E.H., F.A.S.A. and M.H.A.; visualization, M.E.H., F.A.S.A. and M.H.A.; writing—original draft, M.E.H., F.A.S.A. and M.H.A.; writing—review and editing, M.E.H., F.A.S.A. and M.H.A. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

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

Acknowledgments

This work was funded by the University of Jeddah, Jeddah, Saudi Arabia, under grant No. (UJ-24-FR-20823-1). Therefore, the authors thank the University of Jeddah for its technical and financial support. The authors are also grateful to the unknown referees for the many constructive suggestions, which helped to improve the presentation of the paper.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Nosema disease transmission diagram.

Figure 2. Behavior of trajectories of Model (3) for the case where [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], and thus [Forumla omitted. See PDF.] (left), and for the case where [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], and thus [Forumla omitted. See PDF.] (right).

Figure 3. The variable dynamics for several initial values for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], such that [Forumla omitted. See PDF.]. The different colors describe the different initial conditions.

Figure 4. The variable dynamics for several initial values for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], such that [Forumla omitted. See PDF.]. The different colors describe the different initial conditions.

Figure 5. Behavior of trajectories of model (3) for the case where [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], and thus [Forumla omitted. See PDF.] (left), and for the case where [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], and thus [Forumla omitted. See PDF.] (right).

Figure 6. The variable dynamics for several initial values for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], such that [Forumla omitted. See PDF.]. The different colors describe the different initial conditions.

Figure 7. The variable dynamics for several initial values for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], such that [Forumla omitted. See PDF.]. The different colors describe the different initial conditions.

Figure 8. Behavior of trajectories of Model (3) for the case where [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], and thus [Forumla omitted. See PDF.] (left), and for the case where [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], and thus [Forumla omitted. See PDF.] (right).

Figure 9. The variable dynamics for several initial values for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], such that [Forumla omitted. See PDF.]. The different colors describe the different initial conditions.

Figure 10. The variable dynamics for several initial values for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], such that [Forumla omitted. See PDF.]. The different colors describe the different initial conditions.

Table 1

Description of model parameters.

Var.	Signification	Var.	Signification
${\beta }_a$	Transmission rate for the first beehive	${\delta }_a$	Recovery rate for the first beehive
${\beta }_b$	Transmission rate for the second beehive	${\delta }_b$	Recovery rate for the second beehive
$m_a$	Natural death rate for the first beehive	${\mathsf{\Lambda }}_a$	Birth rate for the first beehive
$m_b$	Natural death rate for the second beehive	${\mathsf{\Lambda }}_b$	Birth rate for the second beehive
${\zeta }_a$	Incubation rate for the first beehive	${\eta }_a$	Parasite shedding rate from $A_i$
${\zeta }_b$	Incubation rate for the second beehive	${\eta }_b$	Parasite shedding rate from $B_i$
m	Parasite death rate

Table 2

Constants of the model parameters.

${\mathit{\beta }}_{\mathit{a}0}$	${\mathit{\beta }}_{\mathit{b}0}$	${\mathit{\eta }}_{\mathit{a}0}$	${\mathit{\eta }}_{\mathit{b}0}$	${\mathit{m}}_0$	${\mathit{m}}_{\mathit{a}0}$	${\mathit{m}}_{\mathit{b}0}$	${\mathit{\zeta }}_{\mathit{a}0}$	${\mathit{\zeta }}_{\mathit{b}0}$	${\mathit{\delta }}_{\mathit{a}0}$	${\mathit{\delta }}_{\mathit{b}0}$	$\mathit{\alpha }$	q
$1.4$	$1.5$	$0.34$	$0.33$	$0.27$	$0.3$	$0.4$	$1.35$	$1.42$	$1.15$	$1.12$	4	4
${\mathbf{\Lambda }}_{\mathit{a}\mathbf{1}}$	${\mathbf{\Lambda }}_{\mathit{b}\mathbf{1}}$	${\mathbf{\beta }}_{\mathit{a}\mathbf{1}}$	${\mathbf{\beta }}_{\mathit{b}\mathbf{1}}$	${\mathbf{\eta }}_{\mathit{a}\mathbf{1}}$	${\mathbf{\eta }}_{\mathit{b}\mathbf{1}}$	${\mathit{m}}_{\mathbf{1}}$	${\mathit{m}}_{\mathit{a}\mathbf{1}}$	${\mathit{m}}_{\mathit{b}\mathbf{1}}$	${\mathit{\zeta }}_{\mathit{a}\mathbf{1}}$	${\mathit{\zeta }}_{\mathit{b}\mathbf{1}}$	${\mathit{\delta }}_{\mathit{a}\mathbf{1}}$	${\mathit{\delta }}_{\mathit{b}\mathbf{1}}$
$0.13$	$0.19$	$0.43$	$0.45$	$0.42$	$0.44$	$0.1$	$0.12$	$0.18$	$0.18$	$0.18$	$0.37$	$0.38$

References

1. Brauer, F.; Castillo-Chavez, C.; Mubayi, A.; Towers, S. Some models for epidemics of vector-transmitted diseases. Infect. Dis. Model.; 2016; 1, pp. 79-87. [DOI: https://dx.doi.org/10.1016/j.idm.2016.08.001]

2. Sasmal, S.K.; Ghosh, I.; Huppert, A.; Chattopadhyay, J. Modeling the Spread of Zika Virus in a Stage-Structured Population: Effect of Sexual Transmission. Bull. Math. Biol.; 2016; 80, pp. 3038-3067. [DOI: https://dx.doi.org/10.1007/s11538-018-0510-7]

3. Dénes, A.; Ibrahim, M.A. Global dynamics of a mathematical model for a honeybee colony infested by virus-carrying Varroa mites. J. Appl. Math. Comput.; 2019; 61, pp. 349-371. [DOI: https://dx.doi.org/10.1007/s12190-019-01250-5]

4. Ibrahim, M.A.; Dénes, A. Threshold Dynamics in a Model for Zika Virus Disease with Seasonality. Bull. Math. Biol.; 2021; 83, 27. [DOI: https://dx.doi.org/10.1007/s11538-020-00844-6]

5. El Hajji, M.; Alharbi, N.S.; Alharbi, M.H. Mathematical Modeling for a CHIKV Transmission Under the Influence of Periodic Environment. Int. J. Anal. Appl.; 2024; 22, 6. [DOI: https://dx.doi.org/10.28924/2291-8639-22-2024-6]

6. El Hajji, M. Periodic solutions for chikungunya virus dynamics in a seasonal environment with a general incidence rate. AIMS Math.; 2023; 8, pp. 24888-24913. [DOI: https://dx.doi.org/10.3934/math.20231269]

7. Najim, F.A.A. Mathematical Analysis for a Zika Virus Dynamics in a Seasonal Environment. Int. J. Anal. Appl.; 2024; 22, 71. [DOI: https://dx.doi.org/10.28924/2291-8639-22-2024-71]

8. Almuashi, H.H. Mathematical analysis for the influence of seasonality on Chikungunya virus dynamics. Int. J. Anal. Appl.; 2024; 22, 86. [DOI: https://dx.doi.org/10.28924/2291-8639-22-2024-86]

9. Alshammari, B.S.; Mashat, D.S.; Mallawi, F.O. Mathematical and Numerical Investigations for a Cholera Dynamics with a Seasonal Environment. Int. J. Anal. Appl.; 2023; 21, 127. [DOI: https://dx.doi.org/10.28924/2291-8639-21-2023-127]

10. El Hajji, M.; Alnjrani, R.M. Periodic Behaviour of HIV Dynamics with Three Infection Routes. Mathematics; 2024; 12, 123. [DOI: https://dx.doi.org/10.3390/math12010123]

11. El Hajji, M. Periodic solutions for an “SVIQR” epidemic model in a seasonal environment with general incidence rate. Int. J. Biomath.; 2024; Online ready [DOI: https://dx.doi.org/10.1142/S1793524524500335]

12. El Hajji, M.; Aloufi, M.F.S.; Alharbi, M.H. Influence of seasonality on Zika virus transmission. AIMS Math.; 2024; 9, pp. 19361-19384. [DOI: https://dx.doi.org/10.3934/math.2024943]

13. Ibrahim, M.A.; Dénes, A. A mathematical model for the spread of Varroa mites in honeybee populations: Two simulation scenarios with seasonality. Heliyon; 2022; 8, e10648. [DOI: https://dx.doi.org/10.1016/j.heliyon.2022.e10648] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36148272]

14. Chen, J.; Rincon, J.; DeGrandi-Hoffman, G.; Fewell, J.; Harrison, J.; Kang, Y. Impacts of seasonality and parasitism on honey bee population dynamics. J. Math. Biol.; 2023; 87, 19. [DOI: https://dx.doi.org/10.1007/s00285-023-01952-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37389742]

15. El Hajji, M.; Alzahrani, F.A.S.; Mdimagh, R. Impact of Infection on Honeybee Population Dynamics in a Seasonal Environment. Int. J. Anal. Appl.; 2024; 22, 75. [DOI: https://dx.doi.org/10.28924/2291-8639-22-2024-75]

16. Tokarev, Y.S.; Huang, W.F.; Solter, L.F.; Malysh, J.M.; Becnel, J.J.; Vossbrinck, C.R. A formal redefinition of the genera Nosema and Vairimorpha (Microsporidia: Nosematidae) and reassignment of species based on molecular phylogenetics. J. Invertebr. Pathol.; 2020; 169, 107279. [DOI: https://dx.doi.org/10.1016/j.jip.2019.107279]

17. Abban, S.; Smith, B.; Corona, M.; Cook, S.; Evans, J.; Chen, Y.; Alburaki, M. Prevalence and distribution of Varroa destructor and Nosema spp. in symptomatic honey bee colonies across the USA from 2015 to 2022. Sci. Rep.; 2024; 14, 1726. [DOI: https://dx.doi.org/10.1038/s41598-024-51514-9]

18. Mushayabasa, S.; Losio, A.A.E.; Modnak, C.; Wang, J. Optimal control analysis applied to a two-patch model for Guinea worm disease. Electron. J. Differ. Equ.; 2020; 2020, pp. 1-23. [DOI: https://dx.doi.org/10.58997/ejde.2020.70]

19. den Driessche, P.V.; Watmough, J. Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Math. Biosci.; 2002; 180, pp. 29-48. [DOI: https://dx.doi.org/10.1016/S0025-5564(02)00108-6]

20. LaSalle, J. The Stability of Dynamical Systems; SIAM: Philadelphia, PA, USA, 1976.

21. Freedman, H.I.; So, J.W.H.; Waltman, P. Coexistence in a model of competition in the chemostat incorporating discrete delays. SIAM J. Appl. Math.; 1989; 49, pp. 859-870. [DOI: https://dx.doi.org/10.1137/0149050]

22. Graef, S.P.; Andrews, J.F. Stability and control of anaerobic digestion. J. Water Pollut. Control. Fed.; 1974; 46, pp. 667-682.

23. Eisenberg, M.C.; Shuai, Z.; Tien, J.H.; van den Driessche, P. A cholera model in a patchy environment with water and human movement. Math. Biosci.; 2013; 246, pp. 105-112. [DOI: https://dx.doi.org/10.1016/j.mbs.2013.08.003]

24. Frobenius, G. Uber Matrizen aus Nicht Negativen Elementen; Sitzungsberichte Preussische Akademie der Wissenschaft: Berlin, Germany, 1912; pp. 456-477. [DOI: https://dx.doi.org/10.3931/e-rara-18865]

25. Nakata, Y.; Kuniya, T. Global dynamics of a class of SEIRS epidemic models in a periodic environment. J. Math. Anal. Appl.; 2010; 363, pp. 230-237. [DOI: https://dx.doi.org/10.1016/j.jmaa.2009.08.027]

26. Zhang, F.; Zhao, X. A periodic epidemic model in a patchy environment. J. Math. Anal. Appl.; 2007; 325, pp. 496-516. [DOI: https://dx.doi.org/10.1016/j.jmaa.2006.01.085]

27. Wang, W.; Zhao, X. Threshold dynamics for compartmental epidemic models in periodic environments. J. Dynam. Differ. Equ.; 2008; 20, pp. 699-717. [DOI: https://dx.doi.org/10.1007/s10884-008-9111-8]

28. El Hajji, M.; Alshaikh, D.M.; Almuallem, N.A. Periodic behaviour of an epidemic in a seasonal environment with vaccination. Mathematics; 2023; 11, 2350. [DOI: https://dx.doi.org/10.3390/math11102350]

29. Zhao, X.Q. Dynamical Systems in Population Biology. CMS Books Math./Ouvrages Math; SMC Springer-Verlag: New York, NY, USA, 2003; Volume 16.

30. Al-arydah, M.; Berhe, H.; Dib, K.; Madhu, K. Mathematical modeling of the spread of the coronavirus under strict social restrictions. Math. Methods Appl. Sci.; 2021; [DOI: https://dx.doi.org/10.1002/mma.7965]

31. Ma, J.; Ma, Z. Epidemic threshold conditions for seasonally forced SEIR models. Math. Biosci. Eng.; 2006; 3, pp. 161-172. [DOI: https://dx.doi.org/10.3934/mbe.2006.3.161]

32. Zhang, T.; Teng, Z. On a nonautonomous SEIRS model in epidemiology. Bull. Math. Biol.; 2007; 69, pp. 2537-2559. [DOI: https://dx.doi.org/10.1007/s11538-007-9231-z]

33. Wesley, C.L.; Allen, L.J. The basic reproduction number in epidemic models with periodic demographics. J. Biol. Dyn.; 2009; 3, pp. 116-129. [DOI: https://dx.doi.org/10.1080/17513750802304893]

34. Guerrero-Flores, S.; Osuna, O.; de Leon, C.V. Periodic solutions for seasonal SIQRS models with nonlinear infection terms. Electron. Differ. Equ.; 2019; 2019, pp. 1-13.