1. Introduction

Mosquitoes are small biting insects comprising of the family Culicidae. There are about 3,500 mosquito species in the world (grouped into 41 genera) [4,30,43]. Mosquito species, such as Anopheles, Aedes aegypti, Aedes albopictus and Culex, play significant roles as vectors of some major infectious diseases of humans, such as malaria, yellow fever, Chikungunya, west Nile virus, dengue fever, Zika virus and other arboviruses [4,48]. These diseases are transmitted from human-to-human via an effective bite from an infected female adult mosquito [4,51]. While adult male mosquitoes feed on plant liquids such as nectar, honeydew, fruit juices and other sources of sugar for energy, female mosquitoes, in addition to feeding on sugar sources (for energy), feed on the blood of human and other animals solely to acquire the proteins needed for eggs development [4,30,49,51].

Mosquitoes are the best known disease vector for vector-borne diseases (VBDs) of humans (VBDs account for 17% of the estimated global burden of all infectious diseases) [76,78]. Mosquito is the world's deadliest animal (accounting for more human deaths annually than any other animal), spreading human diseases such as malaria, dengue and yellow fever, which together are responsible for several million deaths and hundreds of millions of cases annually [13,49]. For example, malaria, transmitted by female Anopheles mosquitoes, is endemic in 91 countries, with about 40% of the world's population at risk, causing up to 500 million cases and over 1 million deaths annually [48,49,77]. Similarly, dengue, transmitted by female Aedes mosquitoes, causes over 20 million cases a year in more than 100 countries [48,77].

The life-cycle of the mosquito is completed via four distinct stages, namely: eggs, larva, pupa and adult stages, with the first three largely aquatic [48]. A female mosquito can lay about 100-300 eggs per oviposition [4,48], and this process is temperature dependent [4]. The eggs are laid at a convenient breeding site, usually a swamp or humid area in the aquatic environment (the Anopheles species typically lays their eggs on the surface of the water) and after about 2-3 days, they hatch into larva. Larvae develop through four instar stages [48,4]. At the end of each larval stage, the larvae molt, shedding their skins to allow for further growth (the larvae feed on microorganisms and organic matter in the water) [4]. During the fourth molt, the larvae mature into pupae (the whole process of maturation from larvae to pupae takes 4-10 days [51]). The pupae then develop into adult mosquitoes in about 2-4 days [4,51].

The duration of the entire life-cycle of the mosquito, from egg laying to the emergence of an adult mosquito, varies between 7 and 20 days [51], depending on the ambient temperature of the breeding site (typically a swamp or humid area) and the mosquito species involved [28] (for instance, Culex tarsalis, a common mosquito in California (USA), might go through its life cycle in 14 days at $70^0$F and take only 10 days at $80^0$F [48]). Adult mosquitoes usually mate within a few days after emerging from the pupal stage, after which they go questing for blood meal (required to produce eggs) [4]. Once a blood meal is taken, the female mosquito moves to a convenient breeding site where it lays its eggs [52]. The chances of survival of the female adult mosquitoes depend on temperature and humidity, as well as their ability to successfully obtain a blood meal while avoiding host defenses [4].

The introduction, abundance and distribution of mosquitoes worldwide have been affected by various environmental (climatic) factors such as temperature, humidity, rainfall and wind [2,15,47,54,55,57,58,60,67,79]. These factors have direct effect on different ecological aspects of the mosquito species which includes the oviposition process, development during aquatic stages and the biting rate of mosquitoes [2,47,60]. Furthermore, the oviposition process, development at aquatic stages, emergence of the adult and other developmental processes in the larval habitats of mosquitoes, play a key role in the determination of abundance and distribution of mosquitoes [3,56].

Understanding mosquito population dynamics is crucial for gaining insight into the abundance and dispersal of mosquitoes, and for the design of effective vector control strategies (that is, understanding mosquito population dynamics has important implications for the prediction and assessment of the effects of many vector control strategies [51,52]). The purpose of the current study is to qualitatively assess the impact of temperature and rainfall on the population dynamics of female mosquitoes in a certain region, and taking into consideration the dynamics of the human-vector interaction. This study extends earlier mosquito population biology in literature such as the model in [1], by designing a new temperature-and rainfall-dependent mechanistic mosquito population model that incorporates some more notable features of mosquitoes population ecology such as four stages for larval development and three different dispersal (questing for blood meal, fertilized and resting at breeding site) states of female adult mosquitoes. The non-autonomous model is formulated in Section 2 and its autonomous version is analyzed in Section 3. The full non-autonomous model is analyzed and simulated in Section 4. Since malaria is the world's most important parasitic infectious disease [34], numerical simulations of the model will be carried out using parameters relevant to the population biology of adult female Anopheles mosquitoes in Section 5 (it is worth stating that there are approximately 430 species of the Anopheles mosquitoes, of which 30-40 transmit malaria in humans (i.e., they are vectors) [4,43]).

---

2. Model formulation

This study is based on the formulation and rigorous analysis of a mechanistic model for the dynamics of female Anopheles mosquitoes in a population. The model splits the total immature mosquito population at time $t$ (denoted by $I_M(t)$) into compartments for eggs ($E(t)$), four larval stages ($L_i(t);\ i=1, 2, 3, 4$) and pupae ($P(t)$), so that $I_M(t)= E(t)+\sum\limits_{i=1}^4L_i(t)+P(t)$. Similarly, to account for the mosquito gonotrophic cycle, the population of adult female Anopheles mosquitoes at time $t$ ($A_M(t)$) is sub-divided into compartments for the class of unfertilized female vectors not questing for blood meal and fertilized female mosquitoes that have laid eggs at the mosquitoes breeding site ($V(t)$), the class of fertilized, but not producing, female mosquitoes questing for blood meal ($W(t)$), and the class of fertilized, well-nourished with blood, and reproducing female mosquitoes ($U(t)$), so that $A_M(t)= U(t)+V(t)+W(t)$. Let $N$ represents the amount of nutrients for the larvae (assumed to be constant or uniformly available at the breeding sites). The model is given by the following deterministic system of nonlinear differential equations.

$ $$\begin{split} \frac{d E}{d t}\ =& \ \psi_U(T) \left(1-\dfrac{U}{K_U}\right)U -\left[\sigma_E(R, \hat{T}) +\mu_E(\hat{T})\right] E, \\ \frac{d L_1}{d t} \ =& \ \sigma_E(R, \hat{T}) E -\left[\xi_{1}(N, R, \hat{T}) + \mu_{L}(\hat{T}) + \delta_LL \right] L_1, \\ \frac{d L_i}{d t} \ =& \ \xi_{(i-1)}(N, R, \hat{T})L_{(i-1)} -\left[\xi_{i}(N, R, \hat{T}) + \mu_L(\hat{T}) + \delta_LL\right] L_i \, ; i = 2, 3, 4, \\ \frac{d P}{d t} \ =& \ \xi_{4}(N, R, \hat{T}) L_4 -\left[\sigma_P(R, \hat{T}) + \mu_P(\hat{T})\right] P, \\ \frac{d V}{d t} \ =& \ \sigma_P(R, \hat{T})P + \gamma_U U -\frac{\eta_V H}{H+F}V -\mu_A(T)V, \\ \frac{d W}{d t} \ =& \ \frac{\eta_V H}{H+F}V -\left[\tau_W H + \mu_A(T)\right] W, \\ \frac{d U}{d t} \ =& \ \alpha\tau_W HW -\left[\gamma_U + \mu_A(T)\right] U, \\ L \ =& \ \sum\limits_{i = 1}^{4} L_i. \end{split}$$ $

(1)

In (1), $R = R(t)$, $T = T(t)$, and $\hat{T} = \hat{T}(t)$ denote rainfall, air temperature and water temperature at time $t$, respectively. Typically, a sinusoidal function, such as

$ $$\begin{align*} T(t) = T_0 \biggl[1 + T_1 \ {\rm cos} \biggl( \dfrac{2 \pi}{365} (\omega t + \theta) \biggr) \biggr], \end{align*}$$ $

(where $T_0$ is the mean annual temperature, $T_1$ captures variation about the mean, and $\omega$ and $\theta$ represent, respectively, the periodicity and phase shift of the function) is used to model local temperature fluctuations [2,24] (and similar appropriate time-dependent functions are used to account for rainfall and water temperature variability [1,55]). Thus, the functions $T(t), \hat{T}$ and $R(t)$ are assumed to be continuous, bounded, positive and $\omega$-periodic. Furthermore, the parameters $\psi_U(T), \sigma_E(R, \hat{T})$, $\sigma_P(R, \hat{T})$, $\xi_{i}(N, R, \hat{T})$ (for $i = 1, 2, 3, 4)$, $\mu_E(\hat{T})$, $\mu_L(\hat{T})$, $\mu_P(\hat{T})$ and $\mu_A(T)$ are non-negative, $\omega$-periodic, continuous and bounded functions defined on $[0, \infty)$. The term $\psi_U(T)\biggl(1-\dfrac{U}{K_U}\biggr)$ represents the density-dependent eggs oviposition rate (where $K_U > U(t)$ for all $t \geq 0$ is the environmental carrying capacity of female adult mosquitoes and $\psi_U(T)$ is the temperature-dependent egg deposition rate). Eggs hatch into the first larval stage at a rainfall-and temperature-dependent rate $\sigma_E(R, T)$. Larvae in Stage $i$ mature into Stage $i+1$ at a rate $\xi_i(N, R, T)\, (i=1, 2, 3)$, which is assumed to depend on temperature, rainfall and amount of nutrients. Larvae in Stage 4 ($L_4$ class) mature into pupae at a nutrient-, rainfall-and temperature-dependent rate $\xi_4(N, R, T)$. It should be emphasized that the maturation rates for the larval stages ($\xi_i; \ i = 1, 2, 3, 4)$ are dependent on nutrient, temperature and rainfall because, while nutrients are needed for the growth and development of the larvae, rainfall is required for availability of breeding sites and habitats and favorable temperature values improve development of the larvae [12,33,55,57,58,74]. However, extreme climate conditions such as excessive rainfall washes out the larval breeding sites, such as small stagnant water on yards or lawn, and also too hot or too cold are not favorable for the survival and maturation of the larvae [2,33,55,57,58].

Pupae mature into the female adult mosquitoes of type $V$ at a rainfall-and temperature-dependent rate $\sigma_P(R, \hat{T})$. These female adult mosquitoes quest for blood meal at the human habitat at a rate $\eta_V$ (and become female adult mosquitoes of type $W$) [51]. The term $\dfrac{H}{H + F}$ accounts for the preference of human blood, as opposed to that of other animals [28,32,52] (where $H$ is the population density of humans that are accessible to the mosquitoes (local to the breeding sites of the mosquitoes) and $F$ is a positive constant representing a constant alternative food source for the female adult mosquitoes) [51]. At the human habitat, female adult mosquitoes of type $W$ interact with humans according to standard mass action law, at a constant rate $\tau_W$ [51,52]. This interaction can be successful with probability $\alpha \in [0,1]$, so that questing mosquitoes successfully obtain blood meals and become vectors of type $U$ (at the rate $\alpha\tau_W$) which, in turn, return to become female adult mosquitoes of type $V$ at a rate $\gamma_U$ after laying eggs (see [51,52] for further details on the derivation of the $U-V-W$ component of the model). Furthermore, the parameters $\mu_E(\hat{T}), \, \mu_L(\hat{T}), \, \mu_P(\hat{T}), \, \mu_A(T)$ represent, respectively, the temperature-dependent natural death rate for female eggs, female larvae, female pupae and female adult mosquitoes, and $\delta_LL$ is the density-dependent mortality rate for larvae, accounting for intra and inter-species larval competition for resources (nutrients) and space (Abdelrazec and Gumel [1] and Lutambi et. al. [41] also incorporated density-dependent larval mortality in their models). A flow diagram of the model (1) is depicted in Figure 1, and the model variables and parameters are described in Table 1.

---

2.1. Time-dependent parameters

The functional forms of the nutrient-, rainfall- and temperature-dependent parameters of the model (1) are formulated as follows. This functional forms derived from [2,47,59,60], characterizes the female Anopheles mosquitoes (which transmits malaria in humans). The per-capita rate of deposition of female eggs ($\psi_U(T)$) defined using the quadratic function used in [47], is given by

$ $$\begin{align*} \psi_U(T) &= -0.153T^2 + 8.61T -97.7. \end{align*}$$ $

Similarly, following [47], the per-capita death rate of the female adult mosquitoes ($\mu_A(T)$) is defined as

$ $$\begin{align*} \mu_A(T) = -\ln \left( -0.000828T^2 + 0.0367T + 0.522\right). \end{align*}$$ $

Furthermore, following [60], the per-capita death rate of female eggs ($\mu_E(\hat{T})$), female larvae ($\mu_L(\hat{T})$) and female pupae ($\mu_P(\hat{T})$) are defined, respectively, as

$ $$\begin{align*} \mu_E(\hat{T}) &= \dfrac{1}{1.011 + 20.212[1 + (\frac{\hat{T}}{12.096})^{4.839}]^{-1}}, \\ \mu_L(\hat{T}) &= \frac{1}{8.130 + 13.794[1 + (\frac{\hat{T}}{12.096})^{4.839}]^{-1}}, \\% -\mu_E(\hat{T})}, \\ \mu_P(\hat{T}) &= \frac{1}{8.560 + 20.654[1 + (\frac{\hat{T}}{19.759})^{6.827}]^{-1}}. \end{align*}$$ $

Similarly, following [60], the per-capita maturation rate of eggs (into larvae) and pupae (into female adult mosquitoes) are given by

$ $$\begin{align*} \sigma_i(R, \hat{T}) = \dfrac{(1 -p_i)p_i^{d_i}}{1 -p_i^{d_i}};\ i=\{E, P\}, \end{align*}$$ $

where $d_i = d_i(\hat{T}) > 1\, (i=\{E, P\})$ is the average duration spent in Stage $i$, given by

$ $$\begin{align*} d_i(\hat{T}) = \dfrac{1}{\mu_i(\hat{T})}, \end{align*}$$ $

