How to model honeybee population dynamics: stage structure and seasonality

Western honeybees (Apis Mellifera) serve extremely important roles in our ecosystem and economics as they are responsible for pollinating $ 215 billion dollars annually over the world. Unfortunately, honeybee population and their colonies have been declined dramatically. The purpose of this article is to explore how we should model honeybee population with age structure and validate the model using empirical data so that we can identify different factors that lead to the survival and healthy of the honeybee colony. Our theoretical study combined with simulations and data validation suggests that the proper age structure incorporated in the model and seasonality are important for modeling honeybee population. Specifically, our work implies that the model assuming that (1) the adult bees are survived from the {egg population} rather than the brood population; and (2) seasonality in the queen egg laying rate, give the better fit than other honeybee models. The related theoretical and numerical analysis of the most fit model indicate that (a) the survival of honeybee colonies requires a large queen egg-laying rate and smaller values of the other life history parameter values in addition to proper initial condition; (b) both brood and adult bee populations are increasing with respect to the increase in the {egg-laying rate} and the decreasing in other parameter values; and (c) seasonality may promote/suppress the survival of the honeybee colony.


Introduction
Western honeybee (Apis Mellifera) is an eusocial insect that has an advanced level of social organization. In honeybee colony, the queen produces the offspring and non-reproductive individuals cooperate in caring for the young ones, that forms complex colonies [45]. Honeybees play indispensable and important roles in human life, economy, and agriculture. For example, honeybees not only produce valuable products, such as honey, royal jelly, bee wax and propolis in the market, but also are responsible for pollinating crops such as blueberries, cherries, and almonds, that is worth $215 billion annually worldwide [43]. If there is no honeybee, it likely leads to changes in human diets and a disproportionate expansion of agricultural land in order to fill this shortfall in crop production by volume [35]. Unfortunately, honeybee population has been decreasing globally [43]. In the United States, the total number of honeybee colonies has been reduced approximate 40% to 50%, while in the rest of the world the total number of colonies is reduced by 5% to 10% [32]. The important and critical causes for honeybee colony mortalities include diseases, land-use change, pesticides, pathogens and parasites, and poor beekeeping management [8,9,31,33,43]. The purpose of this article is to explore how we could better model colony population dynamics to help us understand the honeybee colony mortalities.
Honeybee colony itself is a complex adaptive system with its own resilience to disturbances, whose survival depends on its individual quality, its adaptive capacity and its threshold of resilience to pressures [15]. On average, a colony has about 10,000 to 60,000 bees, that consists of a queen (fertile female) who produces all offspring, a few hundred drones (males) and thousands of workers (sterile females). Generally, a queen may lay approximate 1000-2000 eggs per day in the peak period [6]. Due to aging or disability of the queen bee, beekeepers will replace the queen every 1-2 years [6]. Each honeybee goes through four stages of development: egg, larva, pupa and adult [6]. For worker bees, they need 21 days to eclosion to adult bees [10,16,45], and drones need 24 days to mature [10]. Population size at each stage and the related maturation time have huge influences on the colony development and its population dynamics [15]. Needless to say, age is linked to division of labor in honeybees [38]. Young workers In the colony, young workers prefer to perform nursing tasks, while older workers prefer foraging activities. However, colonies can accelerate, delay, or even reverse their recruitment behavior as the internal or external environment changes [18].
Not only the age structure will affect the honeybees colony, but the change of season, temperature, weather, etc. also will influence the honeybees [6,10,20,41]. Through experiments and observations, honeybee population present periodic fluctuations due to different reasons. For instance, we observe great foraging activity during spring, summer and fall but the highest activity during the summer [6]. During spring and summer, pollen and nectar from diverse floras are in great abundance, giving rise to an increase honeybee population. Therefore, given that temperature is one of the main factors in honeybee food availability and thus brood production, honeybee population size is smaller during the winter [10,41]. Thus, the peak of the population is achieved in late June until middle of summer as it starts to decline [37]. The temperature in the colony also will influence honeybee, middle-age honeybees will respond to the heat stress in order to perform [20]. Thus, it is very important to include age structure and the seasonality in studying of honeybee population dynamics and the factors that affect healthy of honeybee colonies. Research has shown that the major problems threatening the survival of honeybee colonies could link to: 1) environmental stressors, such as habitat destruction (urbanization, deforestation, forest fires); 2) parasites and pathogens, such varroosis and virus; 3) genetic variation and vitality, like limited importation [31,33,43]. In order to quantify the problems and consider the difficulty of directly observing the dynamics of bee populations, mathematical models can be a powerful tool to help us understand how the bee population change and predict the fate of the colony.
Mathematical models indeed have been developed to study bee populations dynamics and the related stressors, particularly the effects of pathogens, parasites and nutrient stress factors [1, 2, 5, 11, 22-25, 33, 39, 43]. DeGrandi-Hoffman [10] proposed a first simulation model for honeybee colony dynamics that includes many important factors such weather, egg-laying rate, the age of queen, foraging and brood life cycles. There are some previous work focusing on how the death rate of foragers impacts colony viability [23,24]. Khoury [24] published a compartmental model based on these circumstances. The model includes three states, brood, hive and foragers, and incorporates the recruitment process to study the forager death rate. There is a work [39] that investigated seasonal food availability and transition of hive to foragers. The most recent recent works [29,36] consider the seasonality in the queen egg-laying rate. Messan et al [29] applied seasonality effects into the pollen collection rate that has annual periodicity by the first order harmonic. Ratti et al [36] also agrees that seasonality affects dynamics of honeybee and its parasitic virus. This article [36] incorporated seasonality in varroa treatment control as the treatment is applied with four seasons: spring, summer, fall, and winter [36].
Motivated by the previous work on honeybee population models with age structure [22][23][24] and seasonality [29,36,39], we propose and study honeybee population models with different delay terms to include age structure. We use data to validate our models and explore which model would be more appreciated and the importance of incorporating seasonality in the honeybee population model. More specifically, the objective of our paper is to develop a proper honeybee population dynamical model with age structure to understand important factors for colony survival, and to explore how seasonality may affect the colony dynamics and its survival.
The remaining of the article is shown as follows: In Section 2, we derive two honeybee colony dynamics models that incorporate varied delay terms. In Section 3, we perform rigorous mathematical analysis for those two models and compare their dynamics. In Section 4, we validate our two models with real honeybee data. Our study shows the importance of seasonality and suggests that one of those two proposed models would be more appropriated for studying honeybee population dynamics. In Section 5, we conclude our study. In the last section, we provide detailed proofs of our theoretical results.

