Modeling the variability of temperature on the population dynamics of Anopheles arabiensis

Abstract

Mosquitoes, notorious vectors of numerous diseases transmitted through bites, pose a significant public health threat, particularly in tropical and subtropical regions. Climate change may affect the risk of mosquito-borne diseases. This study investigated how temperature affects the population dynamics of all stages of Anopheles arabiensis mosquitoes, including eggs, larvae, pupae, and adults. We developed and analyzed a mathematical model that incorporates logistic growth and temperature. The well-posedness of the proposed model was proved. We demonstrated that if the non-autonomous basic reproduction number is less than unity, the disease-free equilibrium is locally asymptotically stable. Conversely, if it is greater than unity, there exists at least one positive periodic solution, as established by applying the comparison theorem and the theory of uniform persistence. Furthermore, the model parameters were fitted to real-world data conducted in the Tropical and Infectious Diseases Research Center (TIDRC) in Sekoru, Jimma University. The model presents the population dynamics of both immature and adult An. Arabiensis, which is similar to the observed experimental data obtained from TIDRC. Therefore, our model suggests that using the results obtained, it can be used to predict the impact of various intervention strategies on An. arabiensis distribution.

Introduction

Mosquitoes, notorious vectors of numerous diseases transmitted through bites, pose a significant public health threat, particularly in tropical and subtropical regions. These areas are vulnerable to mosquito-borne diseases such as malaria, yellow fever, Chikungunya, West Nile virus, dengue fever, Zika virus, and other arboviruses [1, 2]. About 17% of the world's infectious disease burden is caused by vector-borne diseases (VBDs), where insects or other vectors transmit infections from one host to another, accounting for approximately 17% of the global infectious disease burden [3].

Mosquitoes belonging to the genus Anopheles transmit Plasmodium spp. that cause malaria in humans. Of the 465 species identified within the genus Anopheles, about 70 are known to transmit malaria [4]. In Africa, three main Anopheles species: Anopheles gambiae sensu stricto (s.s.), An. arabiensis, and An. funestus are primarily responsible for the high rate of malaria transmission [4]. Malaria remains one of the deadliest vector-borne diseases, with an estimated 249 million malaria cases and nearly 608,000 deaths reported in 2022; of these 94% cases and 95% deaths occurred in Africa [5].

Ethiopia has documented 47 species of Anopheles mosquitoes [6], with An. arabiensis, a member species of the An. gambiae complex [7, 8], being the most efficient and widely distributed malaria vector. Anopheles arabiensis is the primary vector transmitting malaria in the country, while An. pharoensis, An. funestus, and An. nili are considered secondary vectors [9]. Complicating matters further, a recent study identified a new invasive species, An. stephensi, in Ethiopia [10], raising concerns about its potential to exacerbate malaria control efforts in the country [11].

Mosquitoes have four stages in their life cycle: egg, larva, pupa, and adult. The first three stages, egg, larva, and pupa, are considered immature and require water for survival [12]. The larval stage, the longest of the aquatic stages, develops through four molting phases ([2] and references therein). After the fourth molt, they enter the pupal stage. Pupae remain on the surface of water and undergo transformation into adult mosquitoes [13]. The entire life cycle, from egg to adult, can take anywhere from 7 to 20 days depending on factors like temperature, species, and the breeding site [14].

Temperature exerts a significant influence on the development rate, mortality rate, and reproduction of Anopheles mosquitoes at both immature and adult stages. As Bayoh and Lindsay (2003) demonstrated, this relationship is nonlinear [15]. For An. gambiae s.s., optimal development rates from one immature stage to the next occur around 28 °C. Adult development rates are highest between 28 and 32 °C, while adult emergence peaks between 22 and 26 °C. Temperatures below 18 °C or above 34 °C result in no adult emergence. Additionally, larval mortality is highest at the upper range of 30–32 °C, with death (rather than adult emergence) accounting for over 70% of the terminal events [16].

Understanding mosquito population dynamics is essential for gaining insights into their abundance and dispersal, which in turn informs the development of effective vector control strategies. A deep comprehension of these dynamics has significant implications for predicting and evaluating the impact of various vector control strategies [14, 17].

Mathematical modeling is a powerful tool for studying infectious disease epidemiology and transmission dynamics, particularly for vector-borne diseases like malaria. By translating complex biological processes such as the parasite life cycle and the mosquito’s life history into a structured framework, these models provide insights that help target interventions more effectively and predict the outcomes of elimination efforts [18,19,20,21]. When this framework is combined with real-world disease data, it enables us to make even more precise predictions about the future trajectory of the disease [22, 23].

A number of mathematical models have been developed to assess the impact of climate variables on mosquito populations. For example, Agusto et al. [24] studied how temperature affects malaria transmission, considering immunity development. Okuneye and Gumel [25] also examined the impact of temperature and rainfall on malaria in age-structured populations, focusing on mosquito dynamics. Yiga et al. [26] applied a human host-mosquito vector model to investigate malaria transmission dynamics. Abiodun et al. [27] developed a climate-based, ordinary differential equation model to analyze how temperature and water availability influence mosquito populations. Parham and Edwin [28] developed a model exploring the relationship between vector ecology and environmental factors. Okuneye et al. [29] studied the weather-driven malaria transmission model with gonotrophic and sporogonic cycles. Additionally, Abdelrazecy et al. [30] utilized a model to examine the potential effects of climate change on vector abundance and global malaria transmission.

Despite extensive research, a significant knowledge gap persists concerning the effects of temperature on the life cycle of An. arabiensis (egg, larva—including the four larval stages (instars), pupa, and adult). This information is crucial for developing accurate mathematical models to predict mosquito population dynamics. Building upon previous work [31], this study aims to develop a model that elucidates the relationship between temperature and An. arabiensis population dynamics.

Methods

Model formulation

The purpose of this section is to formulate a mathematical model of An. arabiensis that incorporates temperature variability. The complete metamorphosis of the Anopheles mosquito involves seven distinct developmental stages: egg, larva (with four instars), pupa, and female adult mosquito (Fig. 1).

Fig. 1