and $p_i = p_i(R, \hat{T})$, (where $0 \leq p_i (R, \hat{T})< 1$) is the daily survival probability of immature mosquitoes in Stage $i$ (assumed to be determined by the mean daily water temperature ($\hat{T}$ ($^\circ$C)) and cumulative daily rainfall ($R$ (mm))), so that

$ p_i(R, \hat{T}) = p_i(R)p_i(\hat{T}), $

(2)

with $p_i(\hat{T}) = \exp{(-\mu_i(\hat{T}))}$ (it should be stated that, in line with [60], the definition of $p_i(R, \hat{T}) = p_i(R)p_i(\hat{T})$ emphasizes the assumed independence of temperature and rainfall with each other). Following [59] (Supplemental Material), the rainfall-dependent daily probability of survival of immature mosquitoes, $p_i(R)$ (with $i = \{E, L_1, L_2, L_3, L_4, P\}$) is given by

$ p_{i}(R) = R(R_{I_M}-R)(4p_{Mi}/R^2_{I_M}), \, \, i = \{E, L_1, L_2, L_3, L_4, P\}, $

(3)

where $p_{Mi}$ is the peak daily survival probability of immature mosquito stage $i$ ($i = E = {\rm eggs};\, $$i = L_j = {\rm larvae\ in\ Stage}\, j; \, $$ j = \{1, 2, 3, 4\};\, $$i = P = {\rm pupae} $) and $R_{I_M} > R(t) > 0$, for all $t$, is the maximum rainfall threshold in the community. The per capita maturation rate of larvae ($\xi_j; \ j=1, 2, 3, 4$), assumed to depend on amount of available nutrients ($N$), rainfall ($R$) and water temperature ($\hat{T}$), is given by

$ $$\begin{align*} \xi_j(N, R, \hat{T}) = \xi_j(N)\xi_j(R, \hat{T});\ j = 1, 2, 3, 4, \end{align*}$$ $

where $\xi_j(N) = e_jN$, with $e_j$ representing the rate of nutrients intake for female larvae in stage $j$. Furthermore, $\xi_j(R, \hat{T})$ is given by

$ $$\begin{align*} \xi_j(R, \hat{T})= \frac{(1 -p_j)p_j^{d_j}}{1 -p_j^{d_j}}; \ j={1, 2, 3, 4}, \end{align*}$$ $

with $p_j$ having similar definition as above, $d_j = d_L(\hat{T}) = 1/\mu_L(\hat{T})$, for $ j={1, 2, 3, 4}$. It is assumed that water and air temperature obey the linear relation (see, for instance, [2,35,57,62,64]), given by $\hat{T} = T + \theta_T$, where $\theta_T > 0$ is a small increment in temperature.

Furthermore, since almost all communities within tropical and sub-tropical regions of the world record temperatures in the range $[16,40]^0$C [11], it is plausible to assume that $16^0$C $\leq T(t), \hat{T}(t) \leq 40^0$C for this study. Using this assumption, the minimum and maximum values of the temperature-dependent parameters of the model (i.e., the per capita rate of deposition of eggs ($\psi_U(T)$), per capita death rate of female eggs ($\mu_E(T)$), per capita death rate of female larvae ($\mu_L(T)$), per capita death rate of female pupae ($\mu_P(T)$) and per capita death rate of female adult mosquitoes ($\mu_A(T)$)) are tabulated in Table 2.

The non-autonomous model (1) is an extension of the autonomous model for the population biology of the mosquito developed in [51,52], by including:

(ⅰ) aquatic stages of the mosquito (i.e., adding the $E, L_1, L_2, L_3, L_4$ and $P$ classes);

(ⅱ) the effect of climate variables (i.e., adding the dependency on temperature and rainfall).

It also extends the model by Lutambi et. al. [41] by, inter alia:

(ⅰ) incorporating the effect of climate variable (temperature and rainfall);

(ⅱ) using logistic growth rate for egg oviposition rate (a constant rate was used in [41]);

(ⅲ) incorporating four larval stages (only one larval class was used in [41]).

Furthermore, the model (1) extends the non-autonomous climate-driven mosquito population biology model developed by Abdelrazec and Gumel [1] by

(ⅰ) including the dynamics of adult female mosquitoes (i.e., including compartments $U$, $V$ and $W$);

(ⅱ) including four larval classes (a single larval class is considered in [1]);

(ⅲ) including dependence on (constant and uniform availability of) nutrients for the four larval stages.

---

2.2. Basic properties

Since, as stated in Section 2, $R = R(t)$, $T = T(t)$ and $\hat{T} = \hat{T}(t)$, the temperature-and rainfall-dependent parameters of the model will, from now on, be expressed as functions of $t$ only. The basic properties of the non-autonomous model (1) will now be explored. Let, from now on, $\mu_M(t) = \min \, \{\mu_E(t), \mu_L(t), \mu_A(t)\}$. It is convenient to define $M(t) = I_M(t) + A_M(t)$, the total number of immature and matured mosquitoes at time $t$. It then follows from (1) that the rate of change of $M(t)$ is given by (where a dot represents differentiation with respect to time $t$)

$ $$\begin{split} \dot{M}\ \ \leq& \ \psi_U(t)\biggl(1 -\frac{U}{K_U}\biggr)U -\delta_LL^2 -\mu_M(t)M -(1-\alpha)\tau_WHW, \\ \leq& \ \psi_U(t)\biggl(1 -\frac{U}{K_U}\biggr)U -\mu_M(t)M, \ \ t>0. \end{split}$$ $

(4)

In order to study the asymptotic dynamics of the mosquito population, subject to fluctuations in temperature and rainfall, we assume that the mosquito population stabilizes at a periodic steady-state. Furthermore, following [40,55], it is assumed that for the time $\omega$-periodic function, $\psi_U(t) \in C^1(0, \mathbb{R}_+)$, there exists a positive number, $h_0$, such that

$ \psi_U(t)\biggl(1 -\frac{U}{K_U}\biggr)U -\mu_M(t)A < 0 \ {\rm \ for\ all}\ A \geq h_0. $

Lemma 2.1. For any $\phi \in \Omega = C([0, \mathbb{R}^9_+])$, the model (1) has a unique non-negative solution through $\phi$, and all solutions are uniformly-bounded.

Proof. Let $\phi : [0, \infty] \rightarrow \mathbb{R}^9_+$ be the vector-valued functions such that $\phi(0) = (E(0), L_1(0), L_2(0), L_3(0), L_4(0), $ $P(0), V(0), W(0), U(0)) $. The system (1) can then be re-written as:

$ \frac{d\phi}{dt} = f(t, \phi(t)), \, t\geq 0 \, \, \, \phi(0) = \phi_0, $

where,

$ $$\begin{align*} f(t, \phi(t)) = \left[\begin {array}{c} \psi_{U}(t) \biggl[1-\dfrac{\phi_9(0)}{K_U}\biggr]\phi_9(0) -\bigl[\sigma_{E}(t) +\mu_{E}(t)\bigr] \phi_1(0) \\ \sigma_{E}(t)\phi_1(0) -\bigl[\xi_{1}(t) + \mu_{L}(t) + \delta_L\phi_L(0) \bigr] \phi_2(0) \\ \xi_{(i-2)}(t)\phi_{(i-1)}(0) -\bigl[\xi_{(i-1)}(t) + \mu_{L}(t) + \delta_L\phi_L(0) \bigr] \phi_{i}(0); i = 3, 4, 5 \\ \xi_{4}(t) \phi_{5}(0) -\bigl[\sigma_{P}(t) + \mu_{P}(t)\bigr]\phi_{6}(0) \\ \sigma_{P}(t)\phi_{6}(0) + \gamma_{U} \phi_{9}(0) -\dfrac{\eta_V H}{H+F}\phi_{7}(0) -\mu_A(t)\phi_{7}(0) \\ \dfrac{\eta_V H}{H+F}\phi_{7}(0) -\bigl[\tau_W H + \mu_A(t)\bigr] \phi_{8}(0) \\ \alpha\tau_W H\phi_{8}(0) -\bigl[\gamma_U + \mu_A(t)\bigr] \phi_{9}(0) \end {array} \right], \end{align*}$$ $

$\phi_L(0) = \sum\limits_{i=2}^{5}\phi_i(0)$, with $\phi_9(0) < K_U$. Thus, for all $\phi \in \Omega$, the function $f(t, \phi(t))$ is continuous and Lipschitzian (with respect to $\phi$ in each compact set in $\mathbb{R} \times \Omega$) [40]. Hence, there is a unique solution of system (1) through $(0, \phi)$. It should be noted that $f_i(t, \pi) \geq 0$ whenever $\pi \geq 0$ and $\pi_i(0) = 0$ [40]. Hence, it follows (from Remark 5.2.1 in [69]) that $\Omega$ is positively-invariant with respect to the model (1).

For the total mosquito population $M(t)$, the rate of change of $M(t)$ satisfies Equation (5). Thus, it follows from the comparison principle [37], that the solution exists for all $t \geq 0$. Moreover,

$ $$\begin{align*} \limsup\limits_{t\rightarrow\infty}(E(t) + L(t) + P(t) + V(t) + W(t) + W(t)) \leq M^*(t), \end{align*}$$ $

where $M^*(t)$ is the unique periodic solution of

$ \dot{M}^* = \psi_U(t)\left(1 -\frac{U}{K_U}\right)U -\mu_M(t)M^*, \ \ t>0, $

(5)

given by,

$ $$\begin{align*} M^*(t) \ =& \ e^{-\int_{0}^{t}\mu_M(s)\, ds} \times \Biggr\{\int_{0}^{t}\biggl[\psi_U(s)\left(1-\frac{U(s)}{K_U}\right)U(s)e^{\int_{0}^{s}\mu_M(\tau)\, d\tau}\biggr]\, ds \\& +\ \dfrac{\int_{0}^{\omega} \psi_U(s)\biggl(1-\dfrac{U(s)}{K_U}\biggr)U(s) \exp{\bigl[\int_{0}^{s}\mu_M(\zeta)\, d\zeta\bigr]}}{\exp{\bigl[\int_{0}^{\omega}\mu_M(s)\, ds\bigr]} -1}\Biggr\}. \end{align*}$$ $

Thus, all solutions of the model (1) are ultimately-bounded [40]. Moreover, it follows from (5) that $\dot{M}^* < 0$ whenever $M^* > h_0$. Hence, all solutions of the model (1) are uniformly-bounded [40,55].

---

3. Analysis of autonomous model

It is instructive to, first of all, analyze the dynamics of the autonomous equivalent of the non-autonomous model (1) to determine whether or not it has some qualitative features that do not exist in the model (1). Consider, now, the non-autonomous model (1) with all rainfall-and temperature-dependent parameters set to a constant (i.e., $\sigma_E(t) = \sigma_E, \, \sigma_P(t) = \sigma_P, \, \xi_{i}(t) = \xi_i, \, \mu_E(t) = \mu_E, \, \mu_L(t) = \mu_L, \, \mu_A(t) = \mu_A$), given by:

$ $$\begin{align*} \frac{d E}{d t} \ =& \ \psi_U \left(1-\dfrac{U}{K_U}\right)U -\left(\sigma_E +\mu_E\right) E, \\ \frac{d L_1}{d t} \ =& \ \sigma_E E -\left[\xi_{1} + \mu_{L} + \delta_LL\right] L_1, \\ \frac{d L_i}{d t} \ =& \ \xi_{(i-1)} L_{(i-1)} -\left[\xi_{i} + \mu_L + \delta_LL \right] L_i \, ; i = 2, 3, 4, \\ \frac{d P}{d t} \ =& \ \xi_{4} L_4 -\left(\sigma_P + \mu_P\right) P, \\ \frac{d V}{d t} \ =& \ \sigma_P P + \gamma_U U -\frac{\eta_V H}{H+F}V -\mu_A V, \\ \frac{d W}{d t} \ =& \ \frac{\eta_V H}{H+F}V -\left(\tau_W H + \mu_A \right) W, \\ \frac{d U}{d t} \ =& \ \alpha\tau_W HW -\left(\gamma_U + \mu_A \right) U, \\ L \ =& \ \sum\limits_{i = 1}^{4} L_i, \end{align*}$$ $

(6)

where, now, $\xi_i = \xi_i(N)=e_iN$.

---

3.1. Asymptotic stability of trivial equilibrium point

In this section, some results for the existence and linear asymptotic stability of the trivial equilibrium point of the autonomous model (6) will be provided. It is convenient to introduce the following parameter groupings (noting that $\alpha < 1$):

$ $$\begin{split} \begin{array}{lll} \small \left\{ \begin{array}{lll} \tau^*_W = \tau_WH, \ \eta^*_V = \dfrac{\eta_VH}{H+F}, \ C_E = \sigma_E + \mu_E, \ C_P = \sigma_P + \mu_P, \\ C_i = \xi_i + \mu_L \ {\rm (for\ } i= 1, 2, 3, 4), \ C_5 = \eta^*_V + \mu_A, \ C_6 = \tau^*_W + \mu_A, \ C_7 = \gamma_U + \mu_A, \\ B = \sigma_E\sigma_P\xi_1\xi_2\xi_3\xi_4, \ C = C_1C_2C_3C_4C_EC_P, \ D = C_5C_6C_7 -\alpha\tau^*_W\eta^*_V\gamma_U > 0. \end{array} \right. \end{array} \end{split}$$ $

(7)

The autonomous model (6) has a trivial equilibrium solution, denoted by $\mathcal{T}_0$, given by

$ $$\begin{align*} \mathcal{T}_0 &= (E^*, L_1^*, L_2^*, L_3^*, L_4^*, P^*, V^*, W^*, U^*)=\bigl(0, 0, 0, 0, 0, 0, 0, 0, 0\bigr). \end{align*}$$ $

The linear stability of $\mathcal{T}_0$ (in $\Omega$) is obtained by using the next generation matrix [23,73] for the system (1). Using the notation in [73], the non-negative matrix $\mathcal{F}$ and the non-singular matrix $\mathcal{V}$, for the new egg deposition terms and the remaining transfer terms, are, respectively, given (at the trivial equilibrium, $\mathcal{T}_0$) by

$ $$\begin{align*} \mathcal{F} = \left[\begin {array}{ccc} \boldsymbol{0}&\boldsymbol{0}&\mathcal{F}_1 \\ \boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0} \\ \boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0} \end {array} \right] \, \, {\rm and \, \, } \mathcal{V} = \left[\begin {array}{ccc} \mathcal{V}_1&\boldsymbol{0}&\boldsymbol{0} \\ \mathcal{V}_2&\mathcal{V}_3&\boldsymbol{0} \\ \boldsymbol{0}&\mathcal{V}_4&\mathcal{V}_5 \end {array} \right], \end{align*}$$ $