Model Derivations
In this section, we focus on modeling of honeybee colony dynamics with age structure. For convenience, we divide the population of honeybee colony into brood and adult bees. Let B(t), H(t) be the population of brood and adult bees in a given hive at time t, respectively. We assume that: A1: The daily egg laying rate of honeybee queen is r with the survival rate of H 2 K+H 2 +αB where the parameter K is the population of adult bee needed for half of the maximum brood survival rate and α represents the regulation effects from brood population B. The term H 2 K+H 2 +αB reflects (1) the cooperative brood care from adult bees that perform nursing and collecting food for brood; and (2) the queen and workers that regulate the actual egg laying/survival rate based on the current available brood population B, which has been supported by the literature work [12,21,29,40].
A2: We assume that both brood and adult bees have constant mortality, d b and d h respectively.
The maturation time from brood B to adult bee H is denoted by τ (τ = 16 for queen, τ = 21 for workers, and τ = 24 for drones [6,23]), thus the maturation rate is termed as follows: The two assumptions above lead to the following non-linear delayed differential equations of honeybee population dynamics (Model (1)): where we assume that the initial condition for H(t) is a nonnegative continuous function when t ∈ [−τ, 0] and B(0) = 0 −τ rH 2 (s)e d b s K+H 2 (s)+αB(s) ds. The biological meaning of each parameter of the proposed model (1) is listed in Table 1. In the case α = 0, the model (1) reduces to the following Notes: Our proposed model (2) (when α = 0 in the model (1)) is a single specie model with brood B and adult H stage where these two stages seem to be decoupled. Thus, we could study the dynamics of Model (2) by exploring the dynamics of H first, then the dynamics of B is totally determined by H. We would see the analytical results in the next section.
In literature (e.g, see [44]), researchers have been using the compartmental models through ODEs to model population dynamics with age and/or stage structure. Motivated by this, we have the following delay model  where the term e −d b τ B(t − τ ) describes the maturation entry rate coming from the juvenile stage with a survival rate e −d b τ during the juvenile period τ .
More specifically, the model (3) above assumes that the adult H(t) matures from the survived brood population at time t − τ . In the following two sessions, we will compare population dynamics of Model (3) and the proposed model (1) to address the importances of deriving proper delay population models due to the outcomes of dynamics and model validations.