Stages of the Anopheles mosquito lifecycle (http://www.malaria.com/questions/mosquito-life-cycle/attachment/mos-lifecycle)

Mosquito egg hatching times differ greatly, ranging from days to months depending on the mosquito type, breeding habitats and egg-laying season. Larvae develop through four stages before pupation, with adult emergence occurring on the water surface [32]. Thus, our model is split into seven compartments at time t: compartments of eggs $\left(E\left(t\right)\right)$, four larvae stages $\left({L}_{i}\left(t\right),i=\text{1,2},\text{3,4}\right)$, pupa $\left(P\left(t\right)\right)$, and female adult mosquitoes $\left({A}_{m}\left(t\right)\right)$.

Furthermore, our model incorporates the effects of environmental temperature on Anopheles arabiensis developmental dynamics.. This is based on entomological data for Anopheles arabiensis obtained from Jimma University Tropical and Infectious Diseases Research Center (JU TIDRC).

Additionally, larvae and pupae require water and nutrients for their development. Employing the logistic growth concept from [31, 33], the carrying capacities for eggs, first-stage larvae, second-stage larvae, third-stage larvae, fourth-stage larvae, and pupae are respectively defined as:

${K}_{E}$,${K}_{{L}_{1}}$,${K}_{{L}_{2}}$,${K}_{{L}_{3}}$,${K}_{{L}_{4}}$, and ${K}_{p}$.

Following the logistic growth concept from [31, 33], we mathematically express this biological phenomenon in the model by introducing the availability of nutrients and the occupation by eggs, larvae, and pupae of the available breeder sites. Thus, the temperature-dependent per capita oviposition rate, the survival rate of eggs to first-stage larvae, the survival rate of first-stage larvae to second-stage larvae, the survival rate of second-stage larvae to third-stage larvae, the survival rate of third-stage larvae to fourth-stage larvae, and the survival rate of fourth-stage larvae to pupae are given by:

$$b\left(t\right)\left(1-\frac{E\left(t\right)}{{K}_{E}}\right){A}_{m}\left(t\right),{\beta }_{E}\left(t\right)\left(1-\frac{{L}_{1}\left(t\right)}{{K}_{{L}_{1}}}\right)E\left(t\right),{\beta }_{{L}_{1}}\left(t\right){\left(1-\frac{{L}_{2}\left(t\right)}{{K}_{{L}_{2}}}\right)L}_{1}\left(t\right),{\beta }_{{L}_{2}}\left(t\right){\left(1-\frac{{L}_{3}\left(t\right)}{{K}_{{L}_{3}}}\right)L}_{2}\left(t\right), {\beta }_{{L}_{3}}\left(t\right)\left(1-\frac{{L}_{4}\left(t\right)}{{K}_{{L}_{4}}}\right){L}_{3}\left(t\right), \text{and }{\beta }_{{L}_{4}}\left(t\right)\left(1-\frac{P\left(t\right)}{{K}_{P}}\right){L}_{4}\left(t\right),\text{respectively}.$$

A system of equations describing the mathematical model of Anopheles arabiensis is presented below. This model is based on the biological descriptions in Table 1, the flow diagram in Fig. 2, and the assumption that temperature varies with time.

$$\left(\text{t}\right).\left.\begin{array}{c}\frac{dE\left(t\right)}{dt}=b\left({T}_{a}\right){\left(1-\frac{E\left(t\right)}{{K}_{E}}\right)A}_{m}\left(t\right)-\left({{\mu }_{E}\left({T}_{w}\right)+\beta }_{E}\left({T}_{w}\right)\right)E\left(t\right), \\ \frac{d{L}_{1}\left(t\right)}{dt}={\beta }_{E}\left({T}_{w}\right)\left(1-\frac{{L}_{1}\left(t\right)}{{K}_{{L}_{1}}}\right)E\left(t\right)- \left({\beta }_{{L}_{1}}\left({T}_{W}\right)+{\mu }_{{L}_{1}}\left({T}_{w}\right)\right){L}_{1}\left(t\right),\\ \begin{array}{c}\frac{d{L}_{2}\left(t\right)}{dt}={\beta }_{{L}_{1}}\left({T}_{w}\right){\left(1-\frac{{L}_{2}\left(t\right)}{{K}_{{L}_{2}}}\right)L}_{1}\left(t\right)- \left({\beta }_{{L}_{2}}\left({T}_{w}\right)+{\mu }_{{L}_{2}}\left({T}_{w}\right)\right){L}_{2}\left(t\right),\\ \frac{d{L}_{3}\left(t\right)}{dt}={\beta }_{{L}_{2}}\left({T}_{w}\right){\left(1-\frac{{L}_{3}\left(t\right)}{{K}_{{L}_{3}}}\right)L}_{2}\left(t\right)- \left({\beta }_{{L}_{3}}\left({T}_{w}\right)+{\mu }_{{L}_{3}}\left({T}_{w}\right)\right){L}_{3}\left(t\right),\\ \begin{array}{c}\frac{d{L}_{4}\left(t\right)}{dt}={\beta }_{{L}_{3}}\left({T}_{w}\right)\left(1-\frac{{L}_{4}\left(t\right)}{{K}_{{L}_{3}}}\right){L}_{3}\left(t\right)- \left({\beta }_{{L}_{4}}\left({T}_{w}\right)+{\mu }_{{L}_{4}}\left({T}_{w}\right)\right){L}_{4}\left(t\right),\\ \frac{dP\left(t\right)}{dt}={\beta }_{{L}_{4}}\left({T}_{w}\right)\left(1-\frac{P\left(t\right)}{{K}_{P}}\right){L}_{4}\left(t\right)-\left({\beta }_{P}\left({T}_{w}\right)+{\mu }_{P}\left({T}_{w}\right)\right)P\left(t\right),\\ \begin{array}{c}\frac{d{A}_{m}\left(t\right)}{dt}={\beta }_{P}\left({T}_{w}\right)P\left(t\right)-{\mu }_{{A}_{m}}\left({T}_{w}\right){A}_{m}\left(t\right)\end{array},\end{array}\end{array}\end{array}\right\}$$

(1)

where ${T}_{a}$ and ${T}_{w}$ are air temperature and water temperature, respectively.

Table 1 Biological descriptions of the parameters in model 1

Fig. 2