where $\boldsymbol{0}$ denotes a zero matrix of order $3$, and

$ $$\begin{align*} \mathcal{F}_1 &= \left[\begin {array}{ccc} {0}&{0}&\psi_{{U}} \\ {0}&{0}&{0} \\ {0}&{0}&{0} \end {array} \right], \mathcal{V}_1 = \left[\begin {array}{ccc} C_E&{0}&{0} \\ -\sigma_E&C_1 &{0} \\ {0}&-\xi_1&C_2 \end {array} \right], \mathcal{V}_2 = \left[\begin {array}{ccc} {0}&{0}&-\xi_2 \\ {0}&{0}&{0} \\ {0}&{0}&{0} \end {array} \right], \\ \mathcal{V}_3 &= \left[\begin {array}{ccc} C_3 &{0}&{0} \\ -\xi_3&C_4 &{0} \\ {0}&-\xi_4&C_P \end {array} \right], \mathcal{V}_4 = \left[\begin {array}{ccc} {0}&{0}&-\sigma_P \\ {0}&{0}&{0} \\ {0}&{0}&{0} \end {array} \right], \mathcal{V}_5 = \left[\begin {array}{ccc} C_5&{0}&-\gamma_U \\ -\eta^*_V&C_6&{0} \\ {0}&-\alpha\tau^*_W&C_7 \end {array} \right]. \end{align*}$$ $

It follows from [73] that the associated vectorial reproduction number of the autonomous model (6) [63], denoted by $\mathcal{R}_0 = \rho(\mathcal{F}\mathcal{V}^{-1})$, is given by (where $\rho$ is the spectral radius of the next generation matrix $\mathcal{F}\mathcal{V}^{-1}$)

$ \mathcal{R}_0 = \dfrac{\alpha\tau^*_W\eta^*_V\psi_UB}{CD}, $

(8)

where $\tau^*_W, \, \eta^*_V, \, B, \, C, \, C_i \ (i=E, P, 1, \ldots, 7)$ and $D$ are as defined in (7). The threshold quantity, $\mathcal{R}_0$, measures the average number of new adult mosquitoes (offspring) produced by one reproductive mosquito during its entire reproductive period [52]. The result below follows from Theorem 2 in [73].

Theorem 3.1. The trivial equilibrium ($\mathcal{T}_0$) is locally-asymptotically stable (LAS) whenever $\mathcal{R}_0 < 1$, and unstable whenever $\mathcal{R}_0 > 1$.

Theorem 3.2. The trivial equilibrium point ($\mathcal{T}_0$) of the autonomous model (6) is GAS whenever $\mathcal{R}_0 \, \leq 1$.

Proof. Consider the Lyapunov function

$ $$\begin{align*} &\mathcal{K}_1 = \ \alpha\tau^*_W\eta^*_V\xi_4\sigma_P\bigl[\sigma_E\xi_1\xi_2\xi_3E + C_E\xi_1\xi_2\xi_3L_1 + C_1C_E\xi_2\xi_3L_2 + C_1C_2C_E\xi_3L_3 \\& + C_1C_2C_3C_EL_4\bigr] + C_1 C_2C_3C_4 C_E \bigl[\sigma_P\eta^*_V\alpha\tau^*_W P + C_P\eta^*_V\alpha\tau^*_W V + C_PC_5\alpha\tau^*_W W \\& + C_PC_5C_6 U\bigr]. \end{align*}$$ $

It is convenient to define

$ $$\begin{align*} S = \alpha\tau^*_W\eta^*_V\xi_4\sigma_P \bigl[C_E\xi_1\xi_2\xi_3L_1 + C_1C_E\xi_2\xi_3L_2 + C_1C_2C_E\xi_3L_3 + C_1C_2C_3C_EL_4 \bigr]. \end{align*}$$ $

Thus, the Lyapunov derivative is given by

$ $$\begin{align*} &\dot{\mathcal{K}}_1 = \ \alpha\tau^*_W\eta^*_V\xi_4\sigma_P\bigl[\sigma_E\xi_1\xi_2\xi_3\dot{E} + C_E\xi_1\xi_2\xi_3\dot{L_1} + C_1C_E\xi_2\xi_3\dot{L_2} + C_1C_2C_E\xi_3\dot{L_3} \\&+ C_1C_2C_3C_E\dot{L_4}\bigr] + C_1C_2C_3C_4C_E \bigl[\sigma_P\eta^*_V\alpha\tau^*_W \dot{P} + C_P\eta^*_V\alpha\tau^*_W\dot{V} + C_PC_5\alpha\tau^*_W \dot{W} \\&+ C_PC_5C_6 \dot{U}\bigr], \\ =&\ \alpha\tau^*_W\eta^*_VB\left[\psi_U \left(1-\frac{U}{K_U}\right)U\right]+ C_1C_2C_3C_4C_E \bigl(C_P\eta^*_V\alpha\tau^*_W\gamma_U U -C_PC_5C_6C_7 U\bigr) \\&-\delta_LLS, \\ =&\ \alpha\tau^*_W\eta^*_VB\psi_U U -CDU -\alpha\tau^*_W\eta^*_VB\psi_U\dfrac{U}{K_U}U -\delta_LLS, \\ =& \ \biggl[CD\bigl(\mathcal{R}_0-1 \bigr)-\alpha\tau^*_W\eta^*_VB\psi_U\dfrac{U}{K_U}\biggr]U - \delta_LLS, \end{align*}$$ $

where $\tau^*_W, \, \eta^*_V, \, B, \, C, \, C_i \ (i=E, P, 1, \ldots, 7)$ and $D$ are as defined in (7). Thus, it follows that, for $\mathcal{R}_0 \leq 1$ in $\Omega$, the Lyapunov derivative $\dot{\mathcal{K}}_1 < 0$. Furthermore, it follows from the LaSalle's Invariance Principle (Theorem 6.4 of [39]) that the maximal invariant set contained in $\{ (E(t), L_1(t), L_2(t), L_3(t), L_4(t), P(t), V(t), W(t), $ $ U(t) ) \in \Omega : \dot{\mathcal{K}}_1=0\}$ is the singleton $\{ \mathcal{T}_{0} \}$. Hence, $\mathcal{T}_0$ is GAS in $\Omega$ whenever $\mathcal{R}_0 \leq1$.

Theorem 3.2 shows that the mosquito population (both immature and mature) will be effectively controlled (or eliminated) if the associated vectorial reproduction threshold, $\mathcal{R}_0$, can be brought to (and maintained at) a value less than or equal to unity.

---

3.2. Existence of non-trivial equilibrium point

The existence of a non-trivial equilibrium of the model (6) will now be explored. Let $\mathcal{T}^{**}_1 = \bigl(E^{**}, L^{**}_1, L^{**}_2, $ $ L^{**}_3, L^{**}_4, P^{**}, V^{**}, W^{**}, U^{**}\bigr)$ represents an arbitrary non-trivial equilibrium of the model (6). Solving for the state variables of the model (6) at $\mathcal{T}_1^{**}$ gives

$ $$\begin{split} E^{**} =& \dfrac{\psi_U}{C_E}\biggl( 1 -\dfrac{U^{**}}{K_U}\biggr)U^{**}, E^{**} = \dfrac{1}{\sigma_E}\bigl(C_1 + \delta_L L^{**}\bigr)L^{**}_1, L^{**}_1 = \dfrac{1}{\xi_1}\bigl(C_2 + \delta_L L^{**}\bigr)L^{**}_2, \\ L^{**}_2 =& \dfrac{1}{\xi_2} \bigl(C_3 + \delta_L L^{**}\bigr)L^{**}_3, L^{**}_3 = \dfrac{1}{\xi_3}\bigl(C_4 + \delta_L L^{**}\bigr)L^{**}_4, \ L_4^{**} = \dfrac{C_PDU^{**}}{\alpha\tau^*_W\eta^*_V\sigma_P\xi_4}, \\ P^{**} =& \dfrac{DU^{**}}{\alpha\tau^{*}_W\eta^{*}_V\sigma_P}, \, V^{**}(U^{**}) = \dfrac{C_6C_7U^{**}}{\alpha\tau^{*}_W\eta^{*}_V}, \ W^{**} = \dfrac{C_7U^{**}}{\alpha\tau^*_W}, \ U^{**} = \dfrac{\alpha\tau^*_W\eta^*_V\sigma_P\xi_4L^{**}_4}{C_PD}. \end{split}$$ $

(9)

It follows from (9) that

$ L^{**}_i \, =\, \dfrac{\bigl(C_{i+1} + \delta_L L^{**}\bigr)L^{**}_{i+1}}{\xi_i}; \, i=1, 2, 3. $

(10)

Multiplying the second, third, fourth and fifth equations of (9), and substituting the expressions for $E^{**}$ and $L^{**}_4$ into the resulting equation (and simplifying), gives

$ B\alpha\tau^{*}_W\eta^{*}_V\psi_U\biggl( 1 -\dfrac{U^{**}}{K_U}\biggr) \, =\, C_EC_PD\prod\limits_{i = 1}^{4}\bigl(C_i + \bar{L}^{**}\bigr). $

(11)

Substituting the expression for $U^{**}$ in (9) into (11), and simplifying, gives (it can be shown that $L_4^{**}> 0$)

$ L_4^{**} =\dfrac{K_UC_PD}{\alpha\tau^*_W\eta^*_V\sigma_P\xi_4} \left[1-\dfrac{C_EC_PD \prod\limits_{i = 1}^{4}\bigl(C_i + \delta_L L^{**}\bigr)} {B\alpha\tau^*_W\eta^*_V\psi_U}\right]. $

(12)

Furthermore, substituting the expressions for $L_i\ (i = 1, 2, 3)$, given in (10), into $L^{**} = \sum\limits_{i=1}^4 L^{**}_i$ gives,

$ L^{**} = \dfrac{1}{\xi_1\xi_2\xi_3}\biggl[\xi_1\xi_2\xi_3 + \xi_1\xi_2\bigl( C_4 + \delta_L L^{**}\bigr) + \xi_1 \prod\limits_{i=3}^{4}\bigl( C_i + \delta_L L^{**}\bigr) + \prod\limits_{i=2}^{4}\bigl( C_i + \delta_L L^{**}\bigr) \biggr] L^{**}_4. $

(13)

Finally, substituting (12) into (13), and simplifying, shows that the non-trivial equilibria of the model (6) satisfy the following polynomial:

$ $$\begin{split} & b_7( L^{**})^7 + b_6( L^{**})^6 + b_5( L^{**})^5 + b_4( L^{**})^4 + b_3 ( L^{**})^3 + b_2 ( L^{**})^2 \\& +b_1 ( L^{**}) + b_0 = 0, \end{split}$$ $

(14)

where the coefficients $b_i\ (i = 0, \ldots, 7)$ are constants, and are given in Appendix A1. It follows from the expressions of $b_i \ (i = 0, \ldots, 7)$ in Appendix A1 that:

(ⅰ) the coefficients $b_i \ (i = 0, \ldots, 7)> 0$ whenever $\mathcal{R}_0 < 1$. Thus, no positive solution exists whenever $\mathcal{R}_0 < 1$. Furthermore, when $\mathcal{R}_0 = 1$, the coefficients $b_i \ (i = 1, \ldots, 7)> 0$ and $b_0 = 0$ (thus, the polynomial has no positive roots for the case when $\mathcal{R}_0 = 1$).