Mathematical Analysis
The state space of the proposed model (1) is X = C([−τ, 0], R + ) × C([−τ, 0], R + ). We first show that the proposed model (1) is positive invariant and bounded in X as the following theorem: Notes: The detailed proof of Theorem 3.1 is in the last section. Theorem 3.1 implies that our proposed model is biologically well defined. The model (1) always has the extinction equilibrium E e = (0, 0) which would be locally or globally stable as stated in the next theorem:  (1). Then it satisfies that the following equations: which gives and Solving the equation (7) gives with H * 1 ≤ H * 2 . Now we have the following proposition: where H * 1 is an increasing function of α, K, d b and d h , and H * 2 is a decreasing function of α, K, d b and d h ; whereas, H * 1 is a decreasing function of r, and H * 2 is an increasing function of r. In the case that has an unique interior equilibrium Notes: Proposition 3.1 implies that one of the necessary conditions for the honeybee colony survival > 1 which requires large values of the queen egg laying rate r, and the smaller values of the maturation time τ and the brood regulation effect α. In addition, Proposition 3.1 indicates that at the interior equilibrium, the ratio of brood B to adult population H is determined by their mortality and maturation time through the equation . Based on simulations and analytical results, the interior equilibrium H * 2 is always locally stable if it exists while H * 1 is locally unstable. If dα < 0. This implies that the brood regulation coefficient α has negative effects on brood and adult population sizes. In the case that α = 0, then if K > 1 holds, the interior equilibria H * i , i = 1, 2 have the following expressions: In order to study the stability of the interior equilibrium E i = (B * , H * ), we start with the characteristic equation of the interior equilibrium E i = B * i , H * i as follows by letting A = We can see the characteristic equation (10) of the interior equilibrium E i = (B * i , H * I ) is very complicated and difficult to analysis. Thus, for convenience, we start with the simpler case by setting α = 0 which is our model (2). K > 1, then Model (2) has two interior equilibria E i where E 1 is always unstable and E 2 is always locally asymptotically stable.
Notes: Theorem 3.3 indicates that the value of the maturation time τ has no effects on the stability of its interior equilibria E i , i = 1, 2 for Model (2). In the case that The following theorem provides results on the interior equilibrium stability of this critical case: Figure 1 shows bifurcation diagrams of Model (1) regarding (a) the queen egg-laying rate (r) (see Figure 1a&1b); (b) the brood regulation effects on reproduction α (see Figure 1c&1d); (c) the half-saturation coefficient K (see Figure 1e&1d) (d) the mortality of brood d b (see Figure 1g&1h); and (e) the mortality of adult d h (see Figure 1i&1j). Those bifurcation diagrams indicates that (1) the colony survival requires the large value of the queen egg-laying rate (r) which leads to the increased brood and adult population as it increases; (2) the large values of α, the half-saturation coefficient K, or any mortality rate d b or d h can lead to the colony collapsing, and both brood and adult population are decreasing with respect to these parameter values; and (3) increasing the value of the adult population mortality can lead to the dramatic decreasing of the adult population.
The another modeling approach with age structure for the honeybee colony: In our model derivation section, we proposed the model (3) below assuming that the adult H(t) matures from the survived brood population at time t − τ .
The model above is motivated from the compartmental ODE model in the literature [44]. We aim to compare the dynamics of Model (3) to the model (1) to address the importance of deriving a proper biological model with age structure.     First, we notice that the extinction equilibrium E e = (0, 0) always exists as for the model (1). However, E e can go through stability switching that leads to an oscillatory solution around E e for Model (3) Theorem 3.5. [Extinction equilibria dynamics] Model (3) always has the extinction equilibrium E e = (0, 0).
Notes: Theorem 3.5 suggests that the smaller value of the brood mortality can destablize the colony dynamics. In addition, it implies that Model (3) is not positive invariant as the extinction equilibrium E e = (0, 0) could have stability switches that lead to oscillatory solution around E e . See  (3), then it satisfies that the following two equations:

Let (B, H) be an interior equilibrium of Model
which gives the brood population at the equilibrium B * = d h e d b τ H * and by solving rH Now define the characteristic equation of the interior equilibrium (B * , H * ) of Model (3) as follows: (3) has two positive interior equilibrium is always stable or occurs stability switching by following cases: , then the stability of E 2 switches just once from stable to unstable as τ increases in [0, τ * ).
, the stability of E 2 can change a finite number of times at most as τ is increased τ ∈ [0, τ * ), and eventually it becomes unstable.
from stable to unstable.  Figure 4 also indicates the importances of initial condition that may lead to the survival or collapsing of the colony.
Comparisons between Model (1) and Model (3): Both models can have up to three equilibria with always the existence of the extinction equilibrium E e . However, the maturation time τ has no effects on the stability of the equilibrium of Model (1) while it could lead to stability switches for Model (3). The consequence is that Model (3) is not positive invariant and could have a oscillatory solution around the extinction equilibrium E e and the interior equilibrium E 2 .
To continue exploring how we should model population dynamics of honeybee with the proper age structure so that we could have a better understanding of important factors contributing the colony survival, we perform bifurcation diagrams on both Model (1) and Model (3).   (1) has only equilibrium dynamics either at the extinction E e or the interior equilibrium E 2 . Honeybee population data shown in Figure 6 seems to exhibit seasonality. By comparing dynamics of Model (1)

Data and Seasonality
In this section, we use the honeybee population data collected by James Harris (1980) [16] to do parameter estimations and model validations. The honeybee population data of two colonies: May 5, 1975 to Oct 22, 1975 andMay 3, 1976 to Dec 5, 1976, is shown in Figure 6 in the left (1975) and right (1976) side, respectively: The brood B (the sum of egg, larvae, and pupa) population is shown by triangle dots while the adult H population is represented by point dots. Based on Figure  6, the initial population of brood is B 0 = 6125 and adult H 0 = 5362 for 1975 while B 0 = 5982 and H 0 = 5362 for 1976. Figure 6 shows seasonality. Mathematical analysis provided in our previous section indicates that Model (3) can have oscillatory solutions under certain conditions while Model  (1) only exhibits simple equilibrium dynamics.
The questions are: (1) Is Model (3) better than Model (1) because it shows oscillatory solutions? Or (2) Is seasonality showed in honeybee population data (see Figure 6) caused by the external factors such as resource?
To address the questions above, we first assume that the queen lays egg is seasonal due to resource constraints. The literature work suggests that food, temperature, weather and oviposition place would affect the queen [3,10,24], thence her egg-laying rate is not constant, and assumed tohave the following expression: where γ indicates the length of seasonality; τ indicates the time length of the juvenile period; ψ indicates the time of the maximum laying rate; and r 0 indicates the baseline egg-laying rate.
Then we perform parameter estimations and model validations based on data shown in Figure  6: We implement the Monte Carlo parameter sweep method as our fitting method to the honeybee   population data to attain parameter estimates [7]. Essentially, we randomly sample hypotheses for the parameters following negative binomial regression with appropriate ranges (see Table.1). For each observed value, we defined the negative binomial probability density function. The mean (µ) is set by the corresponding predictive value, and variance is µ + α * µ 2 , which is α = k −1 . k indicates dispersion parameter, which we set range [0.5,2] [28,34]. Then we calculate the likelihood for negative binomial regression model to get better estimate for parameters [19,27]. Afterwards, we performed data fitting on the above model (model number) and compare the results..
We first assume that the egg laying rate r is constant. The fittings shown in Figure 7 suggest that the assumption of the queen egg laying being constant is not realistic enough. Thus we assume that the egg laying rate is a periodic function r(t) = r 0 * (1 + cos( 2π(t−ψ)  (3) is not positively invariant.
The data fittings assuming that r(t) = r 0 * (1 + cos( 2π(t−ψ) γ )) have better outcomes than the previous ones (Figure 7) by assuming r being constant. By comparison through the negative log likelihood method [19], we could deduce that Model (1) has the best fittings in both scenarios: r being constant in Figure (7) and r being periodic in Figure 8. Thus, based on both theoretical work and model validation, we could conclude that even though Model (3) could have oscillations in its solution, Model (1) with the egg laying rate r being periodic (supported by the best fitting based on data, see Figure 8) should be a better model for us to explore the important factors contributing to the healthy of honeybee colones.

Effects of seasonality:
Here we perform two scenarios that seasonality may promote (see Figure  9) or suppress (see Figure 10) the survival of honeybee colony, respectively. 365 )) whose average is r 0 = 1200, which show that honeybee colony survives.
On the other hand, Figure 10 shows that seasonability can make honeybee colony collapse. 365 )) with r 0 = 1200. In this case, we can see that seasonality may suppress the survival of honeybee colony.

