Mathematical analysis of a weather-driven model for the population ecology of mosquitoes.

A new deterministic model for the population biology of immature and mature mosquitoes is designed and used to assess the impact of temperature and rainfall on the abundance of mosquitoes in a community. The trivial equilibrium of the model is globally-asymptotically stable when the associated vectorial reproduction number (R0) is less than unity. In the absence of density-dependence mortality in the larval stage, the autonomous version of the model has a unique and globally-asymptotically stable non-trivial equilibrium whenever 1 andlt;R0 andlt;RC0 (this equilibrium bifurcates into a limit cycle, via a Hopf bifurcation at R0=RC0). Numerical simulations of the weather-driven model, using temperature and rainfall data from three cities in Sub-Saharan Africa (Kwazulu Natal, South Africa; Lagos, Nigeria; and Nairobi, Kenya), show peak mosquito abundance occurring in the cities when the mean monthly temperature and rainfall values lie in the ranges [22-25]0C, [98-121] mm; [24-27]0C, [113-255] mm and [20.5-21.5]0C, [70-120] mm, respectively (thus, mosquito control efforts should be intensified in these cities during the periods when the respective suitable weather ranges are recorded).


Abstract.
A new deterministic model for the population biology of immature and mature mosquitoes is designed and used to assess the impact of temperature and rainfall on the abundance of mosquitoes in a community. The trivial equilibrium of the model is globally-asymptotically stable when the associated vectorial reproduction number (R 0 ) is less than unity. In the absence of density-dependence mortality in the larval stage, the autonomous version of the model has a unique and globally-asymptotically stable non-trivial equilibrium whenever 1 < R 0 < R C 0 (this equilibrium bifurcates into a limit cycle, via a Hopf bifurcation at R 0 = R C 0 ). Numerical simulations of the weather-driven model, using temperature and rainfall data from three cities in Sub-Saharan Africa (Kwazulu Natal, South Africa; Lagos, Nigeria; and Nairobi, Kenya), show peak mosquito abundance occurring in the cities when the mean monthly temperature and rainfall values lie in the ranges [22 − 25] [20.5 − 21.5] 0 C, [70 − 120] mm, respectively (thus, mosquito control efforts should be intensified in these cities during the periods when the respective suitable weather ranges are recorded).
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]).

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) 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.
In (1), R = R(t), T = T (t), andT =T (t) denote rainfall, air temperature and 60 KAMALDEEN OKUNEYE, AHMED ABDELRAZEC AND ABBA B. GUMEL water temperature at time t, respectively. Typically, a sinusoidal function, such as (where T 0 is the mean annual temperature, T 1 captures variation about the mean, and ω and θ 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),T and R(t) are assumed to be continuous, bounded, positive and ω-periodic. Furthermore, the parameters and µ A (T ) are non-negative, ω-periodic, continuous and bounded functions defined represents the density-dependent eggs oviposition rate (where K U > U (t) for all t ≥ 0 is the environmental carrying capacity of female adult mosquitoes and ψ U (T ) is the temperature-dependent egg deposition rate). Eggs hatch into the first larval stage at a rainfall-and temperaturedependent rate σ E (R, T ). Larvae in Stage i mature into Stage i + 1 at a rate ξ 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 ξ 4 (N, R, T ). It should be emphasized that the maturation rates for the larval stages (ξ 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 σ P (R,T ). These female adult mosquitoes quest for blood meal at the human habitat at a rate η V (and become female adult mosquitoes of type W ) [51]. The term 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 τ W [51,52]. This interaction can be successful with probability α ∈ [0, 1], so that questing mosquitoes successfully obtain blood meals and become vectors of type U (at the rate ατ W ) which, in turn, return to become female adult mosquitoes of type V at a rate γ 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 µ E (T ), µ L (T ), µ P (T ), µ A (T ) represent, respectively, the temperaturedependent natural death rate for female eggs, female larvae, female pupae and female adult mosquitoes, and δ L L 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,47,59,60], characterizes the female Anopheles mosquitoes (which transmits malaria in humans). The per-capita rate of deposition of female eggs (ψ U (T )) defined using the quadratic function used in [47], is given by Similarly, following [47], the per-capita death rate of the female adult mosquitoes (µ A (T )) is defined as µ A (T ) = − ln −0.000828T 2 + 0.0367T + 0.522 . Furthermore, following [60], the per-capita death rate of female eggs (µ E (T )), female larvae (µ L (T )) and female pupae (µ P (T )) are defined, respectively, as Similarly, following [60], the per-capita maturation rate of eggs (into larvae) and pupae (into female adult mosquitoes) are given by is the average duration spent in Stage i, given by and p i = p i (R,T ), (where 0 ≤ p i (R,T ) < 1) is the daily survival probability of immature mosquitoes in Stage i (assumed to be determined by the mean daily water temperature (T ( • C)) and cumulative daily rainfall (R (mm))), so that with p i (T ) = exp (−µ i (T )) (it should be stated that, in line with [60], the definition of p i (R,T ) = p i (R)p i (T ) emphasizes the assumed independence of temperature  Table 2. Range of values of temperature-dependent parameters in the temperature range [16,40] 0 C. and rainfall with each other). Following [59] (Supplemental Material), the rainfalldependent 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 where p M i is the peak daily survival probability of immature mosquito stage i (i = E = eggs; i = L j = larvae in Stage j; j = {1, 2, 3, 4}; i = P = 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 (ξ j ; j = 1, 2, 3, 4), assumed to depend on amount of available nutrients (N ), rainfall (R) and water temperature (T ), is given by where ξ j (N ) = e j N , with e j representing the rate of nutrients intake for female larvae in stage j. Furthermore, ξ j (R,T ) is given by with p j having similar definition as above, d j = d L (T ) = 1/µ L (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 byT = T + θ T , where θ 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 ≤ T (t),T (t) ≤ 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 (ψ U (T )), per capita death rate of female eggs (µ E (T )), per capita death rate of female larvae (µ L (T )), per capita death rate of female pupae (µ P (T )) and per capita death rate of female adult mosquitoes (µ 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: (i) aquatic stages of the mosquito (i.e., adding the E, L 1 , L 2 , L 3 , L 4 and P classes); (ii) 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: (i) incorporating the effect of climate variable (temperature and rainfall); (ii) using logistic growth rate for egg oviposition rate (a constant rate was used in [41]); (iii) 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 (i) including the dynamics of adult female mosquitoes (i.e., including compartments U , V and W ); (ii) including four larval classes (a single larval class is considered in [1]); (iii) 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 T =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, µ , 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) 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 ω-periodic function, ψ U (t) ∈ C 1 (0, R + ), there exists a positive number, h 0 , such that Lemma 2.1. For any φ ∈ Ω = C([0, R 9 + ]), the model (1) has a unique non-negative solution through φ, and all solutions are uniformly-bounded. (1) can then be re-written as: where, is continuous and Lipschitzian (with respect to φ in each compact set in R × Ω) [40]. Hence, there is a unique solution of system (1) through (0, φ). It should be noted that f i (t, π) ≥ 0 whenever π ≥ 0 and π i (0) = 0 [40]. Hence, it follows (from Remark 5.2.1 in [69]) that Ω 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 ≥ 0. Moreover, where M * (t) is the unique periodic solution oḟ given by, Thus, all solutions of the model (1) are ultimately-bounded [40]. Moreover, it follows from (5) thatṀ * < 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., σ where, now, 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 α < 1): The autonomous model (6) has a trivial equilibrium solution, denoted by T 0 , given by The linear stability of T 0 (in Ω) is obtained by using the next generation matrix [23,73] for the system (1). Using the notation in [73], the non-negative matrix F and the non-singular matrix V, for the new egg deposition terms and the remaining transfer terms, are, respectively, given (at the trivial equilibrium, T 0 ) by where 0 denotes a zero matrix of order 3, and It follows from [73] that the associated vectorial reproduction number of the autonomous model (6) [63], denoted by R 0 = ρ(FV −1 ), is given by (where ρ is the spectral radius of the next generation matrix FV −1 ) where τ * W , η * V , B, C, C i (i = E, P, 1, . . . , 7) and D are as defined in (7). The threshold quantity, 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].
Proof. Consider the Lyapunov function It is convenient to define Thus, the Lyapunov derivative is given bẏ where τ * W , η * V , B, C, C i (i = E, P, 1, . . . , 7) and D are as defined in (7). Thus, it follows that, for R 0 ≤ 1 in Ω, the Lyapunov derivativeK 1 < 0. Furthermore, it follows from the LaSalle's Invariance Principle (Theorem 6.4 of [39]) that the maximal invariant set contained in Theorem 3.2 shows that the mosquito population (both immature and mature) will be effectively controlled (or eliminated) if the associated vectorial reproduction threshold, R 0 , can be brought to (and maintained at) a value less than or equal to unity.
Theorem 3.4. The model (6) with δ L = 0 has a unique non-trivial equilibrium whenever 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 densitydependent mortality rate for larvae (i.e., δ L = 0) so that the autonomous model (6) has a unique non-trivial equilibrium (T 1 ) when R 0 > 1. Linearizing the autonomous model (6), with δ L = 0, at T 1 gives: The eigenvalues of the matrix J (T 1 ) satisfy the following polynomial: where τ * W , η * V and C i (i = E, P, 1, . . . , 7) are as defined in (7) and A i (i = 1, . . . , 8) are positive constants given in Appendix A2. It is convenient to re-write the polynomial (16) as where, and, so that, The asymptotic stability of 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 × n complex matrix, and let E k be the sum of all principal minors of A of order k, k ∈ < n >.
Let Ω(A) be the n × n Hurwitz matrix of A and assume that Ω(A) is real. Then A is stable if and only if all leading principal minors of Ω(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 ≥ m. The Bézoutiant defined by a(x) and b(x) is the bilinear form The symmetric matrix (b ik ) n−1 0 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 Theorem 3.7. (Liénard-Chipart) [31] Let f (x) = x n − a n x n−1 − . . . − a 1 be a polynomial with real coefficients, and let a n−1 = −1. Define the polynomials 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 ∈ < n >. We claim the following result.
Using the Routh-Hurwitz Criterion (Theorem 3.5), the principal minors, ∆ k (k = 1, 2, 3), of the associated Hurwitz matrix for G(λ) are Thus, all the roots of G(λ) have negative real part. Hence, all nine roots of F (λ)G(λ) have negative real part. Remark 1. It follows from Lemma 3.9 and Theorem 3.7 that the corresponding Bézout matrix of F (λ)G(λ) is positive-definite [31].
Remark 2. Consider P 9 (λ) = F (λ)G(λ) + CD(R 0 − 2). Then, P 9 (λ) ≤ F (λ)G(λ) whenever 1 < R 0 ≤ 2. Thus, it follows from Lemma 3.9 that all nine roots of P 9 (λ) have negative real part whenever 1 < R 0 ≤ 2 (hence, T 1 is LAS whenever 1 < R 0 ≤ 2). Furthermore, consider the characteristic polynomial P 9 (λ) given in (16). Let A 0 = CD(R 0 − 1), with C and D as defined in (7), so that To apply Theorem 3.7, let and, Thus, it follows from Definition 3.6 that the corresponding Bézout matrix of P 9 (λ), denoted by B h,g (P 9 ), is given by 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 × k sub-matrix of B h,g (P 9 ) is positive). It is convenient to introduce the following notations: } are the entries of the corresponding Bézout matrix of the polynomial F (λ)G(λ) (denoted by B h,g (F G)) and P 9 (λ) (clearly, 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 (F G). That is, in terms of b where C and D are as defined in (7) and K = (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 < R 0 ≤ 2. Furthermore, the Bézout matrix, B h,g (P 9 ), can be re-written as (after row-column operations) where 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 ∆ (P9) k 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 is positive whenever the following inequality holds: and Z 1 is the positive constant such that ∆ such that the k th principal minor ∆ (P9) k is positive whenever 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 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).
and unstable whenever R 0 > R C 0 . The results above (Theorem 3.4 and Theorem 3.10) show that the condition R 0 > 1 defines the existence of a unique non-trivial equilibrium T 1 of the model (6) with δ L = 0 (which is LAS if 1 < R 0 < 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 ≥ 0 and δ L = 0. Then, the non-trivial equilibrium (T 1 ) of the model (6) with δ L = 0 is GAS in Ω \ {T 0 } whenever 1 < R 0 < 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 (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 δ 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 δ L = 0 is ψ U = ψ * U = 107.889493160695073 (so that, ∆ 4 = 0). This is equivalent to R 0 = R C 0 = 4.5573. Therefore, for this particular set of parameter values, the non-trivial equilibrium (T 1 ) is LAS for 1 < R 0 < 4.5573, and unstable whenever R 0 > 4.5573. 3.3. Hopf Bifurcation analysis: Special case. Consider the model (6) with δ L = 0 and R 0 > 1 (so that 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 δ L = 0, evaluated at T 1 , has a pair of purely imaginary eigenvalues. That is, when P 9 (λ) 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].
Proof. To prove Theorem 3.12, it is sufficient to establish the transversality condition [20]. Let ψ U = ψ * U be a bifurcation parameter (and all other parameters of the model (6) are fixed). Then R 0 = R C 0 = 2 + Z 4 . Since Z 4 < Z i ( i = 1, 2, 3), it follows from Theorem 3.10 that ∆ can be re-written as Hence, ∆ (P9) 4 (ψ U ) = 0 if and only if ψ U = ψ * U . Furthermore, it can be verified that d∆ where 'Tr' and 'Adj' denote, respectively, the trace and adjoint of a matrix. Similarly, let µ A be a bifurcation parameter (and all other parameters of the model (6) are fixed). Thus, Theorem 3.12 shows that sustained oscillations are possible, with respect to the autonomous model (6) with δ L = 0, whenever R 0 = 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 (ψ U and µ 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 δ L = 0, which summarizes the main results obtain in Section 3, will be generated in the µ A − ψ U plane as follows: (i) Solving for ψ U from R 0 = 1 gives the following equation for ψ l U (depicted in Figure 4): (ii) Solving for ψ U from ∆ (P9) 4 = 0 (and fixing all parameters of the models (using their values as in Figure 4), except the parameters, µ A and ψ U ) give the following curve ∆ (P9) 4 = 0: where B, C and D are as defined in (7) and Z. The curves l and H (depicted in Figure 4) divide the µ A − ψ U plane into three distinct regions, namely D 1 , D 2 and D 3 , given by: The regions can be described as follows (see also Table 3): (i) Region D 1 : In this region, R 0 ≤ 1. Hence, in this region (note that δ L = 0), the trivial equilibrium (T 0 ) is globally-asymptotically stable (in line with Theorem 3.2). (ii) Region D 2 : Here, 1 < R 0 < R C 0 . Thus, the model has two equilibria, namely the unstable trivial equilibrium (T 0 ) and the locally-asymptotically stable nontrivial equilibrium (T 1 ). The model undergoes a Hopf bifurcation whenever R 0 = R C 0 . (iii) Region D 3 : In this region, R 0 > R C 0 . Thus, the model has the unstable trivial equilibrium (T 1 ), unstable non-trivial equilibrium and a stable limit cycle. Table 3. Stability properties of the solutions of the autonomous model (6). 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.

Unstable Unstable Yes
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 (α), the natural mortality rate of female adult mosquitoes (µ A ) and the natural mortality rate of female larvae (µ 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 (µ L ), the deposition rate of female eggs (ψ U ) and the maturation rate of female larvae from Stage 1 to Stage 2 (ξ 1 ). Furthermore, using the population of  Table 4. Values and ranges of the parameters of the autonomous model (6).
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 (R 0 ) as the response function, the top three PRCC-ranked parameters are the natural mortality rate of female larvae (µ L ), the deposition rate of female eggs (ψ U ) and the natural mortality rate of female adult mosquitoes (µ 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 (α), the natural mortality rate of female adult mosquitoes (µ A ), the natural mortality rate of female larvae (µ L ), the deposition rate of female eggs (ψ U ) and the maturation rate of female larvae (ξ i ). The effect of the aforementioned most dominant parameters (α, ψ U , ξ i , µ L and µ A ) on the population dynamics of the mosquito and the reproduction threshold (R 0 ) is tabulated in Table 6 Table 5. PRCC values for the parameters of the autonomous model (6) using total number of adult mosquitoes of type U , adult mosquitoes of type V , fourth instar larvae (L 4 ), pupae (P ), and R 0 as output. The top (most dominant) parameters that affect the dynamics of the model with respect to each of the six response function are highlighted in bold font. "Notation: a line (−) indicates the parameter is not in the expression for R 0 ".

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 T * 0 = E * , L * 1 , L * 2 , L * 3 , L * 4 , P * , V * , W * , U * = (0, 0, 0, 0, 0, 0, 0, 0, 0) and positive periodic solution(s) denoted by , U * (t) which satisfies the unique periodic system: Control measure Effect on population Effect on vectorial Environmental by model (1) dynamics of mosquitoes reproduction number interpretation R 0 Significant reduction Significant decrease in Significant decrease Personal protection in the value of α: the population size of in the value of R 0 against mosquito (probability of succ-adult mosquitoes of bite plays an imporessfully taking a blo-type U tant role in minimiod meal) zing the size of mosquito population in the community.

Significant reduction Significant decrease in Significant decrease
The removal of in the value of ψ U : the population size of in the value R 0 mosquito breeding (deposition rate of all three adult mosquito (egg laying) sites, female eggs) compartments such as removal of stagnant waters, is an effective control measure against the mosquito population.

Significant reduction Significant decrease in Significant decrease
The removal of in the value of ξ i the population size of in the value R 0 mosquito breeding (maturation rate of all three adult mosquito sites and use of female larvae) compartments larvicides are effective and significant incr-control measures ease of µ L against the mosquito (natural mortality population. rate of female larvae)

Significant increase
Significant decrease in Significant decrease The use of insecticides in the value of µ A : the population size of in the value of R 0 and insecticides treat-(natural mortality adult mosquitoes of ed bednets (ITNs) are rate of female adult type U important control mosquitoes) measures against the mosquito population. Table 6. Control measures obtained from the sensitivity analysis of the model (6).
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)  given, respectively, by where 0 denotes a zero matrix of order 3, and where τ * W , η * V , C i (t) (i = E, P, 1, . . . , 7) are as defined in (7). The linearized version of the model (1), at T * 0 , can be expressed as where x(t) = E(t), L 1 (t), L 2 (t), L 3 (t), L 4 (t), P (t), V (t), W (t), U (t) . Following [75], let Y (t, s), t ≥ s, be the evolution operator of the linear ω-periodic system dy dt = −V (t)y. Thus, for each s ∈ R, the associated 9 × 9 matrix Y (t, s) satisfies [75] dY (t, s) dt where I is the 9 × 9 identity matrix. Suppose that φ(s) (ω-periodic in s) is the initial distribution of new eggs. Thus, F (s)φ(s) is the rate of generation (hatching) of new eggs in the breeding habitat at time s [1,75]. Since t ≥ s, it follows that Y (t, s)F (s)φ(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 (φ(s)) introduced at a prior time s = t, is given by Let C ω be the ordered Banach space of all ω-periodic functions from R to R 9 , which is equipped with maximum norm and positive cone C + ω φ ∈ C ω : φ(t) ≥ 0, ∀ t ∈ R [75]. Define a linear operator L : The vectorial reproduction ratio of the model (22) (R 0T R ) is then given by the spectral radius of L, (i.e., R 0T R = ρ(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.2. The trivial equilibrium T * 0 of the non-autonomous model (1) is GAS in C([0], R 9 + ) whenever R 0T R < 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, R 0T R , 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 nonautonomous 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 ∂X 0 ) X = Ω, Proof. The proof is based on using the approach in [40,55]. Let u(t, φ) be the unique solution of model (1), with u(0, φ) = φ. Let Φ(t)ψ = u(t, ψ) and let P : X → X be the Poincaré map associated with (1) i.e., P(φ) = u(ω, φ) for all φ ∈ 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, P is point dissipative [40]. It then follows from Theorem 1.1.2 in [82] that P admits a global attractor in X.
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 (T ) is taken to beT = 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   (ii) including density dependence for the eggs oviposition process and larval mortality rates; (iii) 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: (i) The trivial equilibrium of the autonomous model (6) is globally-asymptotically stable whenever the associated vectorial reproduction number (R 0 ) is less than unity. For the case when the associated vectorial reproduction number (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., δ L = 0). This unique non-trivial equilibrium is shown to be globally-asymptotically stable under a certain condition (i.e., when 1 < R 0 < R C 0 ) (this equilibrium bifurcates into a limit cycle, via a Hopf bifurcation at R 0 = R C 0 ). (ii) 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 (α), the natural mortality rate of female adult mosquitoes (µ A ), the natural mortality rate of female larvae (µ L ), the deposition rate of female eggs (ψ U ) and the maturation rate of female larvae from Stage 1 to Stage 2 (ξ 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. (iii) 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 R 0T R ) 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 R 0T R > 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: (i) For mean monthly temperature and rainfall values in the range [10,40] Table 9. Monthly mean temperature (in 0 C) and rainfall (in mm) for Nairobi, Kenya [50]. where, with τ * W , η * V , B, C, C i (i = E, P, 1, . . . , 7) and D as defined in (7). (16).

Coefficients of Equation
C l + C E + C P > 0, C j + C E + C P   > 0, where τ * W , η * V , B, C, C i (i = E, P, 1, . . . , 7) and D are as defined in (7). Proof. Consider the model (6) with δ L = 0 and sgn(U * * − U (t)) = sgn(E * * − E(t)) for all t ≥ 0. Furthermore, let 1 < R 0 < R C 0 , so that the unique non-trivial equilibrium T 1 exists and is LAS (in line with Theorem 3.10). Consider, further, the following non-linear Lyapunov function of Goh-Volterra type: where, with τ * W , η * V , B, C, C i (i = E, P, 1, . . . , 7) and D as defined in (7). The following steady state relations will be used to simplify the derivative of the Lyapunov function K 2 .
The Lyapunov derivative of K 2 iṡ Substituting (24) and (25) into (26), and simplifying, giveṡ The first term of (27) is automatically negative in Ω \ {T 0 }, since sgn(U * * − U (t)) = sgn(E * * − E(t)) for all t ≥ 0. Furthermore, since the arithmetic mean exceeds the geometric mean, it follows that the second and third term of (27) are also negative. Hence,K 2 ≤ 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 ≥ 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 ≥ 0. Proof. Consider the non-autonomous model (1) with R 0T R < 1. Using the fact that (and noting that L(t) = 4 i=1 L i (t)), and, , for all t ≥ 0, it follows that the non-autonomous model (1), subject to the aforementioned assumptions, can be re-written as 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 It follows, from Lemma 2.1 in [80], that there exists a positive and bounded ωperiodic function, x(t) = Ē (t),L 1 (t),L 2 (t),L 3 (t),L 4 (t),P (t),V (t),W (t),Ū (t) , such that is a solution of the linearized system (28). Furthermore, it follows from Theorem 2.2 in [75] that R 0T R < 1 if and only if ρ φ F −V (τ ) < 1. Hence, θ is a negative constant. Thus, X(t) → 0 as t → ∞. 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], ((E, L 1 , L 2 , L 3 , L 4 , P, V, W, U )(0)) T ≤ Q * (Ē,L 1 ,L 2 ,L 3 ,L 4 ,P ,V ,W ,Ū )(0) T .