(ⅱ) the polynomial (14) has at least one positive root whenever $\mathcal{R}_0 > 1$ (using the Descartes' Rule of Signs).

These results are summarized below.

Theorem 3.3. The model (6) has at least one non-trivial equilibrium whenever $\mathcal{R}_0 > 1$, and no non-trivial equilibrium whenever $\mathcal{R}_0 \leq 1$.

Furthermore, it is worth stating that, for the special case of the autonomous model (6) with no density-dependent larval mortality (i.e., $\delta_L = 0$), the coefficients $b_i \ (i = 2, \ldots, 7) = 0$ and $b_1 = 1$. Thus, in this case, the polynomial (14) reduces to $L^{**} + b_0 = 0$, so that (where $Q_1$ and $X_6$ are defined in Appendix A)

$ L^{**} = \biggl(1 -\dfrac{1}{\mathcal{R}_0} \biggr)Q_1X_6. $

(15)

Thus, in the absence of density-dependent larval mortality (i.e., $\delta_L = 0$), the model (6) has a unique non-trivial equilibrium (denoted by $\mathcal{T}_1 = (E^{**}, L^{**}_1, L^{**}_2, L^{**}_3, L^{**}_4, $ $ P^{**}, V^{**}, W^{**}, U^{**})$) when $\mathcal{R}_0 > 1$ (the components of this equilibrium can be obtained by substituting (15) into (9)).

Theorem 3.4. The model (6) with $\delta_L = 0$ has a unique non-trivial equilibrium whenever $\mathcal{R}_0 > 1$, and no non-trivial equilibrium otherwise.

---

3.2.1. Asymptotic Stability of Non-Trivial Equilibrium Point: Special Case

Consider the special case of the autonomous model (6) in the absence of density-dependent mortality rate for larvae (i.e., $\delta_L = 0$) so that the autonomous model (6) has a unique non-trivial equilibrium ($\mathcal{T}_1$) when $\mathcal{R}_0 > 1$. Linearizing the autonomous model (6), with $\delta_L = 0$, at $\mathcal{T}_1$ gives:

$ $$\begin{align*} \mathcal{J}(\mathcal{T}_1) = \left[\begin {array}{ccccccccc}-C_E&0&0&0&0&0&0&0&\psi_U\biggl(\dfrac{2}{\mathcal{R}_0}-1\biggr) \\ \sigma_E&-C_1&0&0&0&0&0&0&0 \\ 0&\xi_1&-C_2&0&0&0&0&0&0 \\ 0&0&\xi_2&-C_3&0&0&0&0&0 \\ 0&0&0&\xi_3&-C_4&0&0&0&0 \\ 0&0&0&0&\xi_4&-C_P&0&0&0 \\ 0&0&0&0&0&\sigma_P&-C_5&0&\gamma_U \\ 0&0&0&0&0&0&\eta^*_V&-C_6&0 \\ 0&0&0&0&0&0&0& \alpha\tau^*_W &-C_7 \end {array} \right]. \end{align*}$$ $

The eigenvalues of the matrix $\mathcal{J}(\mathcal{T}_1)$ satisfy the following polynomial:

$ $$\begin{split} P_9(\lambda) =& \ \lambda^9 + A_8\lambda^8 + A_7\lambda^7 + A_6\lambda^6 + A_5\lambda^5 + A_4\lambda^4 + A_3\lambda^3 + A_2\lambda^2 + A_1\lambda \\& + CD(\mathcal{R}_0 -1), \end{split}$$ $

(16)

where $\tau^*_W, \, \eta^*_V$ and $C_i \ (i=E, P, 1, \ldots, 7)$ are as defined in (7) and $A_i \ (i = 1, \ldots, 8)$ are positive constants given in Appendix A2. It is convenient to re-write the polynomial (16) as

$ P_9(\lambda) = F(\lambda)G(\lambda) + CD\bigl(\mathcal{R}_0 -2\bigr), $

(17)

where,

$ F(\lambda) = (\lambda + C_E)(\lambda + C_P)(\lambda + C_1)(\lambda + C_2)(\lambda + C_3)(\lambda + C_4), $

(18)

and,

$ G(\lambda) = \lambda^3 + (C_5 + C_6 + C_7)\lambda^2 + (C_5C_6 + C_5C_7 + C_6C_7)\lambda + D, $

(19)

so that,

$ $$\begin{split} F(\lambda)G(\lambda) =& \ \lambda^9 + A_8\lambda^8 + A_7\lambda^7 + A_6\lambda^6 + A_5\lambda^5 + A_4\lambda^4 + A_3\lambda^3 + A_2\lambda^2 \\& + A_1\lambda + CD. \end{split}$$ $

(20)

The asymptotic stability of $\mathcal{T}_1$ will be explored using the properties of Bézout matrices [31]. Consequently, it is convenient to recall the following four results:

Theorem 3.5.(Routh-Hurwitz)[31]. Let $A$ be an $n \times n$ complex matrix, and let $E_k$ be the sum of all principal minors of $A$ of order $k$, $k \in\, <n>$. Let $\Omega(A)$ be the $n \times n$ Hurwitz matrix of $A$ and assume that $\Omega(A)$ is real. Then $A$ is stable if and only if all leading principal minors of $\Omega(A)$ are positive.

Definition 3.6.(Bézout Matrix)[31]. Let $a(x)$ and $b(x)$ be two polynomials with real coefficients of degree $n$ and $m$ respectively, $n \geq m$. The Bézoutiant defined by $a(x)$ and $b(x)$ is the bilinear form

$ $$\begin{align*} \dfrac{a(x)b(y)-a(y)b(x)}{x-y}\, =\, \sum\limits_{i, k = 0}^{n-1}b_{ik}x^iy^k. \end{align*}$$ $

The symmetric matrix $(b_{ik})_0^{n-1}$ associated with this bilinear form is called the Bézout matrix and is denoted by $B_{a, b}$. Each entry $b_{i, j}$ of $B_{a, b}$ can be computed separately by the entry formula

$ $$\begin{align*} b_{i, j} = \sum\limits_{k=\max (0, i-j)}^{\min (i, n-1-j)} \bigl(b_{i-k}a_{j+1+k} -a_{i-k}b_{j+1+k}\bigr) \ {\rm for\ all\ } i, j\, \leq n. \end{align*}$$ $

Theorem 3.7.(Liénard-Chipart)[31] Let $f(x)= x^n -a_n x^{n-1} -\ldots -a_1$ be a polynomial with real coefficients, and let $a_{n-1} = -1$. Define the polynomials

$ $$\begin{align*}h(u) = -a_1 -a_3u -\cdots, \\ g(u) = -a_2 -a_4u -\cdots. \end{align*}$$ $

The polynomial $f(x)$ is negative stable if and only if the Bézout matrix $B_{h, g}$ is positive definite and $a_i < 0$ for all $i \in \, < n>$.

Theorem 3.8. (Sylvester's Criterion)[27] A real, symmetric matrix is positive definite if and only if all its principal minors are positive.

We claim the following result.

Lemma 3.9. The polynomial $F(\lambda)G(\lambda)$, defined by Equations (18), (19) and (20), is Hurwitz stable (i.e., all its roots have negative real part).

Proof. It follows from the equation for $F(\lambda)$ in (18) that all roots of $F(\lambda)$ are negative. Furthermore, consider $G(\lambda) = 0$ from (19). That is,

$ $$\begin{align*} G(\lambda) = \lambda^3 + (C_5 + C_6 + C_7)\lambda^2 + (C_5C_6 + C_5C_7 + C_6C_7)\lambda + D = 0. \end{align*}$$ $

Using the Routh-Hurwitz Criterion (Theorem 3.5), the principal minors, $\Delta_k\ (k=1, 2, 3), $ of the associated Hurwitz matrix for $G(\lambda)$ are

$ $$\begin{align*} \Delta_1 =&\ C_5 + C_6 + C_7 > 0, \\ \Delta_2 =&\ (C_5 + C_6 + C_7)(C_5C_6 + C_5C_7) + C_6C_7(C_6 + C_7) + \alpha\tau^*_W\eta^*_V\gamma_U > 0, \\ \Delta_3 =&\ D\Delta_2 > 0. \end{align*}$$ $

Thus, all the roots of $G(\lambda)$ have negative real part. Hence, all nine roots of $F(\lambda)G(\lambda)$ have negative real part.

Remark 1. It follows from Lemma 3.9 and Theorem 3.7 that the corresponding Bézout matrix of $F(\lambda)G(\lambda)$ is positive-definite [31].

Remark 2. Consider $P_9(\lambda) = F(\lambda)G(\lambda) + CD(\mathcal{R}_0 -2)$. Then, $P_9(\lambda) \leq F(\lambda)G(\lambda)$ whenever $1 < \mathcal{R}_0 \leq 2$. Thus, it follows from Lemma 3.9 that all nine roots of $P_9(\lambda)$ have negative real part whenever $1 < \mathcal{R}_0 \leq 2$ (hence, $\mathcal{T}_1$ is LAS whenever $1 <\mathcal{R}_0 \leq 2$).

Furthermore, consider the characteristic polynomial $P_9(\lambda)$ given in (16). Let $A_0 = CD(\mathcal{R}_0 -1)$, with $C$ and $D$ as defined in (7), so that

$ $$\begin{align*} P_9(\lambda) = \lambda^9 + A_8\lambda^8 + A_7\lambda^7 + A_6\lambda^6 + A_5\lambda^5 + A_4\lambda^4 + A_3\lambda^3 + A_2\lambda^2 + A_1\lambda + A_0. \end{align*}$$ $

To apply Theorem 3.7, let

$ $$\begin{align*} h(u) = A_0 + A_2u + A_4u^2 + A_6u^3 + A_8u^4, \end{align*}$$ $

and,

$ $$\begin{align*}g(u) = A_1 + A_3u + A_5u^2 + A_7u^3 + u^4. \end{align*}$$ $

Thus, it follows from Definition 3.6 that the corresponding Bézout matrix of $P_9(\lambda)$, denoted by $B_{h, g}(P_9)$, is given by

$ $$\begin{align*} \tiny{ B_{h, g}(P_9) = \left[\begin {array}{cccc} A_1A_2-A_0A_3&A_1A_4-A_0A_5&A_1A_6-A_0A_7&A_1A_8 -A_0 \\ A_1A_4 -A_0A_5&A_3A_4 -A_2A_5 + A_1A_6 -A_0A_7&A_3A_6 -A_2A_7 + A_1A_8 -A_0&A_3A_8 -A_2 \\ A_1A_6 -A_0A_7&A_3A_6 -A_2A_7 + A_1A_8 -A_0&A_5A_6 -A_4A_7 + A_1A_8 -A_2&A_5A_8 -A_4 \\ A_1A_8 -A_0&A_3A_8 -A_2&A_5A_8 -A_4& A_7A_8 -A_6 \end {array} \right].} \end{align*}$$ $

Sylvester's Criterion (Theorem 3.8) can be used to obtain the necessary and sufficient conditions for $B_{h, g}(P_9)$ to be positive-definite. First of all, it is evident that $B_{h, g}(P_9)$ is symmetric. It then suffices to show that the $k^{th}$ leading principal minor of $B_{h, g}(P_9)$ is positive (i.e., to show that the determinant of the upper-left $k \times k$ sub-matrix of $B_{h, g}(P_9)$ is positive). It is convenient to introduce the following notations:

(ⅰ) $b_{i, j}^{(J)}:0\leq i, j \leq 3, J \in \{FG, \, P_9 \}$ are the entries of the corresponding Bézout matrix of the polynomial $F(\lambda)G(\lambda)$ (denoted by $B_{h, g}(FG))$ and $P_{9}(\lambda)$ (clearly, $B_{h, g}(FG) = B_{h, g}(P_9)$ when $A_0 = CD$).

(ⅱ) $\Delta_k^{(P_9)}$ is the $k^{th}$ leading principal minor of Bézout matrix $B_{h, g}(P_9)$.

Therefore, $B_{h, g}(P_9)$ can be re-written (in terms of the entries of the positive-definite Bézout matrix, $B_{h, g}(FG)$. That is, in terms of $b_{i, j}^{(FG)}$) as

$ $$\begin{align*} B_{h, g}(P_9) = \left[\begin {array}{cccc} b^{(FG)}_{0, 0}-CDKA_3&b^{(FG)}_{0, 1}-CDKA_5&b^{(FG)}_{0, 2}-CDKA_7&b^{(FG)}_{0, 3} -CDK \\ b^{(FG)}_{1, 0} -CDKA_5&b^{(FG)}_{1, 1} -CDKA_7&b^{(FG)}_{1, 2} -CDK&b^{(FG)}_{1, 3} \\ b^{(FG)}_{2, 0} -CDKA_7&b^{(FG)}_{2, 1} -CDK&b^{(FG)}_{2, 2}&b^{(FG)}_{2, 3} \\ b^{(FG)}_{3, 0} -CDK&b^{(FG)}_{3, 1}&b^{(FG)}_{3, 2}& b^{(FG)}_{3, 3} \end {array} \right]. \end{align*}$$ $

where $C$ and $D$ are as defined in (7) and $K = (\mathcal{R}_0 -2).$ It follows from Remark 2 that the Bézout matrix, $B_{h, g}(P_9)$, is a positive definite matrix for $1 < \mathcal{R}_0 \leq 2$. Furthermore, the Bézout matrix, $B_{h, g}(P_9)$, can be re-written as (after row-column operations)

$ B_{h, g}(P_9) = \left[\begin {array}{cccc} \Delta^{(P_9)}_1&b^{(P_9)}_{0, 1}&b^{(P_9)}_{0, 2}&b^{(P_9)}_{0, 3} \\ 0&\dfrac{\Delta^{(P_9)}_2}{\Delta^{(P_9)}_1}&b^{(P_9)}_{1, 2}-\dfrac{b^{(P_9)}_{0, 1}b^{(P_9)}_{0, 3}}{b^{(P_9)}_{0, 0}}&b^{(P_9)}_{1, 3}-\dfrac{b^{(P_9)}_{0, 1}b^{(P_9)}_{0, 4}}{b^{(P_9)}_{0, 0}} \\ 0&0&\dfrac{\Delta^{(P_9)}_{3}}{\Delta^{(P_9)}_{2}}&B_1 \\ 0&0&0& \dfrac{\Delta^{(P_9)}_{4}}{\Delta^{(P_9)}_{3}} \end {array} \right], $

(21)

where $B_1 = B_1\bigl(b^{(P_9)}_{i, j}\bigr) \ 0 \leq i, j\leq 3$. However, since the $k^{th}$ leading principal minor of a triangular matrix is the product of its diagonal elements up to row $k$, Sylvester's Criterion is equivalent to finding conditions for which all the diagonal elements of Bézout matrix, $B_{h, g}(P_9)$, in (21) are all positive (i.e., finding the conditions for which $\Delta_{k}^{(P_9)}$ are positive for all $k = 1, 2, 3, 4$) [63]. For example, it can be verified that the first leading principal minor of the matrix $B_{h, g}(P_9)$, given by

$ $$\begin{align*} \Delta^{(P_9)}_{1} =&\ A_1A_2 -CD(\mathcal{R}_0 -1)A_3 = A_1A_2 -CDA_3 -CD(\mathcal{R}_0 -2)A_3 \\ =& \ b^{(FG)}_{0, 0} -CD(\mathcal{R}_0 -2)A_3, \end{align*}$$ $

is positive whenever the following inequality holds:

$ $$\begin{align*} \mathcal{R}_0 < 2 + \dfrac{b^{(FG)}_{0, 0}}{CDA_3} = 2 + Z_1, \end{align*}$$ $

where, $b^{(FG)}_{0, 0} = A_1A_2 -CDA_3 > 0$ is the first leading principal minor of $B_{h, g}(FG)$ and $Z_1$ is the positive constant such that $\Delta^{(P_9)}_{1}$ is positive whenever $\mathcal{R}_0 < 2 + Z_1$. Similarly, we obtain constants $Z_k = Z_k\bigl(b^{(FG)}_{i, j}\, (0 \leq i, j\leq 3), \, A_i\, (1 \leq i\leq 9), C, D\bigr)$ such that the $k^{th}$ principal minor $\Delta^{(P_9)}_{k}$ is positive whenever $\mathcal{R}_0 < 2 + Z_k$, for each $k = 2, 3, 4$ (i.e., $Z_i$ ($i = 1, 2, 3, 4$), is the constant such that $\mathcal{R}_0 < 2 + Z_i$ makes the determinant of the associated matrix of minors of matrix (21) to be positive). Therefore, the result below follows (from the above derivations and Remark 2).

Theorem 3.10 Consider the model (6) with $\delta_L = 0$. The unique non-trivial equilibrium ($\mathcal{T}_1$) is LAS in $\Omega\setminus\{\mathcal{T}_0\}$ whenever

$ $$\begin{align*} 1 < \mathcal{R}_0 < \mathcal{R}^C_0 = 2 + \min \left\{ Z_k: \Delta^{(P_9)}_{k} > 0\ {\rm for \ all\ } k = 1, 2, 3, 4 \right\}, \end{align*}$$ $

and unstable whenever $\mathcal{R}_0 > \mathcal{R}^C_0$.

The results above (Theorem 3.4 and Theorem 3.10) show that the condition $\mathcal{R}_0 > 1$ defines the existence of a unique non-trivial equilibrium $\bigl( \mathcal{T}_1 \bigr)$ of the model (6) with $\delta_L = 0$ (which is LAS if $1 < \mathcal{R}_0 < \mathcal{R}^C_0 $). Thus, it can be deduced that, to maintain a non-trivial mosquito population, each reproducing female mosquito (of type $U$) must produce at least one egg during its entire reproductive life period (see also [51]). In other words, an increase in female mosquitoes of type $U(t)$ leads to a corresponding increase in the number of female eggs laid in the population ($E(t)$. We claim the following result.

Theorem 3.11. Consider a special case of the model (6) with $sgn(E^{**} -E(t)) = sgn(U^{**} -U(t))$ for all $t \geq 0$ and $\delta_L = 0$. Then, the non-trivial equilibrium ($\mathcal{T}_1$) of the model (6) with $\delta_L = 0$ is GAS in $\Omega\setminus\{\mathcal{T}_0\}$ whenever $1 < \mathcal{R}_0 < \mathcal{R}^C_0$.

Proof. The proof of Theorem 3.11, based on using a non-linear Lyapunov function of Goh-Voltera type, is given in Appendix B.

The ecological implication of Theorem 3.11 is that mosquitoes will persist in the community whenever the associated conditions for the global asymptotic stability of the non-trivial equilibrium ($\mathcal{T}_1$) are satisfied. The results of Theorem 3.11 are illustrated numerically, by simulating the compartment of adult mosquitoes of type $U$ in autonomous model (6) with $\delta_L = 0$ using appropriate parameter values (Figure 2). These simulation results show convergence of the solution of $U(t)$ to $U^{**}$ (in line with Theorem 3.11). It is worth mentioning that, for the fixed values of the parameters used in Figure 2, the associated bifurcation point of the model (6) with $\delta_L = 0$ is $\psi_U = \psi^*_U = 107.889493160695073$ (so that, $\Delta_4 = 0$). This is equivalent to $\mathcal{R}_0 = \mathcal{R}^C_0 = 4.5573$. Therefore, for this particular set of parameter values, the non-trivial equilibrium ($\mathcal{T}_1$) is LAS for $1 < \mathcal{R}_0 < 4.5573$, and unstable whenever $\mathcal{R}_0 > 4.5573$.

---

3.3. Hopf Bifurcation analysis: Special case

Consider the model (6) with $\delta_L = 0$ and $\mathcal{R}_0 > 1$ (so that $\mathcal{T}_1$ exists, by Theorem 3.4). Hopf bifurcation can occur (at a fixed value of a chosen bifurcation parameter) when the Jacobian of the system (6) with $\delta_L = 0$, evaluated at $\mathcal{T}_1$, has a pair of purely imaginary eigenvalues. That is, when $P_9(\lambda)$ given in (16), has a pair of purely imaginary roots.

The rank and signature of the Bézout matrix, $B_{h, g}(P_9)$, can be used to evaluate the number of roots with negative real parts [61]. The direct effect of the characteristic polynomial $P_9$ having a pair of purely imaginary eigenvalues is that the rank of the Bézout matrix, $B_{h, g}(P_9)$, is reduced by exactly one [61]. From the stability point of view, this possibility represent the existence of a boundary (Hopf bifurcation) [61]. To prove the existence of Hopf bifurcation, it also suffices to verify the transversality condition [20].

Theorem 3.12. Consider the autonomous model (6) with $\delta_L = 0$. A Hopf bifurcation occurs whenever $\mathcal{R}_0 = \mathcal{R}^C_0 = 2 + Z_4$ or, equivalently, whenever

$ $$\begin{align*} \psi_U = \psi^*_U = \dfrac{CD(2 + Z_4)}{\alpha\tau_W^*\eta_V^*B} = \dfrac{CD\mathcal{R}^C_0}{\alpha\tau_W^*\eta_V^*B}, \end{align*}$$ $

where $Z_4$ is as defined in Theorem 3.10.

Proof. To prove Theorem 3.12, it is sufficient to establish the transversality condition [20]. Let $\psi_U = \psi^*_U$ be a bifurcation parameter (and all other parameters of the model (6) are fixed). Then $\mathcal{R}_0 = \mathcal{R}^C_0 = 2 + Z_4$. Since $Z_4 < Z_i$ ($\, i = 1, 2, 3$), it follows from Theorem 3.10 that $\Delta^{(P_9)}_i > 0$ for all $\, i = 1, 2, 3$. Furthermore, $\Delta^{(P_9)}_4$ can be re-written as

$ $$\begin{align*} \Delta^{(P_9)}_4 = \left| \begin {array}{cccc} b^{(FG)}_{0, 0} -(\alpha\tau^*_W\eta_V^*\psi_UB -2CD)A_3&b^{(FG)}_{0, 1} -(\alpha\tau^*_W\eta_V^*\psi_UB -2CD)A_5&b^{(P_9)}_{0, 2}&b^{(P_9)}_{0, 3} \\ b^{(FG)}_{1, 0} -(\alpha\tau^*_W\eta_V^*\psi_UB -2CD)A_5&b^{(FG)}_{1, 1} -(\alpha\tau^*_W\eta_V^*\psi_UB -2CD)A_7&b^{(P_9)}_{1, 2}&b^{(P_9)}_{1, 3} \\ b^{(FG)}_{2, 0} -(\alpha\tau^*_W\eta_V^*\psi_UB -2CD)A_7&b^{(FG)}_{2, 1} -(\alpha\tau^*_W\eta_V^*\psi_UB -2CD)&b^{(P_9)}_{2, 2}&b^{(P_9)}_{2, 3} \\ b^{(FG)}_{3, 0} -(\alpha\tau^*_W\eta_V^*\psi_UB -2CD)&b^{(P_9)}_{3, 1}&b^{(P_9)}_{3, 2}& b^{(P_9)}_{3, 3} \end {array} \right|. \end{align*}$$ $

Hence, $\Delta^{(P_9)}_4(\psi_U) = 0$ if and only if $\psi_U = \psi^*_U$. Furthermore, it can be verified that

$ $$\begin{align*} \dfrac{d \Delta^{(P_9)}_4(\psi_U)}{d\psi_U}\biggr|_{\psi_U = \psi^*_U} =\ {\rm Tr } \left( {\rm Adj} (B_{h, g}(P_9)(\psi_U))\biggr|_{\psi_U = \psi^*_U} \dfrac{dB_{h, g}(P_9)(\psi_U)}{d\psi}\biggr|_{\psi_U = \psi^*_U}\right) \neq 0, \end{align*}$$ $

where 'Tr' and 'Adj' denote, respectively, the trace and adjoint of a matrix. Similarly, let $\mu_A$ be a bifurcation parameter (and all other parameters of the model (6) are fixed). Thus,

$ $$\begin{align*} \dfrac{d \Delta^{(P_9)}_4(\mu_A)}{d\psi}\biggr|_{\mu_A = \mu_A^*} =\ {\rm Tr } \left( {\rm Adj} (B_{h, g}(P_9)(\mu_A))\biggr|_{\mu_A = \mu_A^*} \dfrac{dB_{h, g}(P_9)(\mu_A)}{d\mu_A}\biggr|_{\mu_A = \mu_A^*}\right), \end{align*}$$ $

for all $ \Delta^{(P_9)}_4(\mu^*_A) = 0$. It can be verified that $\dfrac{d \Delta^{(P_9)}_4(\mu_A)}{d\psi}\biggr|_{\mu_A = \mu_A^*} \neq 0$.

Theorem 3.12 shows that sustained oscillations are possible, with respect to the autonomous model (6) with $\delta_L = 0$, whenever $\mathcal{R}_0 = \mathcal{R}^C_0$. This result, which is numerically illustrated in Figure 3(a), is in line with that reported in [1] (where both logistic and Maynard-Smith-Slatkin eggs-laying functions are used). It is worth mentioning that, in the proof of Theorem 3.12, two bifurcation parameters ($\psi_U$ and $\mu_A$) were considered. The reason is, as noted in [1], that the transversality condition may fail at some points if only one parameter is used (see also [20]). The nature of the Hopf bifurcation property of the model (6) is investigated numerically. The results obtained, depicted in Figure 3(b), show convergence of the solutions to a stable limit cycle. It should be mentioned that the presence of bifurcation was not shown in the dynamics of mosquito population biology model developed in [41].

---

3.3.1. Numerical illustrations

In this section, a bifurcation diagram for the autonomous model (6) with $\delta_L = 0$, which summarizes the main results obtain in Section 3, will be generated in the $\mu_A -\psi_U$ plane as follows:

(ⅰ) Solving for $\psi_U$ from $\mathcal{R}_0 = 1$ gives the following equation for $\psi^l_U$ (depicted in Figure 4):

$ $$\begin{align*} l: \ \ \psi_U = \psi^l_U = \dfrac{CD(\mu_A)}{\alpha\tau^*_W\eta^*_VB}. \end{align*}$$ $

(ⅱ) Solving for $\psi_U$ from $\Delta^{(P_9)}_4 = 0$ (and fixing all parameters of the models (using their values as in Figure 4), except the parameters, $\mu_A$ and $\psi_U$) give the following curve $\Delta^{(P_9)}_4 = 0$:

$ $$\begin{align*} \mathcal{H}:\ \ \psi_U = \psi^*_U = \dfrac{CD(\mu_A)\bigl[2 + Z_4(\mu_A)\bigr]}{\alpha\tau_W^*\eta_V^*B}, \end{align*}$$ $

where $B, C$ and $D$ are as defined in (7) and $Z$. The curves $l$ and $\mathcal{H}$ (depicted in Figure 4) divide the $\mu_A -\psi_U$ plane into three distinct regions, namely $\mathcal{D}_1$, $\mathcal{D}_2$ and $\mathcal{D}_3$, given by:

$ $$\begin{align*} \mathcal{D}_1 \ =&\, \bigl\{ (\mu_A, \psi_U): \, 0 < \psi_U \leq \psi^l_U;\, \mu_A > 0\bigr\}, \\ \mathcal{D}_2 \ =& \, \bigl\{ (\mu_A, \psi_U): \, \psi^l_U < \psi_U < \psi^*_U;\, \mu_A > 0\bigr\}, \\ \mathcal{D}_3 \ =& \, \bigl\{ (\mu_A, \psi_U): \, \psi_U > \psi^{*}_U;\, \mu_A > 0\bigr\}. \end{align*}$$ $

The regions can be described as follows (see also Table 3):

(ⅰ) Region $\mathcal{D}_1$: In this region, $\mathcal{R}_0 \leq 1$. Hence, in this region (note that $\delta_L = 0$), the trivial equilibrium ($\mathcal{T}_0$) is globally-asymptotically stable (in line with Theorem 3.2).

(ⅱ) Region $\mathcal{D}_2$: Here, $1 < \mathcal{R}_0 < \mathcal{R}^C_0$. Thus, the model has two equilibria, namely the unstable trivial equilibrium ($\mathcal{T}_0$) and the locally-asymptotically stable non-trivial equilibrium ($\mathcal{T}_1$). The model undergoes a Hopf bifurcation whenever $\mathcal{R}_0 = \mathcal{R}^C_0$.

(ⅲ) Region $\mathcal{D}_3$: In this region, $\mathcal{R}_0 > \mathcal{R}^C_0$. Thus, the model has the unstable trivial equilibrium ($\mathcal{T}_1$), unstable non-trivial equilibrium and a stable limit cycle.

---

3.4. Uncertainty and sensitivity analysis for autonomous model

Sensitivity analysis determines the effects of parameters on the model outcomes [16]. The effect of these uncertainties, as well as the determination of the parameters that have the greatest influence on the mosquitoes dispersal dynamics (with respect to a given response function), are carried out using an uncertainty and sensitivity analysis [2,14,44,45,46,55]. In particular, following [14], the Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCC) will be used for the autonomous model (6). The range and baseline values of the parameters, tabulated in Table 4, will be used. Appropriate response functions are chosen for these analyses.

Using the population of female adult mosquitoes of type $U$ as the response function, it is shown in Table 5 that the top three PRCC-ranked parameters of the autonomous version of the model are the probability of female adult mosquito of type $W$ successfully taking a blood meal ($\alpha$), the natural mortality rate of female adult mosquitoes ($\mu_A$) and the natural mortality rate of female larvae ($\mu_L$). Similarly, using the population of female adult mosquitoes of type $V$ as the response function, the top three PRCC-ranked parameters are the natural mortality rate of female larvae ($\mu_L$), the deposition rate of female eggs ($\psi_U$) and the maturation rate of female larvae from Stage 1 to Stage 2 ($\xi_1$). Furthermore, using the population of female larvae in stage 4 ($L_4$) and population of female pupae ($P$) as the response functions, it is shown that the same top three PRCC-ranked parameters appeared as in the case when the population of female adult mosquitoes of type $V$ is chosen as the response function for both cases. However, using the vectorial reproduction number of the autonomous version of the model ($\mathcal{R}_0$) as the response function, the top three PRCC-ranked parameters are the natural mortality rate of female larvae ($\mu_L$), the deposition rate of female eggs ($\psi_U$) and the natural mortality rate of female adult mosquitoes ($\mu_A$).

In summary, this study identifies five parameters that dominate the population dynamics and dispersal of the mosquito, namely the probability of female adult mosquito of type $W$ successfully taking a blood meal ($\alpha$), the natural mortality rate of female adult mosquitoes ($\mu_A$), the natural mortality rate of female larvae ($\mu_L$), the deposition rate of female eggs ($\psi_U$) and the maturation rate of female larvae ($\xi_i$). The effect of the aforementioned most dominant parameters ($\alpha, \psi_U, \xi_i, \mu_L$ and $\mu_A$) on the population dynamics of the mosquito and the reproduction threshold ($\mathcal{R}_0$) is tabulated in Table 6.

---

4. Analysis of non-autonomous model

In this section, dynamical properties of the non-autonomous model (1) will be explored. The non-autonomous model (1) has a unique trivial equilibrium point denoted by $\mathcal{T}^*_0 = \bigl(E^{*}, L^{*}_1, L^{*}_2, L^{*}_3, L^{*}_4, P^{*}, V^{*}, $ $W^{*}, U^{*}\bigr) = (0, 0, 0, 0, 0, 0, 0, 0, 0)$ and positive periodic solution(s) denoted by $\mathcal{T}^*_1 = (E^{*}(t), L^{*}_1(t), L^{*}_2(t), L^{*}_3, L^{*}_4(t), $ $P^{*}(t), V^{*}(t), W^{*}(t), U^{*}(t)) $ which satisfies the unique periodic system:

$ $$\begin{align*} \frac{d E^{*}(t)}{d t} \ =& \ \psi_U(t) \left[1-\dfrac{U^{*}(t)}{K_U}\right]U^*(t) -\left[\sigma_E(t) +\mu_E(t)\right] E^{*}(t), \\ \frac{d L_1^{*}(t)}{d t} \ =& \ \sigma_E(t) E^{*}(t) -\left[\xi_{1}(t) + \mu_{L}(t) + \delta_L L^{*}(t) \right] L^{*}_1(t), \\ \frac{d L_i^{*}(t)}{d t} \ =& \ \xi_{(i-1)}(t)L^{*}_{(i-1)}(t) -\left[\xi_{i}(t) + \mu_L(t) + \delta_L L^{*}(t) \right] L^{*}_i(t) \, ;\ i = 2, 3, 4, \\ \frac{d P^{*}(t)}{d t} \ =& \ \xi_{4}(t) L_4^{*}(t) -\left[\sigma_P(t) + \mu_P(t)\right] P^{*}(t), \\ \frac{d V^{*}(t)}{d t} \ =& \ \sigma_P(t)P^{*}(t) + \gamma_U U^{*}(t) -\frac{\eta_V H}{H+F}V^{*}(t) -\mu_A(t)V^{*}(t), \\ \end{align*}$$ $

$ $$\begin{split} \frac{d W^{*}(t)}{d t} \ =& \ \frac{\eta_V H}{H+F}V^{*}(t) -\left[\tau_W H + \mu_A(t)\right] W^{*}(t), \\ \frac{d U^{*}(t)}{d t} \ =& \ \alpha\tau_W HW^{*}(t) -\left[\gamma_U + \mu_A(t)\right] U^{*}(t), \\ L^{*}(t) \ =& \ \sum\limits_{i = 1}^{4} L^{*}_i(t). \end{split}$$ $

(22)

---

4.1. Computation of vectorial reproduction ratio

The vectorial reproduction ratio, associated with the non-autonomous model (6), will be computed using the approach in [5,6,7,8,9,10,75]. The next generation matrix $F(t)$ (of the new eggs deposited) and the $M$-Matrix $V(t)$ (of the remaining transfer terms), associated with the non-autonomous model (6) (linearized at the trivial equilibrium $\mathcal{T}^*_0$), are given, respectively, by

$ $$\begin{align*} F(t) = \left[\begin {array}{ccc} \boldsymbol{0}&\boldsymbol{0}&F_1(t) \\ \boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0} \\ \boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0} \end {array} \right] \, \, {\rm and \, \, } V(t) = \left[\begin {array}{ccc} V_1(t)&\boldsymbol{0}&\boldsymbol{0} \\ V_2(t)&V_3(t)&\boldsymbol{0} \\ \boldsymbol{0}&V_4(t)&V_5(t) \end {array} \right], \end{align*}$$ $