Conclusions
Honeybees have dramatical decreased population over the long-term and each year [43]. As a result, great economic losses and the increase in the price of bee products have adversely affected the market [30,43]. The causes of the decline in the number of honeybees have been the great interests, whether they may directly link to human, environmental or disease [31,33,43]. Some previous work always focus on foragers or recruitment, other works investigated external causes [1, 22-24, 33, 39]. In this study, we focus on modeling proper honeybee population age structure model with model validation using empirical data to obtain better biological understanding of the critical factors that could maintain healthy of honeybee colonies.
We propose two different models with age structure to explore the importance of proper modeling. The first model (1) has an assumption that the adult bees are surviving from eggs while the second model (3) assumes that adult bees are survived from brood stage rather than the egg stage. (3) is not positively invariant. Also our bifurcation diagrams (see Figure 1 & 5) confirmed such different dynamical outcomes. Both theoretical and bifurcation results indicates that the different assumptions can lead to different age structure models with dramatic dynamical outcomes. So which model would be more appealing and biological relevant? Can we say the second model (3) be better as it has oscillatory dynamics that could be supported by seasonality observed in data?

Our theoretical work (see Theorem 3.2, 3.3, 3.5 & 3.6 ) implies that Model (1) and Model (3) have huge differences in their dynamics. Specifically, Model (1) has only equilibrium dynamics and the maturation time doesn't affect its dynamics (see Theorem 3.3 ) while Model (3) can be destabilized by the maturation time along with its life history parameter values (see Theorem 3.6) and Model
Given that the queen reproduction depends on seasonality [3,24,37], this suggests that it is paramount to include seasonality when modeling honeybee population dynamics. To address whether the seasonal pattern observed from data is due to the internal factor such as the maturation time and/or other life history parameters (for example, Model (3) could be a better model for generating seasonal patterns from the internal factors) or the external factor such as the queen egg laying rate that is regulated by the temperature and the resource [3,6,10]. We use data to validate Model (1), and Model (3) by assuming that the queen egg laying rate is constant and seasonal. Our validations on models without seasonality did not have a good fit by comparing to the corresponding models with seasonality. Among all models with seasonality, Model (1) has the best fit (see Figure 8). Our model validations with data suggest that the seasonal pattern observed from data is very likely due to the external factor such the temperature or available resources that may generate periodic dynamics in the queen egg laying rate while the internal factor such as the maturation time doesn't seem to be responsible seasonal pattern observed from data.
Both theoretical and numerical results including model validations suggest that Model (1) with seasonality in the queen egg laying rate seems to be the most fit model for studying honeybee population dynamics with age structure. Theoretical results (Theorem 3.3 ) and bifurcation diagrams ( Figure 1) imply that (1) the survival of honeybee colonies requires a large value of the queen egg-laying rate (r) and smaller values of the other life history parameter values in addition to the proper initial condition; (2) both brood and adult bee population is increasing with respect to the egg-laying rate r and is decreasing with respect to the regulation effects of brood α, the square of half maximum of colony size at which brood survival rate K, and the mortality rates d b , d h ; and (3) seasonality may promote the survival of the honeybee colony (see Figure 9) but also may lead to the colony collapsing (see Figure 10c&10d). In summary, our work suggests that Model (1) with seasonality could be used for our future model that includes more external factors, such as diseases, parasite, food and human activities [23,31,33,43]. Our ongoing work has extended the current model (1) to include parasites.