Transfer diagram: The solid arrows represent the transition from one class to another. The dashed arrow represents the egg-laying of female adult mosquitoes. `

With initial conditions

$$E\left(0\right)\ge 0, {L}_{1}\left(0\right)\ge 0, {L}_{2}\left(0\right)\ge 0, {L}_{3}\left(0\right)\ge 0, {L}_{4}\left(0\right)\ge 0, P\left(0\right)\ge 0, {A}_{m}\left(0\right)\ge 0.$$

(2)

Theorem 1

System (1) with initial condition (2) has a unique, positive solution for all x ∈ ${\Omega }_{m}$= C $\left((0,{\mathbb{R}}_{+}^{7}]\right)$. Additionally, all solutions are uniformly bounded.

Proof

System (1) can be rewritten as:

$$\frac{dx}{dt}=f\left(t,x\left(t\right)\right),t\ge 0, x\left(0\right)={x}_{0}$$

(3)

where: $x(t)=\left(E\left(t\right), {L}_{1}\left(t\right), {L}_{2}\left(t\right), {L}_{3}\left(t\right), {L}_{4}\left(t\right), P\left(t\right), {A}_{m}\left(t\right),\right)\in {\mathbb{R}}_{+}^{7}$ and f is a continuous vector function defined

$$f\left(t,x\left(t\right)\right)=\left[\begin{array}{c}{f}_{1}\left(t,x\left(t\right)\right)\\ {f}_{2}\left(t,x\left(t\right)\right)\\ \begin{array}{c}{f}_{3}\left(t,x\left(t\right)\right)\\ {f}_{4}\left(t,x\left(t\right)\right)\\ \begin{array}{c}{f}_{5}\left(t,x\left(t\right)\right)\\ {f}_{6}\left(t,x\left(t\right)\right)\\ {f}_{7}\left(t,x\left(t\right)\right)\end{array}\end{array}\end{array}\right]=\left[\begin{array}{c}b\left({T}_{a}\right){\left(1-\frac{E\left(t\right)}{{K}_{E}}\right)A}_{m}\left(t\right)-\left({{\mu }_{E}\left({T}_{w}\right)+\beta }_{E}\left({T}_{w}\right)\right)\text{ E}\left(t\right)\\ {\beta }_{E}\left({T}_{w}\right)\left(1-\frac{{L}_{1}\left(t\right)}{{K}_{{L}_{1}}}\right)E\left(t\right)- \left({\beta }_{{L}_{1}}\left({T}_{w}\right)+{\mu }_{{L}_{1}}\left({T}_{w}\right)\right){L}_{1}\left(t\right)\\ \begin{array}{c}{\beta }_{{L}_{1}}{\left({T}_{w}\right)\left(1-\frac{{L}_{2}\left(t\right)}{{K}_{{L}_{2}}}\right)L}_{1}\left(t\right)- \left({\beta }_{{L}_{2}}\left({T}_{w}\right)+{\mu }_{{L}_{2}}\left({T}_{w}\right)\right){L}_{2}\left(t\right)\\ {\beta }_{{L}_{2}}\left({T}_{w}\right){\left(1-\frac{{L}_{3}\left(t\right)}{{K}_{{L}_{3}}}\right)L}_{2}\left(t\right)- \left({\beta }_{{L}_{3}}\left({T}_{w}\right)+{\mu }_{{L}_{3}}\left({T}_{w}\right)\right){L}_{3}\left(t\right)\\ \begin{array}{c}{\beta }_{{L}_{3}}\left({T}_{w}\right)\left(1-\frac{{L}_{4}\left(t\right)}{{K}_{{L}_{4}}}\right){L}_{3}\left(t\right)- \left({\beta }_{{L}_{4}}\left({T}_{w}\right)+{\mu }_{{L}_{4}}\left({T}_{w}\right)\right){L}_{4}\left(t\right)\\ {\beta }_{{L}_{4}}{\left({T}_{w}\right)\left(1-\frac{P\left(t\right)}{{K}_{P}}\right)L}_{4}\left(t\right)-\left({\beta }_{P}\left({T}_{w}\right)+{\mu }_{P}\left({T}_{w}\right)\right)P\left(t\right)\\ \begin{array}{c}{\beta }_{P}\left({T}_{w}\right)P\left(t\right)-{\mu }_{{A}_{m}}\left({T}_{w}\right){A}_{m}\left(t\right).\end{array}\end{array}\end{array}\end{array}\right]$$

(4)

and

$${x}_{0}={\left(\begin{array}{c}E\left(0\right), {L}_{1}\left(0\right), {L}_{2}\left(0\right), {L}_{3}\left(0\right), {L}_{4}\left(0\right), P\left(0\right), {A}_{m}\left(0\right),\end{array}\right)}^{T}\in {R}_{+}^{7}.$$

(5)

Thus, for all x ∈ ${\Omega }_{m}$ the function $f\left(t,x\left(t\right)\right)$ is ${C}^{1}$ continuous and Lipschitzian (with respect to $x$ in each compact set in $({\mathbb{R}}\times {\Omega }_{m}$) [25]. Hence, there is a unique solution of system (1) through $\left(0,x\right)$.

To show the positivity of the solution of the $x(t)=\left(E\left(t\right), {L}_{1}\left(t\right), {L}_{2}\left(t\right), {L}_{3}\left(t\right), {L}_{4}\left(t\right), P\left(t\right), {A}_{m}\left(t\right)\right)$ subject to the initial condition (2) for all $t\ge 0$,we consider T is a number defined as.

$T=sup\left\{t>0:E\left(t\right)>0, {L}_{1}\left(t\right)>0, {L}_{2}\left(t\right)>0, {L}_{3}\left(t\right)>0, {L}_{4}\left(t\right)>0, P\left(t\right)>0, {A}_{m}\left(t\right)>0\right\}$. one can easily conclude that $T>0$. if we again consider the case either $t = T = \infty$ or $T < \infty .$

In earlier case, all solutions are obviously positive. In later case, all solutions are also positive. To prove this, we take the contradiction that at least one of $E\left(T\right),{L}_{1}\left(T\right), {L}_{2}\left(T\right), {L}_{3}\left(T\right), {L}_{4}\left(T\right), P\left(T\right),$ and ${A}_{m}\left(T\right)$ is zero. Let us suppose that $E\left(T\right)=0$.

We can now write the first equation in the model system (1) as follows:

$\frac{dE\left(t\right)}{dt}\ge -\left({\mu }_{E}+{\beta }_{E}\right)E\left(t\right)$,this yields.

$$E\left(T\right)\ge E\left(0\right)exp\left(-{\int }_{0}^{T}\left({\mu }_{E}\left({T}_{w}\right)+{\beta }_{E}\left({T}_{w}\right)\right)ds\right)>0.$$

Thus, $E\left(T\right)>0$ implies that $E\left(t\right)>0$ for any $t>0$. Which contradict our assumption that $E\left(T\right)=0$. Thus, $E\left(t\right)>0$ for all $t>0$.

Similarly, the rest of the Eqs. 2–9 in Eq. (1) can be rewritten as:

$$\frac{d{L}_{1}\left(t\right)}{dt}\ge - \left({\beta }_{{L}_{1}}\left({T}_{w}\right)+{\mu }_{{L}_{1}}\left({T}_{w}\right)\right){L}_{1}\left(t\right),$$

(6)

$$\frac{d{L}_{2}\left(t\right)}{dt}\ge - \left({\beta }_{{L}_{2}}\left({T}_{w}\right)+{\mu }_{{L}_{2}}\left({T}_{w}\right)\right){L}_{2}\left(t\right),$$

(7)

$$\frac{d{L}_{3}\left(t\right)}{dt}\ge - \left({\beta }_{{L}_{3}}\left({T}_{w}\right)+{\mu }_{{L}_{3}}\left({T}_{w}\right)\right){L}_{3}\left(t\right),$$

(8)

$$\frac{d{L}_{4}\left(t\right)}{dt}\ge - \left({\beta }_{{L}_{4}}\left({T}_{w}\right)+{\mu }_{{L}_{4}}\left({T}_{w}\right)\right){L}_{4}\left(t\right),$$

(9)

$$\frac{dP\left(t\right)}{dt}\ge -\left({\beta }_{P}\left({T}_{w}\right)+{\mu }_{P}\left({T}_{w}\right)\right)P\left(t\right),$$

(10)

$$\frac{d{A}_{m}\left(t\right)}{dt}\ge -{\mu }_{{A}_{m}}\left({T}_{w}\right){A}_{m}\left(t\right),$$

(11)

this leads to following solutions,

$${L}_{1}\left(T\right)\ge {L}_{1}\left(0\right)exp\left(- {\int }_{0}^{t}\left({\beta }_{{L}_{1}}\left({T}_{w}\right)+{\mu }_{{L}_{1}}\left({T}_{w}\right)\right)ds\right)>0,$$

(12)

$${L}_{2}\left(T\right)\ge {L}_{2}\left(0\right)exp\left(- {\int }_{0}^{t}\left({\beta }_{{L}_{2}}\left({T}_{w}\right)+{\mu }_{{L}_{2}}\left({T}_{w}\right)\right)ds\right)>0,$$

(13)

$${L}_{3}\left(T\right)\ge {L}_{3}\left(0\right)exp\left(- {\int }_{0}^{t}\left({\beta }_{{L}_{3}}\left({T}_{w}\right)+{\mu }_{{L}_{3}}\left({T}_{w}\right)\right)ds\right)>0,$$

(14)

$${L}_{4}\left(T\right)\ge {L}_{4}\left(0\right)exp\left(- {\int }_{0}^{t}\left({\beta }_{{L}_{4}}\left({T}_{w}\right)+{\mu }_{{L}_{4}}\left({T}_{w}\right)\right)ds\right)>0,$$

(15)

$$P\left(T\right)\ge P\left(0\right)exp\left(-{\int }_{0}^{t}\left({\beta }_{P}\left({T}_{w}\right)+{\mu }_{P}\left({T}_{w}\right)\right)ds\right)>0,$$

(16)

$${A}_{m}\left(T\right)\ge {A}_{m}\left(0\right)exp\left(-{\int }_{0}^{t}\left({\mu }_{{A}_{m}}\left({T}_{w}\right)\right)ds\right)>0.$$

(17)

As we can see from Eqs. (12)–(17), ${L}_{1}\left(T\right)>0, {L}_{2}\left(T\right)>0, {L}_{3}\left(T\right)>0, {L}_{4}\left(T\right)>0, P\left(T\right)>0,$ and ${A}_{m}\left(T\right)>0$, and this contradicts the assumption ${L}_{1}\left(T\right)=0, {L}_{2}\left(T\right)=0, {L}_{3}\left(T\right)=0, {L}_{4}\left(T\right)=0, P\left(T\right)=0,$ and ${A}_{m}\left(T\right)=0$.

Hence, we conclude that ${L}_{1}\left(t\right)>0, {L}_{2}\left(t\right)>0, {L}_{3}\left(t\right)>0, {L}_{4}\left(t\right)>0, P\left(t\right)>0, {A}_{m}\left(t\right)>0$ for all $t>0$.

Consequently, we show that the model system (1) is positively invariant and attracting within ${\mathbb{R}}_{+}^{7}$.

Furthermore,

$${\Omega }_{m}=\left\{\left(E, {L}_{1}, {L}_{2}, {L}_{3}, {L}_{4}, P, {A}_{m}\right)\in {R}_{+}^{7}:E\le {K}_{E},{L}_{1}\le {K}_{{L}_{1}},{L}_{2}\le {K}_{{L}_{2}},{L}_{3}\le {K}_{{L}_{3}},{L}_{4}\le {K}_{{L}_{4}},P\le {K}_{P},{A}_{m}\le \frac{{\beta }_{P}{K}_{P}}{{\mu }_{{A}_{m}}}\right\}.$$

Therefore, we conclude that the proposed model system (1) is mathematically well-posed in the domain ${\Omega }_{m}$ and is epidemiologically sound.

Model analysis

In this section, we present the mosquito-free equilibrium point, the vector reproduction number, stability of the mosquito-free equilibrium, and the existence of a positive periodic solution.

Mosquito-free equilibrium point

Mosquito-free equilibrium point of the model system (1) is obtained by setting the right side of each equation in (1) to zero. This results in the Mosquito-free equilibrium point $\left({X}_{0}\right)$, where ${X}_{0}=\left({0,0},{0,0},{0,0},0\right)$.

The vector reproduction number

To find the vector reproduction number of the non-autonomous model system (1), we employ the method used in Wang and Zhao [34] which extended the framework of Van den Driessche and James Watmough [35].

$$\left.\begin{array}{c}\frac{dE\left(t\right)}{dt}=b\left({T}_{a}\right){\left(1-\frac{E\left(t\right)}{{K}_{E}}\right)A}_{m}\left(t\right)-\left({{\mu }_{E}\left({T}_{w}\right)+\beta }_{E}\left({T}_{w}\right)\right)E(t) , \\ \frac{d{L}_{1}\left(t\right)}{dt}={\beta }_{E}\left({T}_{w}\right)\left(1-\frac{{L}_{1}\left(t\right)}{{K}_{{L}_{1}}}\right)E\left(t\right)- \left({\beta }_{{L}_{1}}\left({T}_{w}\right)+{\mu }_{{L}_{1}}\left({T}_{w}\right)\right){L}_{1}\left(t\right),\\ \begin{array}{c}\frac{d{L}_{2}\left(t\right)}{dt}={\beta }_{{L}_{1}}{\left({T}_{w}\right)\left(1-\frac{{L}_{2}\left(t\right)}{{K}_{{L}_{2}}}\right)L}_{1}\left(t\right)- \left({\beta }_{{L}_{2}}\left({T}_{w}\right)+{\mu }_{{L}_{2}}\left({T}_{w}\right)\right){L}_{2}\left(t\right),\\ \frac{d{L}_{3}\left(t\right)}{dt}={\beta }_{{L}_{2}}\left({T}_{w}\right){\left(1-\frac{{L}_{3}\left(t\right)}{{K}_{{L}_{3}}}\right)L}_{2}\left(t\right)- \left({\beta }_{{L}_{3}}\left({T}_{w}\right)+{\mu }_{{L}_{3}}\left({T}_{w}\right)\right){L}_{3}\left(t\right),\\ \begin{array}{c}\frac{d{L}_{4}\left(t\right)}{dt}={\beta }_{{L}_{3}}\left({T}_{w}\right)\left(1-\frac{{L}_{4}\left(t\right)}{{K}_{{L}_{4}}}\right){L}_{3}\left(t\right)- \left({\beta }_{{L}_{4}}\left({T}_{w}\right)+{\mu }_{{L}_{4}}\left({T}_{w}\right)\right){L}_{4}\left(t\right),\\ \frac{dP\left(t\right)}{dt}={\beta }_{{L}_{4}}{\left({T}_{w}\right)\left(1-\frac{P\left(t\right)}{{K}_{P}}\right)L}_{4}\left(t\right)-\left({\beta }_{P}\left({T}_{w}\right)+{\mu }_{P}\left({T}_{w}\right)\right)P\left(t\right),\\ \begin{array}{c}\frac{d{A}_{m}\left(t\right)}{dt}={\beta }_{P}\left({T}_{w}\right)P\left(t\right)-{\mu }_{{A}_{m}}\left({T}_{w}\right){A}_{m}\left(t\right).\end{array}\end{array}\end{array}\end{array}\right\}$$

(18)

Equation (18) can be written as

$${Z}^{{\prime}\left(t\right)}=\left(F\left(t\right)-V\left(t\right)\right)Z\left(t\right),$$

(19)

where

$Z\left(t\right)={\left(E\left(t\right), {L}_{1}\left(t\right), {L}_{2}\left(t\right), {L}_{3}\left(t\right), {L}_{4}\left(t\right), P\left(t\right), {A}_{m}\left(t\right)\right)}^{T}$, (F(t)), the matrix of new generated mosquitoes, and (V(t)), the matrix of transfer in and out of the compartments.

Linearizing system (18) at the mosquito-free equilibrium ${X}_{0} = (0, 0, 0, {0,0},{0,0})$, We obtain, respectively, the following:

$$F\left(t\right)=\left[\begin{array}{lllllll}0& 0& 0& 0& 0& 0& b\left(t\right)\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right]$$

and

$$V\left(t\right)=\left[\begin{array}{lllllll}{C}_{1}& 0& 0& 0& 0& 0& 0\\ -{\beta }_{E}& {C}_{2}& 0& 0& 0& 0& 0\\ 0& -{\beta }_{{L}_{1}}& {C}_{3}& 0& 0& 0& 0\\ 0& 0& -{\beta }_{{L}_{2}}& {C}_{4}& 0& 0& 0\\ 0& 0& 0& -{\beta }_{{L}_{3}}& {C}_{5}& 0& 0\\ 0& 0& 0& 0& -{\beta }_{{L}_{4}}& \left({\beta }_{P}+{\sigma }_{p}\left(t\right)\right)& 0\\ 0& 0& 0& 0& 0& -{\beta }_{P}& {\sigma }_{V}\left(t\right)\end{array}\right].$$

$${C}_{1}={\beta }_{E}\left({T}_{w}\right){C}_{2}=\left({\mu }_{{L}_{1}}\left({T}_{w}\right)+{\beta }_{{L}_{1}}\left({T}_{w}\right)\right){C}_{3}=\left({\mu }_{{L}_{3}}\left({T}_{w}\right)+{\beta }_{{L}_{3}}\left({T}_{w}\right)\right){C}_{4}=\left({\mu }_{{L}_{4}}\left({T}_{w}\right)+{\beta }_{{L}_{4}}\left({T}_{w}\right)\right){C}_{5}={\beta }_{P}\left({T}_{w}\right),$$

Following [34], let ${\phi }_{M}$ be the monodromy matirx of the linear $\omega$-periodic system

$$Z\dot{\vphantom{a}} = M (t)Z$$

(20)

Furthermore, let $\rho \left({\phi }_{M}\left(\omega \right)\right)$ be the spectra radius of ${\phi }_{M}\left(\omega \right)$ and for all $t \ge s$, let $Y (t, s)$ be the evolution operator of the linear $\omega$-periodic system

$$y\dot{\vphantom{a}} = -V \left(t\right)y.$$

(21)

For each $s \in R$, the $7 \times 7$ matrix $Y (t, s)$ satisfies the equation

$$Y\dot{\vphantom{a}} \left(t, s\right)= -V \left(t\right)Y (t, s), \forall t \ge s, Y (s, s) = I,$$

(22)

where I is the 7 × 7 identity matrix.

Let ${C}_{\omega }$ be the ordered Banach space of all $\omega$-periodic functions from ${\mathbb{R}}$ to ${\mathbb{R}}_{+}^{7}$ which is equipped with the maximum norm $\Vert .\Vert$ and the positive cone

$${C}_{\omega }^{+} = \left\{\varphi \in {C}_{W} : \varphi (t) \ge 0, \forall t\in {\mathbb{R}}\right\}.$$

(23)

Suppose that $\varphi (s)\in {C}_{\omega }$ is the initial distribution of juvenile and adults mosquitoes. Thus $F(s)\varphi (s) \in {C}_{\omega }$ is the distribution of new juvenile mosquitoes in the breeding habits produced by the adult ones which were introduced at time s. Hence, for any $t\ge s$, $Y (t, s) F(s)\varphi (s)$ is the distribution of those mosquitoes which were newly born into the juvenile mosquito compartment at time s and remain alive.

Thus,

$$\psi \left(t\right)={\int }_{-\infty }^{t}Y (t, s) F(s)\varphi (s)ds={\int }_{0}^{\infty }Y (t, t-b) F(t-b)\varphi (t-b)db,b\in \left[0\right.,\left.\infty \right)$$

(24)

is the distribution of new eggs at time t, hatched by all female adult mosquitoes ϕ(s) introduced at the previous time s.

Let $L : {C}_{\omega } \to {C}_{\omega }$ be the linear operator defined by.

$$\left(L\varphi \right)\left(t\right)= {\int }_{0}^{\infty }Y \left(t, t - a\right)F\left(t - a\right)\varphi \left(t - a\right)da, \forall t \in R, \varphi \in {C}_{\omega }.$$

(25)

Then, the vector reproduction ratio is ${R}_{v} := \rho (L),$ the spectral radius of L.

Thus, in order to get the approximate numerical value of ${R}_{v}$ we follow the steps Wang and Zhao [34] follows that is bringing in mind the method in as in [35]. Wang and Zhao [34] described ${\psi }_{\frac{F}{\lambda }-V}\left(t,\lambda \right)$ as the fundamental matrix of the linear ω-periodic equation

$$\frac{dZ}{dt}=\left(\frac{F\left(t\right)}{\lambda }-V\left(t\right)\right)Z,t\in R$$

(26)

with parameter $\lambda \in (0,\infty ).$

As in [34] $Z(t,s, \lambda ), t \ge s, s \in R$, is the evolution operator of the system (26) on ${R}_{v}$,and ${\psi }_{F-V}\left(t\right)=Z(t,0, 1)$, $\forall t \ge 0$.

Stability of mosquito-free equilibrium

Theorem 2 ([34, Theorem 2.2])

Assume that (A1)–(A7) as in [34] hold. Then the following statements are valid:

• a. ${R}_{v}=1 \iff \rho \left({\phi }_{F-V}\left(\omega \right)\right)=1$.

• b. ${R}_{v}<1\iff \rho \left({\phi }_{F-V}\left(\omega \right)\right)<1$.

• c. ${R}_{v}>1 \iff \rho \left({\phi }_{F-V}\left(\omega \right)\right)>1$.

Thus, ${X}_{0}\left(t\right)$ is asymptotically stable if ${R}_{v}<1$, and unstable if ${R}_{v}>1$.

Theorem 3

The mosquito free equilibrium point ${X}_{0}=\left(\text{0,0},\text{0,0},\text{0,0},0\right)$ is locally asymptotically stable if ${R}_{v}<1$.

Proof

The proof the theorem 3 follows directly from theorem 2.

Theorem 4

The mosquito free equilibrium point ${X}_{0}=\left(\text{0,0},\text{0,0},\text{0,0},0\right)$ is globally asymptotically stable if ${R}_{v}<1$.

Proof

We follow the steps as in [31] that is for all $t\ge 0$ and using the assumptions $\left(1-\frac{E\left(t\right)}{{K}_{E}}\right)\le 1$,$\left(1-\frac{{L}_{1}\left(t\right)}{{K}_{{L}_{1}}}\right)\le 1$, $\left(1-\frac{{L}_{2}\left(t\right)}{{K}_{{L}_{2}}}\right)\le 1$, $\left(1-\frac{{L}_{3}\left(t\right)}{{K}_{{L}_{3}}}\right)\le 1$, $\left(1-\frac{{L}_{4}\left(t\right)}{{K}_{{L}_{4}}}\right)\le 1$, $\left(1-\frac{P\left(t\right)}{{K}_{P}}\right)\le 1$ since $E\left(t\right)\le {K}_{E}$, ${L}_{1}\le {K}_{{L}_{1}}$, ${L}_{2}\le {K}_{{L}_{2}}$,${L}_{3}\le {K}_{{L}_{3}}$,${L}_{4}\le {K}_{{L}_{4}}$, $P\le {K}_{P}$, the model system (1) can be rewritten in the form:

$$\left.\begin{array}{c}\frac{d\widetilde{E}\left(t\right)}{dt}\le b\left({T}_{a}\right){A}_{m}\left(t\right)-\left({{\mu }_{E}\left({T}_{w}\right)+\beta }_{E}\left({T}_{w}\right)\right)\widetilde{E}\left(t\right), \\ \frac{d\widetilde{{L}_{1}}\left(t\right)}{dt}\le {\beta }_{E}\left({T}_{w}\right)E\left(t\right)- \left({\beta }_{{L}_{1}}\left({T}_{w}\right)+{\mu }_{{L}_{1}}\left({T}_{w}\right)\right)\widetilde{{L}_{1}}\left(t\right),\\ \begin{array}{c}\frac{d\widetilde{{L}_{2}}\left(t\right)}{dt}\le {\beta }_{{L}_{1}}\left({T}_{w}\right){L}_{1}\left(t\right)- \left({\beta }_{{L}_{2}}\left({T}_{w}\right)+{\mu }_{{L}_{2}}\left({T}_{w}\right)\right)\widetilde{{L}_{2}}\left(t\right),\\ \frac{d\widetilde{{L}_{3}}\left(t\right)}{dt}\le {\beta }_{{L}_{2}}\left({T}_{w}\right){L}_{2}\left(t\right)- \left({\beta }_{{L}_{3}}\left({T}_{w}\right)+{\mu }_{{L}_{3}}\left({T}_{w}\right)\right)\widetilde{{L}_{3}}\left(t\right),\\ \begin{array}{c}\frac{d\widetilde{{L}_{4}}\left(t\right)}{dt}\le {\beta }_{{L}_{3}}\left({T}_{w}\right){L}_{3}\left(t\right)- \left({\beta }_{{L}_{4}}\left({T}_{w}\right)+{\mu }_{{L}_{4}}\left({T}_{w}\right)\right)\widetilde{{L}_{4}}\left(t\right),\\ \frac{d\widetilde{P}\left(t\right)}{dt}\le {\beta }_{{L}_{4}}\left({T}_{w}\right){L}_{4}\left(t\right)-\left({\beta }_{P}\left({T}_{w}\right)+{\mu }_{P}\left({T}_{w}\right)\right)\widetilde{P}\left(t\right),\\ \begin{array}{c}\frac{d\widetilde{{A}_{m}}\left(t\right)}{dt}={\beta }_{P}\left({T}_{w}\right)P\left(t\right)-{\mu }_{{A}_{m}}\left({T}_{w}\right)\widetilde{{A}_{m}}\left(t\right)\end{array}.\end{array}\end{array}\end{array}\right\}$$

(27)

Further, following the work of [34], the equations in (27), with equalities used in place of the inequalities, can be rewritten in terms of the next generation matrices ${F}_{1}\left(t\right)$ and ${V}_{1}\left(t\right)$, as follows:

$${\widetilde{K}}^{{\prime}\left(t\right)}=\left({F}_{1}\left(t\right)-{V}_{1}\left(t\right)\right)\widetilde{K}\left(t\right).$$

(28)

It follows, from theorem 2b, that there exists a positive τ-periodic function,

$$\widetilde{K}\left(t\right)={\left(\widetilde{E}\left(t\right),\widetilde{{L}_{1}}\left(t\right),\widetilde{{L}_{2}}\left(t\right),\widetilde{{L}_{3}}\left(t\right),\widetilde{{L}_{4}}\left(t\right),\widetilde{P}\left(t\right),\widetilde{{A}_{m}}\left(t\right)\right)}^{T},$$

(29)

such that $\widetilde{K}\left(t\right)=v\left(t\right){e}^{rt}$ with $r =\frac{1}{\omega }\text{ln}\rho \left({\phi }_{{F}_{1}-{V}_{1}}\left(\omega \right)\right)$ the solution of the linearized system(28).

Furthermore, the assumption ${R}_{v}<1$ implies that $\rho \left({\phi }_{{F}_{1}-{V}_{1}}\left(\omega \right)\right)<1$ (by theorem 2b). Hence, $r < 0$.Thus, $\widetilde{K}\left(t\right)\to 0$ as $t \to \infty .$ By applying the comparison theorem [36] on system (28), we get

$$\underset{t\to \infty }{\text{lim}}\left(E\left(t\right), {L}_{1}\left(t\right), {L}_{2}\left(t\right), {L}_{3}\left(t\right), {L}_{4}\left(t\right), P\left(t\right), {A}_{m}\left(t\right)\right)=\left(\text{0,0},\text{0,0},\text{0,0},0\right).$$

(30)

Finally, we can conclude that the mosquito-free equilibrium ${X}_{0}$ is globally attractive, which completes the proof.

Definition 1. ([37])

A function $f : X \to X$ is said to be uniformly persistent with respect to $({\text{X}}_{0}, \partial {\text{X}}_{0})$ if there exists $\upeta > 0$ such that $\underset{t\to \infty }{\text{lim}}\text{inf}\left({f}^{n}\left(x\right),\partial {X}_{0}\right)\ge \eta$ for all$x \in {X}_{0}$. If “inf” in this inequality is replaced with “sup”, then f is said to be weakly uniformly persistent with respect to$({\text{X}}_{0}, \partial {\text{X}}_{0})$.

Existence of positive periodic solution.

This section presents the existence of a positive periodic solution of the model system (1) by employing the uniform persistence theorem (see [2] and the references therein). Employing the notions and methods used in Lou and Zhao [38], we define $\left(X,{X}_{0},\partial X\right)$ as follows:

$$X:={R}_{+}^{7},$$

$${X}_{0}=\left\{\left(E\left(t\right), {L}_{1}\left(t\right), {L}_{2}\left(t\right), {L}_{3}\left(t\right), {L}_{4}\left(t\right), P\left(t\right), {A}_{m}\left(t\right)\right)\in X:E>0, {L}_{1}>0, {L}_{2}>0, {L}_{3}>0, {L}_{4}>0, P>0, {A}_{m}>0\right\},$$

$${\partial X}_{0}:=\left\{\left(E\left(t\right), {L}_{1}\left(t\right), {L}_{2}\left(t\right), {L}_{3}\left(t\right), {L}_{4}\left(t\right), P\left(t\right), {A}_{m}\left(t\right)\right)\in X:E {L}_{1}{L}_{2} {L}_{3} {L}_{4}P{A}_{m}=0\right\}.$$

Theorem 6

If ${R}_{v} > 1$, there exists $\beta >0$ such that any solution $\left(E\left(t\right), {L}_{1}\left(t\right), {L}_{2}\left(t\right), {L}_{3}\left(t\right), {L}_{4}\left(t\right), P\left(t\right), {A}_{m}\left(t\right)\right)$ with initial condition $\varphi \in {X}_{0}$ satisfies $\underset{t\to \infty }{\text{lim}}\text{inf}E\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{L}_{1}\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{L}_{2}\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{L}_{3}\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{L}_{4}\left(t\right)\ge \beta$,,$\underset{t\to \infty }{\text{lim}}\text{inf}P\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{A}_{m}\left(t\right)\ge \beta$ and the system (1) has at least one positive periodic solution.

Proof. The proof of theorem 6 is based on using the techniques in [2, 32]. Let $u\left(t,\varphi \right)$ be the unique solution of model (1), with $u\left(0,\varphi \right)=\varphi$. Let $\varphi \left(t\right)\psi =u\left(t,\psi \right)$ and $P:X\to X$ be the Poincare map associated with (1) that is.

$$P\left(\varphi \right)=\psi \left(\omega \right)\varphi =u\left(\omega ,\varphi \right),\forall \varphi \in X.$$

$${P}^{m}\left(\varphi \right)=\phi \left(m\omega \right)\varphi =u\left(m\omega ,\varphi \right),\forall m\in X.$$

Obviously,$u\left(t,\psi \right)\in {X}_{0}$,$\forall t\ge 0$. This implies that both ${X}_{0}$ and ${\partial X}_{0}$ are positively invariant. Since solutions are bounded the Poincare map (P) is point dissipative [38]. This, along with Theorem 1.1.2 in [25], guarantees that P admits a global attractor for P in X.

Thus, it remains to show that the model system is a uniform persistence with respect to $\left({X}_{0},\partial {X}_{0}\right)$. To do this, we first define ${M}_{\partial }$ as:

$${M}_{\partial }=\left\{\varphi \in \partial {X}_{0} : {P}^{n}\left(\varphi \right) \in \partial {X}_{0},\forall n\ge 0\right\}.$$

(31)

Then we first show that: ${M}_{\partial }=\left\{\left(\text{0,0},\text{0,0},\text{0,0},0\right)\right\}$. To achieve this, we need to demonstrate two inclusions: (i) $\left\{\left(\text{0,0},\text{0,0},\text{0,0},0\right)\right\}\subset {M}_{\partial }$ and (ii)${M}_{\partial }\subset \left\{\left(\text{0,0},\text{0,0},\text{0,0},0\right)\right\}$. The first inclusion is self-evident. For the second inclusion, we must show that for any initial condition $\varphi \in \partial {X}_{0}$, $E\left(k\omega \right){l}_{1}\left(k\omega \right){l}_{2}\left(k\omega \right){l}_{3}\left(k\omega \right){l}_{4}\left(k\omega \right) P\left(k\omega \right){A}_{m}\left(k\omega \right)=0$ for all k ≥ 0. We employ proof by contradiction. Suppose $\varphi \in \partial {X}_{0}$ and there exists a positive constant ${k}_{1}$ such that ${\left(E\left({k}_{1}\omega \right),{l}_{1}\left({k}_{1}\omega \right),{l}_{2}\left({k}_{1}\omega \right),{l}_{3}\left({k}_{1}\omega \right){,l}_{4}\left({k}_{1}\omega \right), P\left({k}_{1}\omega \right),{A}_{m}\left({k}_{1}\omega \right)\right)}^{T}>0$, for all ${k}_{1}$. This assumption leads to a contradiction, proving the desired result.

Rewriting $\frac{dE\left(t\right)}{dt}$ from the system of model (1) yields Eq. (32):

$$\frac{dE\left(t\right)}{dt}+\left(\frac{b\left({T}_{a}\right){A}_{m}\left(t\right)}{{K}_{E}}+{\mu }_{E}\left({T}_{w}\right)+{\beta }_{E}\left({T}_{w}\right)\right)E\left(t\right)=b\left({T}_{a}\right){A}_{m}\left(t\right)$$

(32)

Applying the method of constant variation to Eq. (32) results in Eq. (33):

$$E\left(t\right)=exp\left(-{\int }_{0}^{t}\left({K}_{1}\left(s\right)ds\right)\right)\left[E\left(0\right)+{\int }_{0}^{t}b\left(s\right){A}_{m}\left(s\right)exp\left({\int }_{0}^{s}{K}_{1}\left(\tau \right)d\tau \right)ds\right],$$

(33)

where ${K}_{1}\left(t\right)$ is defined as:

$${K}_{1}\left(t\right)=\frac{b\left({T}_{a}\right){A}_{m}\left(t\right)}{{K}_{E}}+{\mu }_{E}\left({T}_{w}\right)+{\beta }_{E}\left({T}_{w}\right).$$

(34)

Similarly, applying the method of constant variation Eqs. (2)–(7) in the system model (1) are equivalent to:

$${L}_{1}\left(t\right)=exp\left(-{\int }_{0}^{t}\left({K}_{2}\left(s\right)ds\right)\right)\left[{L}_{1}\left(0\right)+{\int }_{0}^{t}{\beta }_{E}\left(s\right)E\left(s\right)exp\left({\int }_{0}^{s}{K}_{2}\left(\tau \right)d\tau \right)ds\right]$$

(35)

$${L}_{2}\left(t\right)=exp\left(-{\int }_{0}^{t}\left({K}_{3}\left(s\right)ds\right)\right)\left[{L}_{2}\left(0\right)+{\int }_{0}^{t}{\beta }_{{L}_{1}}\left(s\right){L}_{1}\left(s\right)exp\left({\int }_{0}^{s}{K}_{3}\left(\tau \right)d\tau \right)ds\right]$$

(36)

$${L}_{3}\left(t\right)=exp\left(-{\int }_{0}^{t}\left({K}_{4}\left(s\right)ds\right)\right)\left[{L}_{3}\left(0\right)+{\int }_{0}^{t}{\beta }_{{L}_{2}}\left(s\right){L}_{2}\left(s\right)exp\left({\int }_{0}^{s}{K}_{4}\left(\tau \right)d\tau \right)ds\right]$$

(37)

$${L}_{4}\left(t\right)=exp\left(-{\int }_{0}^{t}\left({K}_{5}\left(s\right)ds\right)\right)\left[{L}_{4}\left(0\right)+{\int }_{0}^{t}{\beta }_{{L}_{3}}\left(s\right){L}_{3}\left(s\right)exp\left({\int }_{0}^{s}{K}_{5}\left(\tau \right)d\tau \right)ds\right]$$

(38)

$$P\left(t\right)=exp\left(-{\int }_{0}^{t}\left({K}_{6}\left(s\right)ds\right)\right)\left[P\left(0\right)+{\int }_{0}^{t}{\beta }_{{L}_{4}}\left(s\right){L}_{4}\left(s\right)exp\left({\int }_{0}^{s}{K}_{6}\left(\tau \right)d\tau \right)ds\right]$$

(39)

$${A}_{m}\left(t\right)=exp\left(-{\int }_{0}^{t}{\sigma }_{m}\left(s\right)ds\right)\left[{A}_{m}\left(0\right)+{\int }_{0}^{t}{\beta }_{P}\left(s\right)P\left(s\right)exp\left({\int }_{0}^{s}{\mu }_{{A}_{m}}\left(\tau \right)d\tau \right)\right]ds$$

(40)

where ${K}_{2}\left(t\right)$, ${K}_{3}\left(t\right)$ , ${K}_{4}\left(t\right)$, ${K}_{5}\left(t\right)$, ${K}_{6}\left(t\right)$ is defined as:

$${K}_{2}\left(t\right)=\frac{{\beta }_{E}\left({T}_{w}\right)E\left(t\right)}{{K}_{{L}_{1}}}+{\beta }_{{L}_{1}}\left({T}_{w}\right).$$

(41)

$${K}_{3}\left(t\right)=\frac{{\beta }_{{L}_{1}}\left({T}_{w}\right){L}_{1}\left(t\right)}{{K}_{{L}_{2}}}+{\beta }_{{L}_{2}}\left({T}_{w}\right),$$

(42)

$${K}_{4}\left(t\right)=\frac{{\beta }_{{L}_{2}}\left({T}_{w}\right){L}_{2}\left(t\right)}{{K}_{{L}_{3}}}+{\beta }_{{L}_{3}}\left({T}_{w}\right),$$

(43)

$${K}_{5}\left(t\right)=\frac{{\beta }_{{L}_{3}}\left({T}_{w}\right){L}_{2}\left(t\right)}{{K}_{{L}_{4}}}+{\beta }_{{L}_{4}}\left({T}_{w}\right),$$

(44)

and

$${K}_{6}\left(t\right)=\frac{{\beta }_{{L}_{4}}\left({T}_{w}\right){L}_{4}\left(t\right)}{{K}_{P}}+{\beta }_{P} \left({T}_{w}\right),$$

(45)

respectively.

If we consider a specific starting point ${k}_{1}\omega$ for the Eqs. ($12)$ and $(15),$ we find that all the elements $\left(E\left(t\right),{L}_{1}\left(t\right),{L}_{2}\left(t\right),{L}_{3}\left(t\right),{L}_{4}\left(t\right), P\left(t\right),{A}_{m}\left(t\right)\right)$ will be positive for all $t\ge {k}_{1}\omega$. This contradicts the fact that $\partial {X}_{0}$ is positively invariant and hence ${M}_{\partial }\subset \left\{\left(\text{0,0},\text{0,0},\text{0,0},0\right)\right\}$. So from the two inclusions, ${M}_{\partial }=\left\{\left(\text{0,0},\text{0,0},\text{0,0},0\right)\right\}$.

The set ${M}_{\partial }=\left\{\left(\text{0,0},\text{0,0},\text{0,0},0\right)\right\}$ reveals that there exists a fixed point $\varepsilon =\left(\text{0,0},\text{0,0},\text{0,0},0\right)$ of the Poincare map P and acyclic in ${M}_{\partial }$, every solution in ${M}_{\partial }$ approaches$\varepsilon$. Moreover, Lemma 3.5 implies that $\varepsilon$ is an isolate invariant set in X and ${W}^{s}(\varepsilon )\cap {X}_{0}= \varnothing .$ Therefore, by Theorem 3.1.1 in [40], we finally obtain that P is uniformly persistent with respect to $\left({X}_{0},\partial {X}_{0}\right)$. So, the periodic semiflow $\varphi (t)$ is also uniformly persistent. This implies that there exists a positive constant $\beta$ such that any solution $\left(E\left(t\right), {L}_{1}\left(t\right), {L}_{2}\left(t\right), {L}_{3}\left(t\right), {L}_{4}\left(t\right), P\left(t\right), {A}_{m}\left(t\right)\right)$ with initial condition $\left(E\left(0\right), {L}_{1}\left(0\right), {L}_{2}\left(0\right), {L}_{3}\left(0\right), {L}_{4}\left(0\right), P\left(0\right), {A}_{m}\left(0\right)\right) \in {X}_{0}$ satisfies.

$\underset{t\to \infty }{\text{lim}}\text{inf}E\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{L}_{1}\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{L}_{2}\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{L}_{3}\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{L}_{4}\left(t\right)\ge \beta$, $\underset{t\to \infty }{\text{lim}}\text{inf}P\left(t\right)\ge \beta$,$\underset{t\to \infty }{\text{lim}}\text{inf}{A}_{m}\left(t\right)\ge \beta$.

This result is guaranteed by Theorem 1:3:6 in [39], which ensures the existence of at least one periodic solution for system (1). This solution, denoted by ${u}_{p}\left(t,{\varphi }^{*}\right)$, has an initial state ${\phi }^{*}$ that belongs to a specific set ${X}_{0}$. Because of this property of ${\phi }^{*}$, and according to Eqs. (35) to (44), all the components of ${u}_{p}\left(t,{\varphi }^{*}\right)$ will be positive for any value of t. In other words, this periodic solution has strictly positive values for each variable. This finding of a positive periodic solution completes the proof of the theorem.

Entomological data

This study used data on An. arabiensis mosquitoes' development from egg to adult to identify the effects of temperature. The experiment was conducted under two conditions: controlled insectary settings as a control, and a separate simulated environment where both temperature and humidity were allowed to fluctuate.

Mosquito egg hatching experiment

Four groups of 100 mosquito eggs were each placed on wet filter paper in separate petri dishes. After 24 h, the eggs were rinsed with distilled water and the filter paper was moved to hatching trays. The experiment tracked the number of eggs that hatched into first-stage larvae (${L}_{1}$) and those that remained unhatched, along with the hatching times.

Larval rearing and development

Following the egg hatching, four groups of first-stage larvae (${L}_{1}$), each containing approximately 100 individuals, placed into separate plastic containers (35 × 25 × 10 cm) filled to a depth of 4 cm with distilled water under controlled environmental conditions and in simulated conditions.

The larvae were provided a diet of finely ground tetrafin fish food or cat food and cerifam. The survival rates of the larvae were monitored as they progressed through each developmental stage (${L}_{1}$ to ${L}_{2}$, ${L}_{2}$ to ${L}_{3}$, ${L}_{3}$ to ${L}_{4}$, and ${L}_{4}$ to pupae). Temperature and humidity were also recorded throughout the study period. Additionally, the number of larvae that died at each stage (mortality rate) was documented.

Pupation and adult emergence

Each day, pupae were handpicked (collected) using a dropper. The time taken for larvae to develop into pupae and the total number collected were documented. The pupae were carefully transferred to individual small cups for monitoring. The experiment then tracked the number of days it took for adult mosquitoes to emerge from the pupae. Additionally, the number of larvae that failed to pupate (un-pupated individuals) and the number that died during the pupal stage (mortality rate) were recorded.

Adult mosquito housing and monitoring

After emerging from pupae, adult mosquitoes were kept in an individual metal cages. Each cage was assigned a unique identifier to track its corresponding replicate group and collection date. To ensure their survival, all adult mosquitoes received a sugar solution daily as a source of nutrition. The experiment tracked the daily mortality rates and lifespan of individual adult mosquitoes.

This comprehensive experiment provided valuable data on mosquito development, from egg hatching to adult emergence. The data allows to study factors that influence mosquito growth, survival rates, and overall development stages.

Estimation of parameters in the model

Feasibility of estimating entomological parameters used in model (1) such as egg development rate, larval instar mortality rate and development rate, pupae mortality rate and development rate, and adult mortality rate from a unique experiment involving whole aquatic phase parameters based on the real-world data.

To modify the entomological parameters based on temperature, we choose polynomial function of degree n,

$${P}_{i}\left(T\right)=\sum_{i=0}^{n}{a}_{i}{T}^{i}={a}_{n}{T}^{n}+{a}_{n-1}{T}^{n-1}+\cdots +{a}_{1}T+{a}_{0}.$$

(46)

In this equation, T represents the temperature in degrees Celsius, and the coefficients ${a}_{i}$ (where i = 0, 1,…,n) are determined using the linear least squares estimation method [40], which minimizes ${R}^{2}$,

$${R}^{2}=\sum_{j=1}^{N}\left[\frac{{P}_{n}\left({T}_{j}\right)-{\varphi }_{j}}{{\sigma }_{j}}\right],$$

(47)

where ${\varphi }_{j}$ is the mortality or transition rate at temperature ${T}_{j}$, with ${\sigma }_{j}$ being the corresponding error, and N is the total number of observed rates.

The fitted cubic function:

$${B}_{E}\left({T}_{a}\right)= 0.0188{\left({T}_{a}\right)}^{3}+ 0.0202{\left({T}_{a}\right)}^{2} + 0.0289{T}_{a}+ 0.1422$$

(48)

The fitted cubic function:

$${\beta }_{{L}_{1}}\left({T}_{w}\right) =0.0018{\left({T}_{w}\right) }^{3} + 0.0135{\left({T}_{w}\right) }^{2}+ 0.0411{T}_{w}+ 0.1126$$

(49)

The fitted cubic function:

$${\beta }_{{L}_{2}}\left({T}_{w}\right) = 0.0070{\left({T}_{w}\right) }^{3}+ 0.0219{\left({T}_{w}\right) }^{2}+ 0.0566{T}_{w}+ 0.1489$$

(50)

The fitted cubic function:

$${\beta }_{{L}_{3}}\left({T}_{w}\right) =0.0728{\left({T}_{w}\right) }^{3} + 0.0651{\left({T}_{w}\right) }^{2}-0.0224{T}_{w}+ 0.2167$$

(51)

The fitted cubic function:

$${\beta }_{{L}_{4}}\left({T}_{w}\right) =0.0239{\left({T}_{w}\right)}^{3}+ 0.0228{\left({T}_{w}\right)}^{2} + 0.0032{T}_{w}+0.1683$$

(52)

The fitted cubic function:

$${\beta }_{P}\left({T}_{w}\right) =0.0032{\left({T}_{w}\right)}^{2} -0.1737{T}_{w}+ 2.4345$$

(53)

The fitted curve:

$${\mu }_{E}\boldsymbol{ }\left({T}_{w}\right)= \boldsymbol{ }0.0036{\left({T}_{w}\right)}^{2}-\boldsymbol{ }0.2034{T}_{w}+3.1554$$

(54)

The fitted curve:

$${\mu }_{{L}_{1}}\boldsymbol{ }\left({T}_{w}\right)= \boldsymbol{ }0.0019{\left({T}_{w}\right)}^{2}\boldsymbol{ }-0.1064{T}_{w}+\boldsymbol{ }1.5162$$

(55)

Polynomial equation:

$${\mu }_{{L}_{2}}\boldsymbol{ }\left({T}_{w}\right)= 0.0022 {\left({T}_{w}\right)}^{2}\boldsymbol{ }-0.1281{T}_{w} +1.8257$$

(56)

Polynomial equation:

$${\mu }_{{L}_{3}}\boldsymbol{ }\left({T}_{w}\right)= 0.0021 {\left({T}_{w}\right)}^{2}\boldsymbol{ }-0.1174{T}_{w}+\boldsymbol{ }1.6715$$

(57)

Polynomial equation:

$${\mu }_{{L}_{4}}\boldsymbol{ }\left({T}_{w}\right)= 0.0016{\left({T}_{w}\right)}^{2}\boldsymbol{ }-0.1019{T}_{w}+\boldsymbol{ }1.6473$$

(58)

Fitted function:

$${\mu }_{p}\boldsymbol{ }\left({T}_{w}\right)= 0.0086 {\left({T}_{w}\right)}^{2}\boldsymbol{ }-0.4296 {T}_{w}+ 5.4571$$

(59)

Fitted function:

$${\mu }_{{A}_{m}}\boldsymbol{ }\left({T}_{w}\right)= 0.0016 {\left({T}_{w}\right)}^{2}\boldsymbol{ }-0.1019{T}_{w}+ 1.6473$$

(60)

Numerical simulations

In this section, fitted parameters in the model (Eqs. 48–60) were used for different temperature ranges (14.55 °C–34.35 °C), using MATLAB computer software and some parameter values was taken from [31, 33] to illustrate our analytic results.

Results

This section of the paper presents the results from both the experimental and numerical analyses of the impact of temperature on all developmental stages of An. arabiensis.

The results from the experimental data are presented in Figs. 3, 4, 5, 6, 7, 8, 9. Figures 3, 4, 5, 6, 7, 8 show the development rates for all Anopheles arabiensis developmental stages (eggs, larval stages 1–4, and pupa) across the temperature range of 14.55 °C to 34.35 °C. Figure 9 visualizes these development rates in a single graph. Figures 10, 11, 12, 13, 14, 15, 16, 17 illustrate the mortality rates for all Anopheles arabiensis developmental stages (eggs, larval stages 1–4, pupa, and adults) plotted against temperature across the range of 14.55 °C to 34.35 °C. Figure 17 visualizes these mortality rates in a single graph.

Fig. 3

A cubic polynomial regression model illustrating the relationship between temperature and the predicted number of egg development rate

Fig. 4

A cubic polynomial regression model illustrating the relationship between temperature and the predicted number of ${L}_{1}$ development rate

Fig. 5

A cubic polynomial regression model illustrating the relationship between temperature and the predicted number of ${L}_{2}$ development rate

Fig. 6

A cubic polynomial regression model illustrating the relationship between temperature and the predicted number of ${L}_{3}$ development rate

Fig. 7

Third-degree polynomial regression analysis of ${L}_{4}$ development rate as a function of temperature

Fig. 8

The predicted number of pupa development rate as a function of temperature with a polynomial function of degree two

Fig. 9

Development Rate of All An.Arbisansis Stages at 14.55 $^\circ{\rm C}$ to 34.35 $^\circ{\rm C}$ for the Experimental Data

Fig. 10

The predicted egg mortality rate as a function of temperature is fitted by a two-degree polynomial function

Fig. 11

The predicted ${L}_{1}$ mortality rate as a function of temperature is fitted by a second-degree polynomial function

Fig. 12

Relationship between Temperature and ${L}_{2}$ Mortality Rate

Fig. 13

The predicted number of ${L}_{3}$ mortality rate as a function of temperature with a polynomial function of degree two

Fig. 14

The predicted number of ${L}_{4}$ mortality rate as a function of temperature with a polynomial function of degree two

Fig. 15

The predicted number of pupa mortality rate as a function of temperature with a polynomial function of degree two

Fig. 16

The fitting of female adult mosquitoe mortality rate as a function of temperature

Fig. 17

Mortality Rate of All Anopheles arabiensis Stages Against Temperature (14.55 °C to 34.35 °C)

As one can see from Figs. 9 and 17, the development of An. arabiensis was experimentally observed to cease below 14.55 °C or above 34.35 °C. The estimated lower and upper developmental thresholds for An. arabiensis were approximately 14.55 °C and 34.35 °C, respectively.

Furthermore, the temperature at which all developmental stages of An. arabiensis reached their peak development rate was at 32.57 °C, as shown in Fig. 9.

The total mortality rate of egg to adult of An. arbiensis as in Fig. 17 shows that pupa mortality rates increase from 14.55 °C to 15.62 °C, peaking again at 15.62 °C and reaching a minimum at 28.37 °C and 29.65 °C. In contrast, adult mortality rate is nonexistent between 14.55 °C and 17.57 °C, rises from 27.55 °C to 31.55 °C, declines from 31.55 °C to 33.35 °C, and reaches its maximum at 34.35 °C. Notably, the mortality rates of all immature stages (larval stages 1–4 and pupa) and adult stages are comparable within the temperature range of 26.62 °C to 28.37 °C.

Furthermore, all developmental stages of Anopheles arabiensis peak in mortality rate between 33.35 °C and 34.35 °C.

Figure 18A–G illustrate the population dynamics of Anopheles arabiensis across different life stages (egg, larval stages 1–4, pupa, and adult) at varying temperatures (14.55 °C–34.35 °C). The x-axis of each graph represents time (in days), while the y-axis represents population size. These figures provide a visual representation of how temperature influences the growth and survival of Anopheles arabiensis populations.

Fig. 18

A–G represent graphs showing the population changes of Anopheles arabiensis, specifically for larval stages 1, 2, 3, and 4, eggs, pupae, and female adults, over time at different temperatures

Numerical analysis revealed that immature Anopheles arabiensis life stages peaked in abundance at temperatures above 29.65 °C. Conversely, as clearly shown in Fig. 18A–G, temperatures below 14.55 °C negatively impacted all developmental stages of An. arabiensis.

Discussion

In this paper, we developed and analyzed a mathematical model incorporating logistic growth and temperature to determine the impact of temperature on each developmental stage of An. arabiensis: eggs, the four larval stages (four instars), pupae, and adults.

First, we established the epidemiological and mathematical soundness of the proposed model. Then, we demonstrated that if the non-autonomous basic reproduction number is less than unity, the disease-free equilibrium is locallly asymptotically stable. Conversely, if it is greater than unity, there exists at least one positive periodic solution, as established by applying the comparison theorem and the theory of uniform persistence.

The dynamics of both immature and adult Anopheles arabiensis populations were assessed across a range of temperatures: 14.55 °C, 15.62 °C, 16.62 °C, 17.57 °C, 18.53 °C, 19.42 °C, 20.6 °C, 21.52 °C, 22.5 °C, 23.65 °C, 24.72 °C, 25.3 °C, 26.62 °C, 27.55 °C, 28.37 °C, 29.65 °C, 31.55 °C, 32.57 °C, 33.35 °C, and 34.35 °C.

Previous studies have examined the relationship between temperature and the development rate of other Anopheles species ([41] and references therein), as well as other mosquito vectors [42]. Bayoh and Lindsay [15] indicated that increasing environmental temperature on the developmental stages of Anopheles mosquitoes results in faster development rates, peaking at 28 °C in their study. Experimental data, as in [41], revealed that the fastest development rate occurs at 31 °C, and as in [42], the maximum development rate of An. arabiensis occurs at higher temperatures of 32 °C, which is in excellent agreement with our findings of 32.57 °C, as shown in Fig. 10.

Figures 9 and 17 indicate that all developmental stages of Anopheles arabiensis ceased development below 14.55 °C and above 34.35 °C. This result aligns well with findings in [15, 41,42,43,44].

The rate of mortality for Anopheles arabiensis larvae in stages 1–4 (or four instars) was highest at temperatures between 15.62 °C and 16.62 °C and above 31.55 °C, as shown in Fig. 17. This aligns well with the finding that Anopheles gambiae larvae mortality rates peak at temperatures below 18 °C and above 32 °C, as reported in [16]. Similar trends have been observed in Culex tarsalis, where mortality increases below 12 °C and at 32 °C or above [45,46].

To numerically validate our analytical results, we obtained parameter values from [31, 33] and computed the remaining parameters through model fitting based on real-world data from the Tropical and Infectious Diseases Research Center (TIDRC) in Sekoru, Jimma University.

As clearly indicated in all Fig. 18A–G, temperatures below 14.55 °C negatively impact An. arabiensis, aligning with the findings reported in [44] and our entomological data collected from JU TIDRC. Additionally, immature mosquito populations are more sensitive to temperatures at 29.65 °C, in good agreement with [44], whereas adult female mosquitoes are more sensitive above 34.35 °C.

Conclusions

In this study, we used both entomological data and a mathematical model to analyze the impact of temperature variability on each developmental stage of An. arabiensis. The basic reproduction number is computed by employing the method used in Wang and Zhao [34], which extended the framework of Van den Driessche and James Watmough [35]. If the basic reproduction number is less than unity, the disease-free equilibrium is locally asymptotically stable. Conversely, if it is greater than unity, there exists at least one positive periodic solution, as established by applying the comparison theorem and the theory of uniform persistence.This study contributes to the existing knowledge on the relationship between temperature and each developmental stage of An. arabiensis, utilizing both entomological actual data conducted at TDRC and a mathematical model. The findings provide valuable information for modeling vector population dynamics in the context of climate variability and designing tailored malaria vector control strategies. Future studies could further refine our understanding of An. arabiensis population dynamics by incorporating additional factors such as humidity.

Limitations of the model

It is essential to acknowledge the inherent limitations of mathematical models. They often serve as approximations of complex real-world phenomena, inevitably introducing some degree of inaccuracy. Furthermore, the parameters employed in these models are typically derived from empirical observations and experiments, often relying on numerical methods and software. This process can introduce uncertainty into the model's inputs. Consequently, readers should exercise caution when interpreting the results, considering the potential discrepancies between model predictions and actual outcomes.