where $\boldsymbol{0}$ denotes a zero matrix of order $3$, and

$ $$\begin{align*} F_1(t) =& \left[\begin {array}{ccc} {0}&{0}&\psi_U(t) \\ {0}&{0}&{0} \\ {0}&{0}&{0} \end {array} \right], V_1(t) = \left[\begin {array}{ccc} C_E(t)&{0}&{0} \\ -\sigma_E(t)&C_1(t) &{0} \\ {0}&-\xi_1(t)&C_2(t) \end {array} \right], \\ V_2(t) = &\left[\begin {array}{ccc} {0}&{0}&-\xi_2(t) \\ {0}&{0}&{0} \\ {0}&{0}&{0} \end {array} \right], V_3(t) = \left[\begin {array}{ccc} C_3(t) &{0}&{0} \\ -\xi_3(t)&C_4(t) &{0} \\ {0}&-\xi_4(t)&C_P(t) \end {array} \right], \\ V_4(t) = &\left[\begin {array}{ccc} {0}&{0}&-\sigma_P(t) \\ {0}&{0}&{0} \\ {0}&{0}&{0} \end {array} \right], V_5 = \left[\begin {array}{ccc} C_5(t)&{0}&-\gamma_U \\ -\eta^*_V&C_6(t)&{0} \\ {0}&-\alpha\tau^*_W&C_7(t) \end {array} \right], \end{align*}$$ $