Proofs
Proof of Theorem 3.1 Proof. First, we look at the following equation that describes the population of adult bees: Since H(t) is a nonnegative continuous function during the time t ∈ [−τ, 0], the equation (15) implies that By deduction on intervals [(n − 1)τ, nτ ], n ≥ 1], we could show that H(t) ≥ 0.
By integration, we could set K+H 2 (s)+αB(s) ds. Thus, the equation (16) provides an explicit mathematical expression of B(t) which is nonnegative for all time t ≥ 0. Thus, the state space X of the proposed model (1) is positive invariant.
To show the boundedness of the model, define V = B + H, then we have

Proof of Theorem 3.2
Proof. We linearize the Model 1: For extinction equilibrium, the matrix (17) gives: and from this we obtain the following eigenvalues: Thus, we can conclude that E e is always locally asymptotically stable. (1), then (B(t), H(t)) is bounded and positive for all t > 0 by Theorem 3.1. Define lim sup t→∞ H(t) = H ∞ . By Lemma B in Appendix A, there exists sequence {t n } ↑ ∞ such that lim n→∞ H(t n ) = H ∞ , and lim n→∞ H (t n ) = 0, and same for B(t). For any > 0, there exists N such that for n > N , H(t n − τ ) ≤ H ∞ + holds. Thus, according to Model (1), for n > N we have

Now we show its global stability as follows. Let (B(t), H(t)) be a solution of Model
Let n→∞, we get It follows by the arbitrariness of that Therefore, H ∞ = 0, and hence from the positivity of solution we have lim t→∞ H(t) = 0. Furthermore, we obtain the following limiting equation which implies that lim t→∞ B(t) = 0. Therefore, E e = (0, 0) is globally asymptotically stable.

Proof of Theorem 3.3
Proof. Let (B * , H * ) be an equilibrium of Model (1). Then the linearization of the proposed model (1) at the equilibrium (B * , H * ) is shown as follows: In the case of the interior equilibrium (B * , H * ) = E i , from Equation (5), we have .
The characteristic equation (10) evaluated at the positive interior equilibrium (B * , H * ) when α = 0 gives the following expression: . This implies that the stability of the interior equilibrium (B * , H * ) is determined by the eigenvalues of the following equation evaluated at Therefore, we can obtain the following inequalities 2 ) 2 < 0. By applying Theorem 4.7 from Hal Smith [42], we can conclude that the interior equilibrium E 1 is always unstable and E 2 is always locally asymptotically stable for any delay τ > 0.