where $\tau^*_W, \, \eta^*_V, \, C_i(t) \ (i=E, P, 1, \ldots, 7)$ are as defined in (7). The linearized version of the model (1), at $\mathcal{T}^*_0$, can be expressed as

$ $$\begin{align*} \dfrac{dx(t)}{dt} = \bigl[F(t)-V(t)\bigr]x(t) \end{align*}$$ $

where $x(t) = \bigl(E(t), L_1(t), L_2(t), L_3(t), L_4(t), P(t), V(t), W(t), U(t) \bigr)$. Following [75], let $Y(t, s)$, $t\geq s$, be the evolution operator of the linear $\omega$-periodic system $\dfrac{d y}{d t} = -V(t)y$. Thus, for each $s \in \mathbb{R}$, the associated $9 \times 9$ matrix $Y(t, s)$ satisfies [75]

$ $$\begin{align*} \frac{d Y(t, s)}{d t} = -V(t)Y(t, s)\, \, \forall \, \, t\geq s, \, \, Y(s, s) = I. \end{align*}$$ $

where $I$ is the $9 \times 9$ identity matrix.

Suppose that $\phi(s)$ ($\omega$-periodic in $s$) is the initial distribution of new eggs. Thus, $F(s)\phi(s)$ is the rate of generation (hatching) of new eggs in the breeding habitat at time $s$ [1,75]. Since $t \geq s$, it follows that $Y(t, s)F(s)\phi(s)$ represents the distribution of new eggs at time $s$, and became adult at time $t$. Hence, the cumulative distribution of new eggs at time $t$, produced by all female adult mosquitoes ($\phi(s)$) introduced at a prior time $s = t$, is given by

$ $$\begin{align*} \Psi(t) = \int_{-\infty}^t Y(t, s)F(s)\phi(s)ds = \int_0^{\infty} Y(t, t-a)F(t-a)\phi(t-a)da. \end{align*}$$ $

Let $\mathbb{C}_{\omega}$ be the ordered Banach space of all $\omega$-periodic functions from $\mathbb{R}$ to $\mathbb{R}^9$, which is equipped with maximum norm and positive cone $\mathbb{C}_{\omega}^{+} \bigl\{\phi \in \mathbb{C}_\omega : \phi(t) \geq 0, \forall \, \, t \in \mathbb{R}\bigr\}$ [75]. Define a linear operator $L: \mathbb{C}_{\omega} \rightarrow \mathbb{C}_{\omega}$ [75]

$ $$\begin{align*} (L\phi)(t) = \int_0^{\infty} Y(t, t-a)F(t-a)\phi(t-a)\, da\, \, \forall \, \, t \in \mathbb{R}, \phi \in \mathbb{C}_{\omega}. \end{align*}$$ $

The vectorial reproduction ratio of the model (22) $(\mathcal{R}_{0TR})$ is then given by the spectral radius of $L$, (i.e., $\mathcal{R}_{0TR} = \rho(L)$). It can be verified that system (6) satisfy the assumptions $A1 -A7$ in [75]. Hence, the result below follows from Theorem 2.2 in [75].

Theorem 4.1. The trivial equilibrium $\left(\mathcal{T}^*_{0}\right)$, of the non-autonomous model (6), is LAS in $C([0], \mathbb{R}_{+}^{9})$ if $\mathcal{R}_{0TR} < 1$, and unstable if $\mathcal{R}_{0TR} > 1$.

The global asymptotic stability of the trivial equilibrium $\mathcal{T}^*_{0}$ of the model is now considered.

Theorem 4.2. The trivial equilibrium $\mathcal{T}^*_{0}$ of the non-autonomous model (1) is GAS in $C([0], \mathbb{R}_{+}^{9})$ whenever $\mathcal{R}_{0TR} < \, 1$.

The proof of Theorem 4.2, based on using comparison theorem [69], is given in Appendix B. The epidemiological implication of Theorem 4.2 is that the mosquito population (both immature and mature) can be effectively controlled (or eliminated) if the associated vectorial reproduction threshold, $\mathcal{R}_{0TR}$, can be brought to (and maintained at) a value less than or equal to unity.

---

4.2. Existence of non-trivial positive periodic solution

In this section, the possibility of the existence of a non-trivial positive periodic solution for the non-autonomous system (1) will be explored using uniform persistence theory [40,72,81,82] (see also [55]). Following and using notations in, Lou and Zhao [40], it is convenient to define the following sets ($X$, $X_0$ and $\partial X_0$)

$ $$\begin{align*} X =&\ \Omega, \\% C([0], \mathbb{R}_{+}^{9}), \\ X_0 =&\ \bigl\{\phi = (\phi_1, \phi_2, \phi_3, \phi_4, \phi_5, \phi_6, \phi_7, \phi_8, \phi_9) \in X : \phi_i(0) > 0\ {\rm for\ all \ } i \in [1,9]\bigr\}, \\ \partial X_0 =&\ X\setminus X_0 = \bigl\{\phi \in X : \phi_i(0) = 0 \, \, {\rm for \, some}\, \, i \in [1,9]\bigr\}. \end{align*}$$ $

Theorem 4.3. Consider the non-autonomous model (2.1). Let $\mathcal{R}_{0TR} > 1$. The model has at least has at least one positive periodic solution, and there exists a $\varphi > 0$ such that any solution $u(t, \phi)$ of the model with initial value $\phi \in X_0$ satisfies

$ $$\begin{align*} \liminf\limits_{t\rightarrow\infty}(E, L_1, L_2, L_3, L_4, P, V, W, U)(t) \geq (\varphi, \varphi, \varphi, \varphi, \varphi, \varphi, \varphi, \varphi, \varphi). \end{align*}$$ $

Proof. The proof is based on using the approach in [40,55]. Let $u(t, \phi)$ be the unique solution of model (1), with $u(0, \phi) = \phi$. Let $\Phi(t)\psi = u(t, \psi)$ and let $\mathcal{P} : X \rightarrow X$ be the Poincaré map associated with (1) i.e., $\mathcal{P}(\phi) = u(\omega, \phi)$ for all $\phi \in X$. Then, using similar approach as in Lemma (2.1), it is easy to see that $X_0$ is a positively invariant compact set. Hence, since solutions of model (1) are uniformly (ultimately) bounded, $\mathcal{P}$ is point dissipative [40]. It then follows from Theorem 1.1.2 in [82] that $\mathcal{P}$ admits a global attractor in X.

Thus, it suffices to show that model (1) is uniformly-persistent with respect to $(X_0, \partial X_0)$ [40]. It is convenient to define

$ $$\begin{split} K_\partial =&\ \bigl\{\phi \in \partial X_0 : \mathcal{P}^n(\phi) \in \partial X_0 {\rm\, \, for\, \, } n \geq 0\bigr\}, \\ D_1 =&\ \bigl\{\phi \in X : \phi_i(0) = 0\ {\rm for\ all \ } i \in [1,9]\bigr\}, \\ \partial X_0\setminus D_1 =&\ \bigl\{\phi \in X : \phi_i(0) \geq 0\ {\rm for\ some \ } i \in [1,9]\bigr\}. \end{split}$$ $

(23)

We claim that $K_\partial = D_1$ [40,55]. This claim can be proved by, first of all, seing that for any $\psi \in D_1$, $u_i(t, \psi) = 0$ for all $i \in [1,9]$, (hence, $D_1 \subset K_\partial$). Furthermore, for any $\psi \in \partial X_0 \setminus D_1$, we can choose $\psi_i(0) > 0$ for all $i = 1, \ldots, 9$, so that $u(t, \psi) \in X_0$. This implies that for any $\psi \in \partial X_0 \setminus D_1$, there exist some $n$ with $n\omega > t_0$, such that $\mathcal{P}^n(\phi) \notin \partial X_0$. Hence, $K_\partial \subset D_1$. This concludes the proof of the claim.

Thus, from (23), it can be verified that $\mathcal{P}^n(\phi)$, $n \geq 0$ contains exactly one trivial equilibrium:

$ $$\begin{align*} \mathcal{T}^*_0 = \bigl(0, 0, 0, 0, 0, 0, 0, 0, 0\bigr). \end{align*}$$ $

Hence, the set $\mathcal{A} := \{ \mathcal{T}^*_0\}$ is a compact and isolated invariant set for the Poincaré map $P$ in $K_\partial$ [40] and $\cup_{\phi \in K_\partial}\omega(\phi) = {\mathcal{A}}$ [40]. Furthermore, $\{\mathcal{A}\}$ does not form a cycle in $K_\partial$ (and hence in $\partial X_0$). In addition, since $M^*(t)$ is a positive periodic solution with respect to $X_0$, then it follows from the proof of Theorem 3.2 in [40] that there exists a $\epsilon > 0$ such that

$ $$\begin{align*} \limsup\limits_{t\rightarrow\infty}|\Phi(n\omega)\phi -\mathcal{T}_0| \geq \epsilon {\rm \ \, for \ \, all\ \, } \phi \in X_0. \end{align*}$$ $

Thus, it follows that $\mathcal{A}$ is an isolated invariant set for $\mathcal{P}$ in X and $W^s(\mathcal{A}) \cap X_0 = \emptyset$ where $W^s(\mathcal{A})$ is the stable manifold of $\mathcal{A}$ for $\mathcal{P}$ [40]. Hence, every trajectory in $K_\partial$ converges to $\mathcal{A}$, and $\{\mathcal{A}\}$ is acyclic in $K_\partial$ [82]. It then follows from Theorem 1.3.1 in [82] that $\mathcal{P}$ is uniformly persistent with respect to $X_0$. Thus, from Theorem 3.1.1 in [82], the periodic semiflow $\Phi(t): X \rightarrow X$ is also uniformly persistent to $X$ [40] where $\Phi(t)\psi = u_t(\psi)$. It then follows from Theorem 4.5 in [42] (see also Theorem 4.6 in [71] and Theorem 3.1 in [83]) that the model (1) admits a positive $\omega$-periodic solution $\mathcal{T}^*_1 = \Phi(t)\phi^*$ with $\phi^* \in X_0$.

It follows, from Theorem 4.5 in [42] (see also Theorem 2.1 in [84]), that $\mathcal{P}: X_0 \rightarrow X_0$ has a compact global attractor, denoted by $\mathcal{A}_0$. Hence, $\mathcal{A}_0$ is invariant for $\mathcal{P}$ (that is, $\mathcal{A}_0 = \mathcal{P}(\mathcal{A}_0) = \Phi(\omega)A_0$). Furthermore, using the notation in [55,84], let $\mathcal{A}^*_0 := \bigcup\limits_{t \in [0, \omega]}\Phi(t)\mathcal{A}_0$. Then, $\psi_i(0) > 0$ for all $\psi \in \mathcal{A}^*_0, \, i = 1, \ldots, 9$ [40]. Since $X_0$ is invariant, it follows that $\Phi(t)X_0 \subset X_0$. Thus, $\mathcal{A}^*_0 \subset X_0$ and $\limsup\limits_{t \rightarrow \infty} d(\Phi(t)\phi, \mathcal{A}^*_0) = 0$ for all $\phi \in X_0$ [40,84]. Also, it follows, by the continuity of $\Phi(t)\phi$ for $(t, \phi) \in [0, \infty) \times X_0$ and the compactness of $[0, \tau] \times \mathcal{A}_0$ [84], that $\mathcal{A}^*_0$ is compact in $X_0$ [40,84]. Thus, $\inf\limits_{\phi \in \mathcal{A}^*_0}d(\phi, \partial X_0) = \min\limits_{\phi \in \mathcal{A}^*_0} d(\phi, \partial X_0) > 0$ [40,84]. Consequently, there exists $\varphi > 0$ such that

$ $$\begin{align*} \liminf\limits_{t\rightarrow \infty}\, {\rm min } &\bigl(E(t, \phi), L_{1}(t, \phi), L_2(t, \phi), L_{3}(t, \phi), L_4(t, \phi), P(t, \phi), V(t, \phi), W(t, \phi), \\&U(t, \phi) \bigr) = \liminf\limits_{t\rightarrow\infty}d(\phi, \partial X_0) \geq \varphi, {\ \, \rm for \ all \ } \phi \in X_0. \end{align*}$$ $

In particular, $\liminf \limits_{t\rightarrow\infty}\min (\Phi(t)\phi^*) \geq \varphi$. Hence, $u_i(t, \phi) > 0$, $i = 1, \ldots, 9$ for all $t \geq 0$.

---

5. Numerical simulations

The non-autonomous model (6) is simulated to illustrate the effect of the two climate variables (temperature and rainfall) on the population dynamics of adult mosquitoes in a community. Suitable functional forms for the temperature-and rainfall-dependent functions, relevant to Anopheles mosquitoes (mostly given in [2,47,55,60]) as defined in Section 2.1, will be used. For these simulations, water temperature $(\hat{T})$ is taken to be $\hat{T} = T + 3^{0}$C. Furthermore, the simulations are carried out using the parameter values in Table 4 (with a fixed nutrient value of $N = 100000$).

The combined effect of mean monthly temperature and rainfall is assessed by simulating the non-autonomous model using various mean monthly temperature and rainfall values in the range $[16,40]^0$C and $[90-120]$ mm, respectively (the temperature ranges for most tropical and sub-tropical regions of the world lie within this temperature range [11]). The results obtained (as measured in terms of the total number of female adult mosquitoes), depicted in Figure 5, show that the total mosquito population (of a typical community with the aforementioned temperature and rainfall ranges) is maximized when the mean monthly temperature and rainfall values lie in the range $[20-25]^0$C and $[105-115]$ mm, respectively. Furthermore, simulations were carried out using weather (temperature and rainfall) data for three cities in Africa, namely, KwaZulu-Natal, South-Africa (Southern Africa); Lagos, Nigeria (Western Africa) and Nairobi, Kenya (Eastern Africa) (see profiles Tables 7, 8 and 9, respectively). While the peak mosquito abundance for KwaZulu-Natal (Figure 6a) and Lagos (Figure 6b) occur when temperature and rainfall values lie in the range $[22-25]^0$C, $[98-121]$ mm (occurring during the months of January, March, April, November and December) and $[24-27]^0$C $[113-255]$ mm (occurring during the months of May, July, August, September and October) respectively, the peak mosquito abundance for Nairobi (Figure 6c) occurs for temperature and rainfall ranges $[20.5-21.5]^0$C and $[50-120]$ mm (occurring during the months of January, February, March and April).

---

6. Conclusions