Proof of Theorem 3.4
Proof. According to Proposition 3.1, we know that Model (1) has a unique interior E = (B * , K . In order to study its stability, we define the following matrices: Let L(λ) to be represented as follows which has λ = 0 as one of its eigenvalues. By applying Lemma A in Appendix A to (19), except for the root λ = 0, all roots of (19) has negative real parts for all 0 ≤ τ < ∞. Since (19) has a simple zero eigenvalue, we need to use the center manifold and normal form theory in [13] to obtain the ,H = H − √ K, and still denoteB = B,H = H, then the system (1) becomes Let Λ = {0}. From normal form theory in [13], there exists a one-dimension ordinary differential equation which has the same dynamical property as (20) near 0. Rewriting (20) aṡ where z t (θ) = z(t + θ) ∈ C := ([−τ, 0], R 2 + ), C is the phase space with the norm |φ| = max −τ ≤θ≤0 |φ(θ)|, and where U and V are given in (18), and Then we have follows By the adjoint theory of FDE, the phase space C can be decomposed as Taking the base Ψ of adjoint space P * of P satisfies Ψ, Φ = 1, where ·, · is a bilinear function defined in C * × C by where The projection mapping π : BC→P : leads to the decomposition BC = P ⊕ Kerπ. Decomposing u in Equation (22) in the form of u = Φx + y where x ∈ R and y ∈ Q 1 := Kerπ ∩ D(A) = Q ∩ C 1 . Then (22) is equivalent to the system   ẋ = Ψ(0)F (Φx + y), Thus, the local invariable manifold of (20) at 0 with the tangency with the space P satisfies y(θ) = 0, the flow on this manifold is given by the following one-dimension ODĖ This implies that the zero solution of (23) is stable. Therefore, the interior equilibrium E = (B * , H * ) is locally asymptotically stable for all τ > 0. The proof is complete.

Proof of Theorem 3.5
Proof. We linearize the equations of Model 3: The matrix (24) evaluated at the extinction equilibrium E e gives the characteristic equation . Clearly, one characteristic root is λ = −d h < 0, others are the roots of the following equation When there is no delay, i.e., τ = 0, (25) has only a negative characteristic root λ == −d b , Model 3 is asymptotically stable at (0, 0). Moreover, for every τ ≥ 0, (25) has no nonnegative real root.
We assume λ = iw, w > 0, is a root of (25) for some τ > 0. Then, we have which gives It is clear that (27) has a positive real root if and only if e −d b τ > d b , i.e., τ < τ * := 1 d b ln 1 d b , and 0 < d b < 1. Notice that cos(wτ ) < 0, sin(wτ ) > 0, there is a unique θ, π 2 < θ < π, such that wτ = θ makes both equation of (26) hold. Then, if e −d b τ > d b , we get a set of values of τ for which there are imaginary roots: where w is given by (28) and From (25), we have Thus, Here, Clearly, In order to get the stability of the equilibrium (0, 0), we first claim: • For all τ ≥ max{0, τ * }, (0, 0) is asymptotically stable.
In fact, we rewrite (25) as the form Thus, by applying Theorem 4.7 from Hal Smith, (0, 0) is local asymptotically stable. Now, we consider the following two cases.
2 . This case is divided into two subcases: (i) d b ≥ 1 and (ii) In this subcase, the above claim implies that (0, 0) is asymptotically stable for all τ ≥ 0. This also can be proved as follows. For all τ ≥ 0, e −2d b τ − d 2 b < 0 holds, and hence the equation (27) has no positive real root. This implies that the equilibrium (0, 0) has no stability switch as τ increases in [0, ∞). Since (0, 0) is asymptotically stable at τ = 0, then it remains asymptotically stable for all τ ≥ 0.
In order to determine the local stability of the interior equilibrium E i , i = 1, 2 when τ ∈ (0, τ * ), we proceed as follows Kuang's book Chapter 3 [26].
number of times at most as τ is increased τ ∈ [0, τ * ), and eventually it becomes unstable.
This scenario is again divided into following two cases: From (53), ϕ(0) > 0 and hence for all τ ∈ [0, τ * ), c(τ ) > 0. In this case, we also have τ * ≤ 1 2d b ln , similar to the arguments above, we get that no matter the stability of E 2 can change a finite number of times at most as τ is increased in [0, τ * ).
The proof is completed. Therefore, H * 1 is monotonically increasing by d h , and H * 2 is monotonically decreasing by d h .
Finally, let us see d b . Then Therefore, H * 1 is monotonically increasing by d b , and H * 2 is monotonically decreasing by d b . 0