This study presents a new mathematical model for the population biology of the mosquito (the world's deadliest animal, which accounts for 80% of vector-borne diseases of humans). Some of the notable features of the new model are:

(ⅰ) incorporating four developmental stages of the mosquito larvae ($L_1, L_2, L_3, L_4$);

(ⅱ) including density dependence for the eggs oviposition process and larval mortality rates;

(ⅲ) including the dispersal states of female adult mosquitoes($U, V$ and $W$).

The model, which takes the form of a non-autonomous deterministic system of non-linear differential equations, is used to assess the impact of temperature and rainfall on the population dynamics of the mosquito. The main theoretical and epidemiological findings of this study are summarized below:

(ⅰ) The trivial equilibrium of the autonomous model (6) is globally-asymptotically stable whenever the associated vectorial reproduction number $(\mathcal{R}_0)$ is less than unity. For the case when the associated vectorial reproduction number $(\mathcal{R}_0)$ exceeds unity, the autonomous model (6) has at least one non-trivial equilibrium. Furthermore, it is shown that the autonomous model has a unique non-trivial equilibrium for the special case with no density-dependent larval mortality (i. e., $\delta_L = 0$). This unique non-trivial equilibrium is shown to be globally-asymptotically stable under a certain condition (i.e., when $1< \mathcal{R}_0<\mathcal{R}_0^C$) (this equilibrium bifurcates into a limit cycle, via a Hopf bifurcation at ${\mathcal R}_0={\mathcal R}_0^C$).

(ⅱ) Uncertainty and sensitivity analysis of the autonomous version of the model shows that the top five parameters that have the most influence on the dynamics of the model (with respect to various response functions) are the probability of female adult mosquito of type $W$ successfully taking a blood meal ($\alpha$), the natural mortality rate of female adult mosquitoes ($\mu_A$), the natural mortality rate of female larvae ($\mu_L$), the deposition rate of female eggs ($\psi_U$) and the maturation rate of female larvae from Stage 1 to Stage 2 ($\xi_1$). Hence, this study suggests that the population of adult mosquito in a community can be effectively-controlled using mosquito-reduction strategies, as well as personal protection against mosquito bites.

(ⅲ) The trivial periodic solution of the non-autonomous model (1) is shown to be globally-asymptotically stable, whenever the spectral radius of a certain linear operator (denoted by $\mathcal{R}_{0TR}$) is less than unity. Furthermore, it is shown, using uniform persistence theory, that the non-autonomous model (1) has at least one positive periodic solution whenever $\mathcal{R}_{0TR} > 1$.

Numerical simulations of the non-autonomous model, using relevant functional forms (given in Section 2.1) and parameter values associated with the Anopheles species (which causes malaria in humans), show the following:

(ⅰ) For mean monthly temperature and rainfall values in the range $[10,40]^0$C and $[90-120]$ mm, respectively, the peak mosquito abundance lie in the range $[20-25]^0$C and $[105-115]$ mm, respectively.

(ⅱ) For mean monthly temperature and rainfall data for three cities in Africa, namely, KwaZulu-Natal, South-Africa; Lagos, Nigeria and Nairobi, Kenya (Tables 7, 8 and 9). The peak mosquito abundance for KwaZulu-Natal (Figure 6a) and Lagos (Figure 6b) occur when temperature and rainfall values lie in the range $[22-25]^0$C, $[98-121]$ mm (occurring during the months of January, March, April, November and December) and $[24-27]^0$C, $[113-255]$ mm (occurring during the months of May, July, August, September and October) respectively. The peak mosquito abundance for Nairobi (Figure 6c) occurs for temperature and rainfall ranges $[20.5-21.5]^0$C and $[50-120]$ mm (occurring during the months of January, February, March and April).

---

Acknowledgments

One of the authors (ABG) is grateful to National Institute for Mathematical and Biological Synthesis (NIMBioS) for funding the Working Group on Climate Change and Vector-borne Diseases (VBDs). NIMBioS is an Institute sponsored by the National Science Foundation, the U.S. Department of Homeland Security, and the U.S. Department of Agriculture through NSF Award #EF-0832858, with additional support from The University of Tennessee, Knoxville.

---

7. Appendices

---

7.1. Appendix A: Coefficients of Equations (14) and (16)

---

7.1.1. Coefficients of Equation (14)

$ $$\begin{align*} b_0\ =&\ \biggl(\dfrac{1}{\mathcal{R}_0} -1\biggr)Q_1X_6, \\ \end{align*}$$ $

$ $$\begin{align*} b_1 \ =&\ \biggl[\biggl(\dfrac{1}{\mathcal{R}_0}-1\biggr)Q_1X_3 + Q_2 X_5 X_6 \biggr]\delta_L + 1, \\ b_2 \ =&\ \biggl[\biggl(\dfrac{1}{\mathcal{R}_0}-1\biggr)Q_1X_2 + Q_2\bigl(X_5 X_3 + X_4 X_6\bigr)\biggr] \left(\delta_L\right)^2, \\ b_3 \ =&\ \biggl[\biggl(\dfrac{1}{\mathcal{R}_0}-1\biggr)Q_1 + Q_2\bigl( X_1 X_6 + X_2 X_4 + X_3 X_4 \bigr)\biggr]\left(\delta_L\right)^3, \\ b_4 \ =&\ Q_2\bigl(X_1 X_3 + X_2X_4 + X_5 + X_6 \bigr)\left(\delta_L\right)^4, \\ b_5 \ =&\ Q_2 \bigl(X_1 X_2 + X_3 + X_4\bigr)\left(\delta_L\right)^5, \\ b_6 \ =&\ Q_2 \left( X_1 + X_2\right) \left(\delta_L\right)^6, \\ b_7 \ =&\ Q_2\left(\delta_L\right)^7, \end{align*}$$ $

where,

$ $$\begin{align*} Q_1 =& \ \dfrac{\sigma_EC_PDK_U}{\alpha\tau^*_W\eta^*_VB}, \ \ Q_2 = \dfrac {K_{{U}}(C_PD)^{2}\sigma_EC_E}{(\alpha\tau^*_W\eta^*_VB)^{2}\psi_U}, \ \ X_1 = C_1+C_2+C_3+C_4, \\ X_2 =& \ C_2+C_3+C_4+\xi_1, \ \ X_3 = C_2C_3 + C_2C_4 + C_3C_4 + C_3\xi_1+C_4\xi_1 + \xi_1\xi_2, \\ X_4 =& \ C_1C_2+C_1C_3+C_1C_4+C_2C_3+C_2C_4+C_3C_4, \ \ X_5 =\ C_1C_2C_3 + C_1C_2C_4 \\& + C_1C_3C_4 + C_2C_3C_4, \ \ X_6 = C_2C_3C_4+C_3C_4\xi_1+C_4\xi_1\xi_2+\xi_1\xi_2\xi_3, \end{align*}$$ $

with $\tau^*_W, \, \eta^*_V, \, B, \, C, \, C_i \ (i=E, P, 1, \ldots, 7)$ and $D$ as defined in (7).

---

7.1.2. Coefficients of Equation (16)

$ $$\begin{align*} A_1 =& D\bigl( C_1 C_2 C_3 C_4 C_E + C_1 C_2 C_3 C_4 C_P + C_1 C_2 C_3 C_E C_P + C_1 C_2 C_4 C_E C_P \\&+ C_1 C_3 C_4 C_E C_P + C_2 C_3 C_4 C_E C_P\bigr) + C\bigl(C_5C_6 + C_5C_7 + C_6C_7\bigr) >0, \end{align*}$$ $

$ $$\begin{split} A_2 =& \sum\limits_{i=1}^{2}C_i\sum\limits_{j=i+1}^{3}C_j\sum\limits_{k=j+1}^{4}C_j\sum\limits_{l=k+1}^{5}C_l\sum\limits_{m=l+1}^{6}C_m\sum\limits_{n=m+1}^{7}C_n\left(C_E + C_P + \sum\limits_{q=n+1}^{7}C_q\right) \\&+ C_EC_P\biggl(\sum\limits_{i=1}^{3}C_i\sum\limits_{j=i+1}^{4}C_j\sum\limits_{k=j+1}^{5}C_k\sum\limits_{l=k+1}^{6}C_l\sum\limits_{m=l+1}^{7}C_m -\alpha\tau^*_W\eta^*_V\gamma_U\sum\limits_{i}^{3}C_i \\&\sum\limits_{j=i+1}^{4}C_j\biggr) -\alpha\tau^*_W\eta^*_V\gamma_U\sum\limits_{i=1}^{2}C_i\sum\limits_{j=i+1}^{3}C_j\sum\limits_{k=j+1}^{4}C_k\left( \sum\limits_{l=k+1}^{4}C_l + C_E + C_P\right) > 0, \\ A_3 =& \sum\limits_{i=1}^{3}C_i\sum\limits_{j=i+1}^{4}C_j\sum\limits_{k=j+1}^{5}C_j\sum\limits_{l=k+1}^{6}C_l\sum\limits_{m=l+1}^{7}C_m\left(C_E + C_P + \sum\limits_{n=m+1}^{7}C_n\right) \\&+ C_EC_P\left(\sum\limits_{i=1}^{4}C_i\sum\limits_{j=i+1}^{5}C_j\sum\limits_{k=j+1}^{6}C_k\sum\limits_{l=k+1}^{7}C_l -\alpha\tau^*_W\eta^*_V\gamma_U\sum\limits_{i}^{4}C_i\right) \\& -\alpha\tau^*_W\eta^*_V\gamma_U\sum\limits_{i=1}^{3}C_i\sum\limits_{j=i+1}^{4}C_j\left( \sum\limits_{k=j+1}^{4}C_k + C_E + C_P\right) > 0, \\ A_4 =& \sum\limits_{i=1}^{4}C_i\sum\limits_{j=i+1}^{5}C_j\sum\limits_{k=j+1}^{6}C_j\sum\limits_{l=k+1}^{7}C_l\left(C_E + C_P + \sum\limits_{m=l+1}^{7}C_m\right) \\&+ C_EC_P\left(\sum\limits_{i=1}^{5}C_i\sum\limits_{j=i+1}^{6}C_j\sum\limits_{k=j+1}^{7}C_k -\alpha\tau^*_W\eta^*_V\gamma_U\right) \\&-\alpha\tau^*_W\eta^*_V\gamma_U\sum\limits_{i=1}^{3}C_i\left( \sum\limits_{j=i+1}^{4}C_j + C_E + C_P\right) > 0, \\ A_5 =& \sum\limits_{i=1}^{5}C_i\sum\limits_{j=i+1}^{6}C_j\sum\limits_{k=j+1}^{7}C_j\left(C_E + C_P + \sum\limits_{l=k+1}^{7}C_l\right) + C_EC_P\sum\limits_{i=1}^{6}C_i\sum\limits_{j=i+1}^{7}C_j \\&-\alpha\tau^*_W\eta^*_V\gamma_UC > 0, \\ A_6 =& \sum\limits_{i=1}^{6}C_i\sum\limits_{j=i+1}^{7}C_j\left(C_E + C_P + \sum\limits_{k=j+1}^{7}C_k\right) + C_EC_P\sum\limits_{i=1}^{7}C_i -\alpha\tau^*_W\eta^*_V\gamma_U > 0, \\ A_7 =& \sum\limits_{i=1}^{7}C_i\left( C_E + C_P + \sum\limits_{j=i+1}^{7}C_j \right) > 0, \\ A_8 =& C_E + C_P + \sum\limits_{i=1}^{7}C_i > 0, \end{split}$$ $

where $\tau^*_W, \, \eta^*_V, \, B, \, C, \, C_i \ (i=E, P, 1, \ldots, 7)$ and $D$ are as defined in (7).

---

7.2. Appendix B: Proof of Theorem 3.11

Proof. Consider the model (6) with $\delta_L = 0$ and $sgn(U^{**} -U(t)) = sgn(E^{**} -E(t))$ for all $t \geq 0$. Furthermore, let $1 < \mathcal{R}_0 < \mathcal{R}^C_0$, so that the unique non-trivial equilibrium $\mathcal{T}_1$ exists and is LAS (in line with Theorem 3.10). Consider, further, the following non-linear Lyapunov function of Goh-Volterra type:

$ $$\begin{align*} \mathcal{K}_2 =& \ E -E^{**}\ln E + d_1(L_1 -L_1^{**}\ln L_1) + d_2(L_2 -L_2^{**}\ln L_2) \\& + d_3(L_3 -L_3^{**}\ln L_3) + d_4(L_4 -L_4^{**}\ln L_4) + d_5(P -P^{**}\ln P) \\& + d_6(V -V^{**}\ln V) +\ d_7(W -W^{**}\ln W) + d_8(U -U^{**}\ln U), \end{align*}$$ $

where,

$ $$\begin{split} & d_1 = \dfrac{C_E}{\sigma_E}, \ d_2 = \dfrac{C_1C_E}{\xi_1\sigma_E}, \ d_3 = \dfrac{C_2C_1C_E}{\xi_2\xi_1\sigma_E}, \ d_4 = \dfrac{C_3C_2C_1C_E}{\xi_3\xi_2\xi_1\sigma_E}, \\&d_5 = \dfrac{C_4C_3C_2C_1C_E}{\xi_4\xi_3\xi_2\xi_1\sigma_E}, \ d_6 = \dfrac{C}{B}, \ d_7 = \dfrac{C_5C}{\eta^*_VB}, \ d_8 = \dfrac{C_6C_5C}{\alpha\tau^*_W\eta^*_VB}, \end{split}$$ $

(24)

with $\tau^*_W, \, \eta^*_V, \, B, \, C, \, C_i \ (i=E, P, 1, \ldots, 7)$ and $D$ as defined in (7). The following steady state relations will be used to simplify the derivative of the Lyapunov function $\mathcal{K}_2$.

$ $$\begin{split} & C_EE^{**} = \psi_U \left(1-\dfrac{U^{**}}{K_U}\right)U^{**}, C_1 L^{**}_1 = \sigma_E E^{**}, \ C_i L^{**}_i = \xi_{(i-1)}L^{**}_{(i-1)} ;\, \ i = 2, 3, 4, \\ & C_P P^{**} = \xi_{4} L_4^{**}, C_5V^{**} = \sigma_P P^{**} + \gamma_U U^{**}, \ C_6 W^{**} = \eta^*_VV^{**}, C_7 U^{**} = \alpha\tau^*_W W^{**}. \end{split}$$ $

(25)

The Lyapunov derivative of $\mathcal{K}_2$ is

$ $$\begin{split} &\dot{\mathcal{K}}_2 = \ \biggl(1 -\dfrac{E^{**}}{E} \biggr)\left[\psi_U \left(1-\dfrac{U}{K_U}\right)U-C_EE\right] + d_1\biggl(1 -\dfrac{L_1^{**}}{L_1} \biggr)\bigl[\sigma_EE-C_1L_1 \bigr] \\& + d_2 \biggl(1 -\dfrac{L_2^{**}}{L_2} \biggr) \bigl[\xi_1L_1-C_2L_2 \bigr] + d_3 \biggl(1 -\dfrac{L_3^{**}}{L_3} \biggr) \bigl[\xi_2L_2-C_3L_3 \bigr] + d_4 \biggl(1 -\dfrac{L_4^{**}}{L_4} \biggr)\\& \bigl[\xi_3L_3-C_4L_4 \bigr] + d_5 \biggl(1 -\dfrac{P^{**}}{P} \biggr) \bigl[\xi_4L_4-C_PP \bigr] + d_6 \biggl(1 -\dfrac{V^{**}}{V} \biggr)\bigl[\sigma_PP-C_5V \\& + \gamma_UU \bigr] + d_7 \biggl(1 -\dfrac{W^{**}}{W} \biggr) \bigl[\eta^*_VV-C_6W \bigr] + d_8 \biggl(1 -\dfrac{U^{**}}{U} \biggr) \bigl[\alpha\tau^*_WW-C_7U \bigr], \end{split}$$ $

(26)

Substituting (24) and (25) into (26), and simplifying, gives

$ $$\begin{split} \dot{\mathcal{K}}_2 =& \ -\dfrac{\psi_UU}{EK_U}(E^{**} -E)(U^{**} -U) + \gamma_Ud_6U^{**}\biggl(3 -\dfrac{UV^{**}}{U^{**}V} -\dfrac{VW^{**}}{V^{**}W} -\dfrac{U^{**}W}{UW^{**}} \biggr) \\& + C_EE^{**} \biggl( 9 -\dfrac{L_1^{**}E}{L_1E^{**}} -\dfrac{L_2^{**}L_1}{L_2L_1^{**}} -\dfrac{L_3^{**}L_2}{L_3L_2^{**}} -\dfrac{L_4^{**}L_3}{L_4L_3^{**}} -\dfrac{P^{**}L_4}{PL_4^{**}} -\dfrac{V^{**}P}{VP^{**}} \\& -\dfrac{W^{**}V}{WV^{**}} -\dfrac{U^{**}W}{UW^{**}} -\dfrac{E^{**}U}{EU^{**}} \biggr) . \end{split}$$ $

(27)

The first term of (27) is automatically negative in $\Omega\setminus\{\mathcal{T}_0\}$, since $sgn(U^{**} -U(t)) = sgn(E^{**} -E(t))$ for all $t \geq 0$. Furthermore, since the arithmetic mean exceeds the geometric mean, it follows that the second and third term of (27) are also negative. Hence, $\dot{\mathcal{K}}_2 \leq 0$. The proof is concluded as in the proof of Theorem 3.2.

The assumption $sgn(U^{**} -U(t)) = sgn(E^{**} -E(t))$ for all $t \geq 0$ can be justified by noting the fact that, in order to maintain a non-trivial equilibrium for the adult mosquito population, it is necessary that each reproducing female mosquito must produce at least one egg during its entire reproductive life period. Thus, if the population of the reproduction female adult mosquitoes $(U(t))$ increases, so does the population of eggs laid $(E(t))$ for all $t \geq 0$.

---

7.3. Appendix C: Proof of Theorem 4.2

Proof. Consider the non-autonomous model (1) with $\mathcal{R}_{0TR}<1$. Using the fact that (and noting that $ L(t) = \sum\limits_{i = 1}^{4} L_i(t)$),

$ $$\begin{align*}\psi_U\biggl( 1 -\dfrac{U(t)}{K_U} \biggr) U(t) \leq \, \psi_U U(t) \ ( {\rm\, since \ } K_U > U(t) \ {\rm for \ all \ } t \geq 0), \end{align*}$$ $

and,

$ $$\begin{align*} C_i(t) + \delta_LL(t) \, \geq\, C_i(t), {\rm \ for \ all\ } t \geq 0, \end{align*}$$ $

it follows that the non-autonomous model (1), subject to the aforementioned assumptions, can be re-written as

$ $$\begin{split} \frac{d E}{d t} \ \leq&\ \psi_U U -C_E(t) E, \\ \frac{d L_1}{d t} \ \leq&\ \sigma_E(t) E -C_1(t) L_1, \\ \frac{d L_i}{d t} \ \leq&\ \xi_{(i-1)}(t)L_{(i-1)} -C_i(t) L_i \ ;\ i = 2, 3, 4, \\ \frac{d P}{d t} \ =& \ \xi_{4}(t) L_4 -C_P(t) P, \\ \frac{d V}{d t} \ =& \ \sigma_P(t) P + \gamma_U U(t) -C_5(t) V, \\ \frac{d W}{d t} \ =& \ \eta^*_V V -C_6(t) W, \\ \frac{d U}{d t} \ =& \ \alpha\tau^*_W W -C_7(t) U. \end{split}$$ $

(28)

Following [75], the equations in (28), with equalities used in place of the inequalities, can be re-written in terms of the next generation matrices $F(t)$ and $V(t)$, as follows

$ \frac{d X(t)}{d t} = [F(t)-V(t)]X(t). $

(29)

It follows, from Lemma 2.1 in [80], that there exists a positive and bounded $\omega$-periodic function, $x(t) = \bigl(\bar{E}(t), \bar{L}_1(t), \bar{L}_2(t), \bar{L}_3(t), \bar{L}_4(t), \bar{P}(t), \bar{V}(t), \bar{W}(t), \bar{U}(t) \bigr)$, such that

$ $$\begin{align*} X(t)=e^{\theta t}x(t), \ {\rm with}\ \, \theta = \dfrac{1}{\omega}\ln\rho\bigl[\phi_{F-V}(\omega)\bigr], \end{align*}$$ $

is a solution of the linearized system (28). Furthermore, it follows from Theorem 2.2 in [75] that $\mathcal{R}_{0TR}<1$ if and only if $\rho\bigl[\phi_{F-V}(\tau)\bigr]<1$. Hence, $\theta$ is a negative constant. Thus, $X(t) \rightarrow 0$ as $t \rightarrow \infty$. Thus, the unique trivial solution of the linear system (28), given by $X(t) = 0$, is GAS [40,55,66]. For any non-negative initial solution $(E, L_1, L_2, L_3, L_4, P, V, W, U)(0))^T$ of the system (29), there exists a sufficiently large $Q^* > 0$ such that [55,66],

$ $$\begin{align*} \left((E, L_1, L_2, L_3, L_4, P, V, W, U)(0)\right)^T \leq Q^*\left((\bar{E}, \bar{L}_1, \bar{L}_2, \bar{L}_3, \bar{L}_4, \bar{P}, \bar{V}, \bar{W}, \bar{U})(0)\right)^T. \end{align*}$$ $

Thus, it follows, by comparison theorem [37,69], that

$ $$\begin{align*} \bigl(E(t), L_1(t), L_2(t), L_3(t), L_4(t), P(t), V(t), W(t), U(t)\bigr) \leq Q^*X(t) {\rm \, \, for\ all\ \, } t > 0, \end{align*}$$ $

where $Q^*X(t)$ is also a solution of (29).

Hence,

$ $$\begin{align*} \bigl(E(t), L_1(t), L_2(t), L_3(t), L_4(t), P(t), V(t), W(t), U(t)\bigr) \rightarrow \bigl(0, 0, 0, 0, 0, 0, 0, 0, 0\bigr), \end{align*}$$ $

as $t \rightarrow \infty$.

---