220
Views
0
CrossRef citations to date
0
Altmetric
Research Article

Mathematical model of mosquito population dynamics with constants and periodic releases of Wolbachia-infected males

, &
Article: 2372581 | Received 09 Oct 2023, Accepted 20 Jun 2024, Published online: 08 Jul 2024

Abstract

In this article, we propose three mathematical models of mosquito dynamics. The first model describes the dynamics of wild mosquitoes in combination with mechanical control and larvicide. For this model, we derive a threshold parameter N, called the basic number of offspring, and show that the trivial equilibrium is globally asymptotically stable when N 1, while if N>1, the non-trivial equilibrium is globally asymptotically stable. In the second model, we evaluated a constant release of Wolbachia-infected males in the wild mosquito population. Our analysis identified a critical threshold for releasing Wolbachia-infected male mosquitoes Micrit to reduce or eliminate the wild mosquito population. In the third model, we assessed wild mosquito dynamics by periodically releasing Wolbachia-infected male mosquitoes. The analysis of this model produced two thresholds, T the threshold on the release period and Λ the threshold on the release quantity. The existence of trivial periodic solutions is obtained, and the corresponding local and global stability conditions are proved by Floquet's theory and Lyapunov's stability theorem, respectively. We also prove the existence conditions for non-trivial periodic solutions and their local stability. Finally, we perform some numerical simulations to illustrate our theoretically pertinent results.

2020 Mathematics Subject Classifications:

1. Introduction

Mosquito-borne diseases are caused by bacteria, viruses or parasites transmitted by mosquitoes. Every year, nearly 700 million people suffer from mosquito-borne diseases, leading to more than a million deaths and a major economic burden in endemic countries [Citation1]. In recent years, with global warming, ease of transportation and environmental degradation, the geographic environment suitable for mosquitoes has expanded, and various mosquito-borne diseases have continuously spread on a large scale. Anopheles, Aedes and Culex are three major mosquito genera, which comprise approximately 3,500 species, are found beyond the tropical and subtropical regions of the world. Different kinds of mosquitoes can transmit different pathogens. For example, Anopheles primarily transmit malaria and filariasis, Aedes can transmit Zika, dengue, yellow fever, and chikungunya, while Culex can transmit West Nile, Japanese encephalitis, and filariasis [Citation2]. At present, there are no effective vaccines or therapeutic drugs to combat some mosquito-borne diseases, so controlling mosquito vectors has proven to be a suitable method to prevent the spread of transmitted diseases by mosquitoes [Citation3–5]. However, excessive use of traditional methods to suppress the mosquito population is generally unhealthy or unsustainable for humans and the environment. For example, heavy use of chemical insecticides can lead to insecticide resistance and toxicological accumulation [Citation3]. Thus, the search for new technologies aimed at breaking the cycle of transmission of diseases transmitted by mosquitoes has been applied in recent years. Notably, the Sterile Insect Technique (SIT) which is a method of biological insect control and requires the release of a large number of sterile insects into the wild. This technique was developed in the 1940s and 1950s and adopted for Anopheles, Culex, and Aedes mosquito control in the early 1970s and is now being utilized and accepted by some mosquito control programs. The released mosquitoes are preferably sterile males that compete with wild males to mate with females of the natural populations. After mating with a sterile male, females produce no offspring, thus reducing the next generation's population. This is an environmentally friendly control method involving mass-rearing and sterilization by radiation and other methods. In addition, we have the Cytoplasmic Incompatibility (CI) caused by endosymbiotic bacteria, Wolbachia [Citation6]. Wolbachia is present in up to 60% of insect species including some mosquitoes. This bacterium induce cytoplasmic incompatibility (CI) in mosquitoes, which causes early embryonic mortality resulting from mating between Wolbachia-infected males and females which are either uninfected or infected with a different Wolbachia strain [Citation7–9].

Mathematical models are useful tools for finding the best strategies in the implementation of different types of vector control interventions and for assessing their impact on the disease propagation [Citation10–14]. Some authors had studies only one type of vector control interventions, such as insecticide spraying [Citation15, Citation16], Sterile Insect Technique or SIT, [Citation17–22], Wolbachia infestation [Citation23–27] and introduction of natural predators of mosquito larvae [Citation28, Citation29]. Others had focused on combinations of various types of vector control, namely, the combined effects of insecticide spraying with personal protection measures [Citation30] or educational campaigns [Citation31], interventions that complement the insecticide spraying with the use of larvicides [Citation19, Citation21, Citation32, Citation33].

Our model is inspired by the work of [Citation18, Citation33–35], by taking into account the total adult mosquito population in the mating process. Then, we introduce mechanical and larvicide control in the dynamics of the aquatic stage mosquitoes, combined with a constant and impulsive release of Wolbachia-infected males. Our main objective is to evaluate the impact of Wolbachia-infected male releases combined with mechanical and larvicidal control on the growth dynamics of wild mosquitoes. To achieve our objectives, we considered a simple model of the growth dynamics of mosquitoes combined with mechanical and larvicidal control and obtained a threshold N. We then disrupted the dynamics of wild mosquitoes by releasing Wolbachia-infected males in a constant and impulsive manner.

The rest of the paper is organized as follows. In Section 2, we propose and analyse a mathematical model for the population dynamics of wild mosquitoes combined with mechanical control and larvicide treatment. This model was extended by considering a constant release of Wolbachia-infected male in Section 3. In addition, we present a qualitative analysis of this model. Specifically, we calculate the equilibria and study their stabilities. In Section 4, we establish an impulse model by periodically releasing Wolbachia-infected male mosquitoes. We then perform a dynamic analysis of the model to obtain the existence and stability of the periodic solution. In Section 5, we perform some numerical simulations to illustrate the results obtained. Finally, the conclusion of this paper is given in the last section.

2. Mathematical model of wild mosquito dynamics with mechanical and larvicide control

2.1. Mathematical formulation of the wild mosquito model

The life cycle of a mosquito consists of two main stages: aquatic (egg, larva, pupa) and adult (males and females). A few days after being fertilized, a female mosquito may deposit a few dozen eggs possibly divided between several breeding sites. Once laid, the eggs of some species can resist up to several months before hatching. After stimulation, the eggs hatch to give rise to larvae that will develop in the water and reach the stage of pupae. The larval stage lasts from a few days to a few weeks.

For the mathematical description, we consider two stages: the aquatic stage and adult stage. We denote by A the aquatic stage (egg, larva, pupa) and we divided the adult stage in three sub-compartments: M wild males, F young females not yet laying eggs, Fs fertilized and egg-laying females. We consider that a female to be in the F compartment from its emergence from aquatic stage until her gonotrophic cycle has began, that is the time of mating and taking the first blood meal, which takes typically 3−4 days. The vector control terms added to the model are the application of larvicide (chemical control) and the removal of potential breeding sites (mechanical control). The description of the control terms is as follows:

  • θ is defined as the parameter associated with the effectiveness of mechanical control, which consists of destroying the breeding sites.

  • ω is the mortality rate caused by the use of the larvicide that will affect the aquatic stage of the mosquito population.

The life cycle is described using the diagram in Figure . The variables and parameters of the system are detailed in Tables  and , respectively, where we indicate each variable and each parameter of the system giving its physical meaning and its quasi-dimensional unity [Citation36–39].

Figure 1. Wild mosquito flow chart.

Figure 1. Wild mosquito flow chart.

Table 1. Types of state variables and their quasi-dimensional unit which is either a for aquatic life stage and f for adult mosquitoes or terrestrial forms of the mosquito.

Table 2. Parameters of model (Equation1), their description and corresponding quasi-dimension.

Thus, the corresponding mathematical model is given by system (Equation1): (1) {A˙(t)=ρϕ(1A(t)θK)Fs(t)(νA+ω+μA)A(t),M˙(t)=(1r)λνAA(t)μMM(t),F˙(t)=νAA(t)(β+μF)F(t),F˙s(t)=βM(t)M(t)+F(t)F(t)μFsFs(t),(1) with non-negative initial condition, (2) A(0)=A0,M(0)=M0,F(0)=F0,Fs(0)=Fs0.(2) For, the variation of the population of male mosquitoes M, the inflow is (1r)λνAA, which represents the number of mosquitoes that develop to male adults per unit of time and the outflow is represented by μMM, which expresses the male mosquitoes that die per unit of time. For the variation of young females not yet laying eggs F, the inflow is νAA, which indicates the number of mosquitoes that develop to young female adults. The parameter 1/(β+μF) can represents the average time given by the length of a female's first gonotrophic cycle of a female, i.e. the interval between immediately after mating and the first blood meal [Citation18]. As mating is a complex process, as indicated in [Citation18, Citation40], the male mosquito can mate practically throughout his life. A female mosquito needs successful mating to reproduce throughout her life [Citation18, Citation40]. Note that mosquitoes locate themselves in space and time to ensure they are available to mate. Thus, we consider that the total number of mosquitoes in the aerial phase is expressed by the term M + F, and the fraction M/(M+F) represents the probability that a young females meets a male in the adult phase. Then, MF/(M+F) is the number of encounters between males and young females. So, βMF/(M+F) is the number of matings between young females and males per unit of time. μFF (resp. μsFs) is the number of young females (resp. fertilized and egg-laying females) mosquitoes that die per unit of time. ρϕFs represents the number of eggs per unit of time laid by fertilized and egg-laying females without considering the number of aquatic mosquitoes that the medium can support. (1A/θK) indicates the probability that there is a place available for more aquatic mosquitoes and (ω+μA)A which represents the aquatic mosquitoes that die per unit of time.

2.2. Mathematical analysis of the wild mosquito model

In this section, we present some theoretical results of system (Equation1). First, we present the basic properties of the solutions of system (Equation1). These are useful for showing that system (Equation1) is well posed mathematically and biologically. Indeed, note that system (Equation1) can be written in matrix form as follows: (3) X˙=P(X)X,(3) where X=(A,M,F,Fs)T and P(X)=((ρϕFsθK+νA+ω+μA)00ϕ(1r)λνAμM00νA0(β+μF)000βMM+FμFs)The right-hand side of system (Equation1) is continuous and infinitely differentiable on R+4. Therefore, it is locally Lipschitz on R+4, which guarantees the existence and uniqueness of the solution of system (Equation1) associated with the initial conditions X(0)=X0R+4.

We have the following result.

Proposition 2.1

The set Ω={XR+4:AθK, M(1r)λνAθKμM, FνAθKβ+μF,FsνAβθKμFs(β+μF)}is positively invariant for system (Equation1).

Proof.

We prove it in two steps:

Step 1:

We show that the solution variables (A(t),M(t),F(t),Fs(t)) of system (Equation1) corresponding to initial conditions such that A(0)>0,M(0)>0,F(0)>0 and Fs(0)>0 are non-negative. Let t1=sup{t>0/ u[0,t]:A(u)>0,M(u)>0,F(u)>0,Fs(u)>0}.The initial conditions above in addition to the continuity of all functions A, M, F and Fs guarantee the existence of t1.

If t1=+, then all solutions of system (Equation1) are positive. Suppose t1<+, then there exists at least one solution A(t),M(t),F(t) or Fs(t) that is equal to zero at t1. Suppose, for example that A(t1)=0, then the first equation of system (Equation1) yields A˙=ρϕFs(ρϕFsθK+(νA+ω+μA))Awhich may be rewritten as follows ddt[A(t)exp(0tρϕFs(v)θK+ϑdv)]=ρϕFs(t)exp(0tρϕFs(v)θK+ϑdv),with ϑ=νA+ω+μA.Integrating the above inequality from 0 to t1 gives A(t1)=[A(0)+0t1ρϕFs(t)exp(0tρϕFs(v)θK+ϑdv)dt]×exp(0t1ρϕFs(v)θK+ϑdv)>0,which is in contradiction with A(t1)=0. Similarly, one can show that M(t1)>0,F(t1)>0 and Fs(t1)>0. Thus, A(t)>0,M(t)>0,F(t)>0 and Fs(t)>0 for all t>0.

Step 2:

From the first equation of system (Equation1), we will show that ∀t>0,A(t)θK.

Indeed, suppose that there exists t0>0 such that A(t0)>θK. By continuity, there exists t1<t0 such that A(t1)=θK and ∀t]t1,t0], we have A(t)>θK. Thus, if A(t1)=θK, from the first equation of system (Equation1), we have (4) A˙(t1)=(νA+ω+μA)θK<0.(4) Then, there exists t2]t1,t0[ such that A(t2)<θK which is absurd because A(t)>θK, ∀t]t1,t0]. Thus, for all t>0, we have A(t)θK.

Using the fact that A(t)θK,∀t>0, the second equation in system (Equation1) gives (5) M˙(t)(1r)λνAθKμMM(t).(5) Applying the Gronwall's theorem to Equation (Equation5), one has that M(t)(1r)λνAθKμM+(M(0)(1r)λνAθKμM)eμMt,where M(0) represents the initial value of M(t). It then follows that lim suptM(t)(1r)λνAθKμM,therefore M(t)(1r)λνAθKμM for all t>0.Using the fact that A(t)θK,∀t>0, the third equation in system (Equation1) gives (6) F˙(t)νAθK(β+μF)F(t).(6) Applying the Gronwall's theorem to Equation (Equation6), one has that F(t)νAθKβ+μF+(F(0)νAθKβ+μF)e(β+μF)t,where F(0) represents the initial value of F(t). It then follows that lim suptF(t)νAθKβ+μF,therefore F(t)νAθKβ+μF for all t>0.Now, using the fact that F(t)νAθKβ+μF for all t>0, the fourth equation in system (Equation1) gives (7) F˙s(t)βνAθKβ+μFμFsFs(t).(7) Applying the Gronwall's theorem to Equation (Equation7), one has that Fs(t)βνAθKμFs(β+μF)+(Fs(0)βνAθKμFs(β+μF))eμFst,

where Fs(0) represents the initial value of Fs(t). It then follows that lim suptFs(t)βνAθKμFs(β+μF),therefore Fs(t)βνAθKμFs(β+μF) for all t>0.

Thus, Ω is a positively invariant set under the flow of system (Equation1). Hence, it is sufficient to consider the dynamics of model (Equation1) in Ω. This concludes the proof.

Now, we will derive some quantitative analysis of system (Equation1). Let us consider the following threshold N=(1r)νAβρϕμFs(νA+ω+μA)[(1r)(β+μF)+rμM].The introduction of the threshold parameter N is intricately connected to the equilibria analysis of system (Equation1). This parameter quantifies the average expected number of new adult female offspring produced by a single female vector throughout its lifetime, often denoted as the basic offspring number [Citation18, Citation22, Citation34, Citation41]. It quantifies the reproductive potential of individuals in a population and is crucial for evaluating the sustainability of the mosquito population. The threshold N holds particular importance in Population Dynamics and its significance extends to population biology, as it effectively encapsulates the reproductive potential of individuals. This, in turn, offers valuable insights into the mosquito population's capacity to sustain stability over time. Furthermore, stability analysis emphasizes the importance of N in determining the behaviour of the system around equilibrium points. If N is less than or equal to 1, it indicates a critical threshold where the population struggles to sustain itself in the environment. Conversely, if N is greater than 1, the mosquito population persists in the environment. The comparison of N to 1 serves as a biological metric to assess the potential of the mosquito population to persist or decline. This comparison is rooted in the mathematical relationships governing the equilibrium conditions, highlighting the direct connection between mathematical analysis and biological interpretation.

Lemma 2.2

System (Equation1) always has the trivial equilibrium point E0=(0,0,0,0).

(i)

If N1, then system (Equation1) has no other equilibria.

(ii)

If N>1 then system (Equation1) admits a unique non-trivial equilibrium point E1=(A,M,F,Fs), defined as follows: (8) {A=μMM(1r)λνA,M=(1r)(N1)λνAθKμMN,F=rμMM(1r)(β+μF),Fs=rμMβMμFs[(1r)(β+μF)+rμM].(8)

Proof.

The equilibrium point of system (Equation1) is found by solving the following system: {A˙(t)=0,M˙(t)=0,F˙(t)=0,F˙s(t)=0.{ρϕ(1AθK)Fs(νA+ω+μA)A=0,(1r)λνAAμMM=0,νAA(β+μF)F=0,βMM+FFμFFs=0.It is clear that (0,0,0,0) satisfies the previous equations. The second equilibrium point is obtained quickly by solving the system with, for example, the substitution method [Citation42, Citation43].

Consider a n-dimensional autonomous differential system: (9) x˙=f(x),x(0)=x0,(9) where f is a given vector function, i.e. f=(f)i, with fi:R+nR. Then, we recall the following general definitions [Citation18, Citation44]:

Definition 2.3

System (Equation9) is called cooperative if for every i,j{1,2,,n} such that ij the function fi(x1,x2,,xn) is monotone increasing with respect to xj.

For cooperative systems, like system (Equation9), the global asymptotic stability of an equilibrium can be studied using the following theorem:

Theorem 2.4

Let a,bΩ such that a<b, [a,b]Ω and f(b)0f(a); where [a,b]={xRn,aixibi,i=1,2,,n}.Then, system (Equation9) defines a (positive) dynamical system on [a,b]. Moreover, if [a,b] contains a unique equilibrium p then p is globally asymptotically stable on [a,b].

Now, we use the fact that system (Equation1) is a cooperative system [Citation18, Citation44] to prove the global asymptotic stability of the trivial and non-trivial equilibrium points (E0=(0,0,0,0) and E1=(A,M,F,Fs)). The dynamic of system (Equation1) is summarized in the following theorem:

Theorem 2.5

Assume that the initial condition (A(0),M(0),F(0),Fs(0))Ω.

(i)

When N1, the trivial equilibrium E0 is globally asymptotically stable, which means that the mosquito population will decrease until extinction, whatever the initial population.

(ii)

When N>1, the trivial equilibrium E0 is unstable and the positive equilibrium E1 is globally asymptotically stable, which means that the mosquito population will persist.

Proof.

It suffices to verify the assumptions of Theorem 2.4.

  1. When N1, system (Equation1) has only the trivial equilibrium E0. By taking a=E0 and b=(θK,2(1r)λνAθKμM,2νAθKβ+μF,2νAβθKμFs(β+μF)),one has that f(a)=0andf(b)=((νA+ω+μA)θK(1r)λνAθKrνAθK2r2λμMνAβθK(β+μF)[(1r)(β+μF)+rμM])0.Thus, according to Theorem 2.4, the trivial equilibrium E0 is globally asymptotically stable on [a,b] hence on Ω, when N1.

  2. When N>1, there exists ϵ>0 such that N>1+ϵ. Let Aϵ sufficiently small such that Aϵϵ,Mϵ=(1r)λνAμMAϵ,Fϵ=νAβ+μFAϵ,F=(νA+ω+μA)NρϕAϵ.Let bϵ=(Aϵ,Mϵ,Fϵ,F). Then, from the right-hand side of system (Equation1) and the fact that N>1+ϵ and θK1, we deduce: f(bϵ)=(ρϕ(1AϵθK)F(νA+ω+μA)Aϵ(1r)λνAAϵμMMϵνAAϵ(β+μF)FϵβMϵMϵ+FϵFϵμFsF)(ϵ(νA+ω+μA)(11+ϵθK)000)0.Therefore, according to Theorem 2.4, equilibrium E1 is globally asymptotically stable on [bϵ,b]. Moreover, since bϵ can be selected to be smaller than any x>0, we have that E1 is asymptotically stable on Ω with basin of attraction Ω~=Ω∖{(0,0,0,0)}. This also implies that E0 is unstable.

3. Mathematical model of mosquito population dynamics with constant release of Wolbachia-infected males

3.1. Formulation of the mathematical model

In this section, we formulate an extension of model (Equation1). Indeed, we introduce the compartment Mi which represents the population of Wolbachia-infected male mosquitoes. Afterwards, we introduce a release quantity of Wolbachia-infected male mosquitoes Λ and the natural mortality of infected males, which is noted μi. Then, we divide the population of adult females into three compartments, young females not yet laying eggs F, females fertilized by males M noted Fs and the females fertilized by Wolbachia-infected male mosquitoes Mi noted Fst. A female from compartment F must mate with a wild male M to be able to pass into the compartment of fertilized females Fs and be able to deposit her eggs. Similarly, a female from compartment F must mate with a Wolbachia-infected male mosquitoes Mi to be able to move to the compartment Fst and be able to deposit her eggs but the eggs do not hatch. We consider that μst represents the natural mortality rate of Wolbachia-infected females. Thus, the compartmental diagram is shown in Figure .

Figure 2. Flow chart of the interaction of Wolbachia-infected male mosquitoes and wild mosquitoes.

Figure 2. Flow chart of the interaction of Wolbachia-infected male mosquitoes and wild mosquitoes.

The flow of F to the compartments Fs and Fst depends mainly on the number of encounters of young females not yet laying eggs with wild males and Wolbachia-infected male mosquitoes respectively, and the corresponding mating rates between them. We consider that after a release of Wolbachia-infected male mosquitoes, the total population of mosquitoes in the aerial phase can be expressed by M+Mi+F. The fraction M/(M+Mi+F) represents the probability that a young females not yet laying eggs meets a wild male and MF/(M+Mi+F) is the number of encounters between wild males and young females not yet laying eggs. Thus, βMF/(M+Mi+F) is the number of matings between young females and wild males. Similarly, the fraction Mi/(M+Mi+F) represents the probability that a young females not yet laying eggs meets a Wolbachia-infected male and MiF/(M+Mi+F) is the number of encounters between Wolbachia-infected male and young females not yet laying eggs. Thus, βMiF/(M+Mi+F) is the number of matings between young females and Wolbachia-infected male mosquitoes.

The mathematical equations that integrate the assumptions given above are: (10) {A˙(t)=ρϕ(1A(t)θK)Fs(t)(νA+ω+μA)A(t),M˙(t)=(1r)λνAA(t)μMM(t),F˙=νAA(t)(β+μF)F(t),F˙s(t)=βM(t)M(t)+Mi(t)+F(t)F(t)μFsFs(t)F˙st(t)=βMi(t)M(t)+Mi(t)+F(t)F(t)μstFst(t),M˙i(t)=ΛμiMi(t).(10) Since variable Fst is not present in the first four equations of system (Equation10), thus, according to [Citation22] the dynamics of system (Equation10) is equivalent to that of the following system: (11) {A˙(t)=ρϕ(1A(t)θK)Fs(t)(νA+ω+μA)A(t),M˙(t)=(1r)λνAA(t)μMM(t),F˙(t)=νAA(t)(β+μF)F(t),F˙s(t)=βM(t)M(t)+Mi(t)+F(t)F(t)μFsFs(t),M˙i(t)=ΛμiMi(t).(11)

3.2. Mathematical analysis of the model

In this section, the existence, uniqueness, positivity and boundedness of the solutions are demonstrated in the same way as the results obtained in the mathematical study of system (Equation1). Note that the dynamics of Wolbachia-infected males mosquitoes, i.e. the last equation of system (Equation11) is independent of other variables, but it can impact the dynamics of others. We can verify that its equilibrium Mi=Λ/μi is locally stable, and that for any initial value Mi(0), (12) Mi(t)=Mi+eμit(Mi(0)Mi)Mi,as t.(12) Hence equilibrium Mi is globally asymptotically stable. Moreover, we can obtain the limiting system of system (Equation11) as follows (13) {A˙(t)=ρϕ(1A(t)θK)Fs(t)(νA+ω+μA)A(t),M˙(t)=(1r)λνAA(t)μMM(t),F˙(t)=νAA(t)(β+μF)F(t),F˙s(t)=βM(t)M(t)+Mi+F(t)F(t)μFsFs(t).(13) Note that, the dynamics of system (Equation11) is qualitatively equivalent to their limiting systems (Equation13) [Citation19, Citation22, Citation45]. Thus, the principal results for system (Equation11) can be obtained by the same method as for system (Equation13).

3.2.1. The trivial equilibrium and its stability analysis

The trivial equilibrium, of system (Equation11) is obtained by setting the right-hand sides of the equations in system (Equation11) to zero. This is given by: E2=(A,M,F,Fs,Mi)=(0,0,0,0,Mi),corresponding to the state where wild mosquitoes are absent and there is only a constant population of Wolbachia-infected males.

Now, we analyse the stability conditions of the trivial equilibrium point. To do this, we calculate the Jacobian matrix of system (Equation11), given by J=((ρϕFsθK+νA)00ρϕ(1AθK)0(1r)λνAμM000νA0(β+μF)000βF(F+Mi)(M+Mi+F)2βM(M+Mi)(M+Mi+F)2μFsβMF(M+Mi+F)20000μi)

Theorem 3.1

The trivial equilibrium of system (Equation11), given by E2=(0,0,0,0,Mi), is locally asymptotically stable.

Proof.

The Jacobian matrix of system (Equation11) at the trivial equilibrium point E2 has the following form JE2=((νA+ω+μA)00ρϕ0(1r)λνAμM000νA0(β+μF)00000μFs00000μi)The eigenvalues of JE2 at the equilibrium point E2=(0,0,0,0,Mi) are (νA+ω+μA), μM, (β+μF), μFs and μi which are all negative, therefore E2 is always locally asymptotically stable.

Theorem 3.2

The trivial equilibrium of system (Equation11), given by E2, is globally asymptotically stable whenever N1 and unstable if N>1.

Proof.

Consider the Lyapunov function for system (Equation11): (14) V=μFs[(1r)(β+μF)+rμM](1r)rρϕβA+M+F+(1r)(β+μF)+rμM(1r)rβFs(14) with Lyapunov derivative given by V˙=μFs[(1r)(β+μF)+rμM](1r)rρϕβA˙+M˙+F˙+(1r)(β+μF)+rμM(1r)rβFs˙=μFs[(1r)(β+μF)+rμM](1r)rρϕβ[ρϕ(1AθK)Fs(νA+ω+μA)A]+[(1r)λνAAμMM]+[νAA(β+μF)F]+(1r)(β+μF)+rμM(1r)rβ[βMM+Mi+FFμFsFs]=[λνAμFs[(1r)(β+μF)+rμM](νA+ω+μA)(1r)rρϕβ]AμMM+(1r)(β+μF)+rμM(1r)rMM+Mi+FF(β+μF)FμFs[(1r)(β+μF)+rμM](1r)rρϕβAFsθK.Then, using the fact that there exists a ϵ]0,1] such that MM+Mi+Fϵ, and taking ϵ=(1r)r(β+μF)(1r)(β+μF)+rμM, we have the following inequality: V˙[λνAμFs[(1r)(β+μF)+rμM](νA+ω+μA)(1r)rρϕβ]AμFs[(1r)(β+μF)+rμM](1r)rρϕβAFsθKμMM(1r)rρϕβμFs[(1r)(β+μF)+rμM](νA+ω+μA)(N1)AμFs[(1r)(β+μF)+rμM](1r)rρϕβAFsθKμMM(1r)rρϕβμFs[(1r)(β+μF)+rμM](νA+ω+μA)(N1)A.Since all the model parameters are non-negative, it follows that V˙0 for N1. The largest compact invariant set in {(A,M,F,Fs,Mi)R+5:V˙=0} is the singleton {E2}. It follows from the LaSalle Invariance Principle [Citation19, Citation46, Citation47] that every solution of system (Equation11) with initial conditions in R+5 converge to the trivial equilibrium E2 as t. Thus (A,M,F,Fs,Mi)(0,0,0,0,Mi) when t for N1. So, From the LaSalle principle we deduce the attractiveness of E2. Since E2 is locally asymptotically stable when N1, we deduce that it is not only attractive, but it is also globally asymptotically stable. This ends the proof.

3.2.2. The non-trivial equilibrium and its stability analysis

3.2.2.1 Existence of non-trivial equilibrium

In this subsection, we consider that the basic offspring number N is greater than unity. By setting the right-hand sides of the equations of system (Equation11) to zero and expressing the variables in terms of M, we obtain the non-trivial equilibrium. Thus, the non-trivial equilibrium point E3=(A,M,F,Fs,Mi) satisfy the following relations: {A=μMM(1r)λνA,F=rμMM(1r)(β+μF),Fs=rμMβ(M)2μFs[((1r)(β+μF)+rμM)M+(1r)(β+μF)Mi],where M is a solution of the second degree equation (15) g(M)=a1(M)2+a2M+a3=0,(15) with coefficients a1=rμMβρϕμFsθK(νA+ω+μA)((1r)(β+μF)+rμM),a2=1(1r)νAβρϕμFs(νA+ω+μA)[(1r)(β+μF)+rμM]=1N,a3=(1r)(β+μF)Mi((1r)(β+μF)+rμM).It should be noted that the coefficients a1 and a3 of the quadratic Equation (Equation15) are always positive, however the coefficient a2 is negative if the basic offspring number N is greater than unity. It follows that g(M) has two or no positive roots if and only if N>1 and a224a1a3>0 or N>1 and a224a1a3<0, respectively. Then, the conditions for biological existence of the non-trivial equilibrium points become (16) N>1andMi<rλ2νA2βρϕθK(1r)4μMμFs(νA+ω+μA)(β+μF)(11N)2.(16) Let Micrit=rλ2νA2βρϕθK(1r)(11N)24μMμFs(νA+ω+μA)(β+μF).Then, we have (N1)2Micrit=4μMμFs(νA+ω+μA)(β+μF)N2rλ2νA2βρϕθK(1r).So condition (Equation16) becomes: N>1 and Mi<Micrit . Thus, the solutions are given by M=λνAθK(1r)(N1)2μMN[11MiMicrit],andM+=λνAθK(1r)(N1)2μMN[1+1MiMicrit].Moreover, under conditions (Equation16), we have N1>0 and MiMicrit<1,thus M>0 and M+>0.Therefore, under conditions (Equation16), system (Equation11) has two positive equilibria E3 and E3+, corresponding to M and M+.

Thus, the following result have been established:

Theorem 3.3

Assume that N>1.

(i)

If Mi>Micrit , then system (Equation11) has no positive equilibrium point.

(ii)

If Mi<Micrit , then system (Equation11) admits two positive equilibria points E3±=(A±,M±,F±,Fs±,Mi)where {A±=μMM±(1r)λνA,M±=νAθK(1r)λ2μMN×((N1)±4μMμFs(νA+ω+μA)(β+μF)N2(MicritMi)rλ2νA2βρϕθK(1r)),F±=rμMM±(1r)(β+μF),Fs±=rμMβ(M±)2μFs[((1r)(β+μF)+rμM)M±+(1r)(β+μF)Mi].

Remark 3.1

Notice that if Mi=Micrit , then E3 and E3+ collapse to an equilibrium E3 with M=θKνAλ(1r)(N1)2μMN,which provides the minimum threshold condition N=1+2(1r)(β+μF)Mi+(1r)(β+μF)Mi[(1r)(β+μF)(1+Mi)+rμM](1r)(β+μF)+rμMfor the existence of the non-trivial equilibria, with N>1.

Remark 3.2

According to Theorem 3.3, the constant release system (Equation11) has a threshold Micrit for the release amount Mi. When the release amount Mi>Micrit thus this leads to an elimination of the population of wild mosquitoes.

3.2.2.2 Stability of the non-trivial equilibria

In this subsection, we assume that N>1 and Mi<Micrit . Thus, we will analyse the stability of the equilibria E3 and E3+. For these points it is clear that μi is an eigenvalue of the Jacobian JE3. To determine the characteristic equation for the other four eigenvalues, we use the following equations obtained at equilibrium of system (Equation11). ρϕFsθK+νA+ω+μA=ρϕβMFμFsθK(M+Mi+F)+νA+ω+μA,ρϕ(1AθK)=μFs(νA+ω+μA)(M+Mi+F)AβMF,νAρϕ(1AθK)βM(M+Mi)(M+Mi+F)2=μFs(β+μF)(νA+ω+μA)(M+Mi)M+Mi+F,νAρϕλ(1r)(1AθK)βF(F+Mi)(M+Mi+F)2=μFsμM(νA+ω+μA)(F+Mi)M+Mi+F.After some manipulation it can be seen that the eigenvalues are the roots of the polynomial (17) Q(x)=x4+b1x3+b2x2+b3x+b4,(17) with b1=ρϕβMFμFsθK(M+Mi+F)+νA+ω+μA+μM+μF+μFs+β,b2=(μM+μF+μFs+β)[ρϕβMFμFsθK(M+Mi+F)+νA+ω+μA]+μM(μF+μFs+β),b3=ρϕβ((μM+μFs)(β+μF)+μMμFs)μFsθK(M+Mi+F)+μM(β+μF)(μM+νA+ω+μA)+μMμFs(νA+ω+μA)(μMM+(β+μF)F)M+Mi+F,b4=rμM2ρϕβ((M)2μFsθK(1r)(νA+ω+μA)(β+μF)MirμMρϕβ)θK(1r)(M+Mi+F).By the Routh-Hurwitz criteria, the roots of a polynomial of order four have negative real parts if and only if bi>0, with i=1,,4 and (b1b2b3)b3>b12b4. It is clear that for the polynomial (Equation17) the coefficients b1,b2 and b3 are positive. Moreover, after some calculations, it can be seen that the last condition is fulfilled for all positive M. Therefore, the stability of both M and M+ is given by the sign of the coefficient b4, which is positive if and only if (18) f(M)=(M)2μFsθK(1r)(νA+ω+μA)(β+μF)MirμMρϕβ>0.(18) Notice that f(M) has a unique positive root M1=μFsθK(1r)(νA+ω+μA)(β+μF)MirμMρϕβ.Now, we evaluate the polynomial (Equation15) at M1. Indeed, we have, g(M1)=2(1r)(β+μF)Mi(1r)(β+μF)+rμM(N1)μFsθK(1r)(νA+ω+μA)(β+μF)MirμMρϕβ<0,according to the biological existence conditions (Equation16). The inequality g(M1)<0 implies M<M1<M+ due to the positively defined coefficient of (M)2. Since f(M) is such that f(M)<0 for M<M1 and f(M)>0 for M>M1. Thus, we have b4<0 for M<M1 and b4>0 for M+>M1. Therefore, E3 is always unstable and E3+ is locally asymptotically stable.

This result is stated below:

Theorem 3.4

Assume that N>1 and Mi<Micrit . The non-trivial equilibrium point E3 is always unstable and the non-trivial equilibrium point E3+ is locally asymptotically stable.

4. Mathematical model of mosquito population dynamics with impulsive release of Wolbachia-infected males

4.1. Model formulation and preliminary results

The most common approach of mosquito release often assumes that the release rate is constant, i.e. that the release of Wolbachia-infected male mosquitoes is uninterrupted. In reality, the mosquito release process is not constant. Each release of Wolbachia-infected male mosquitoes results in an instantaneous change in mosquito numbers, which can result in a pulsed change in the number of Wolbachia-infected male mosquitoes. Consequently, an impulsive system is a more reasonable description of the whole release process. Thus, we assume that Wolbachia-infected male mosquitoes releases are carried out at fixed intervals T and that the quantity of mosquitoes released is the same each time. Model of mosquito dynamics by the release of Wolbachia-infected male mosquitoes, is given by the following impulsive system: (19) {M˙i(t)=μiMi(t),tnT Mi(nT+)=Mi(nT)+Λ,t=nTwith n=1,2,,NT.(19) The last condition implies the release of Wolbachia-infected male mosquitoes every T period and at most NT times, Dumont and Tchuenche [Citation33]. Altogether, this leads to the efficient number of Wolbachia-infected male mosquitoes given by Λ. In addition, nT+ denotes the moment after pulse releases at time nT. The items Mi(nT+) denote the mosquito densities of Wolbachia-infected male mosquitoes ones at time nT after pulse releases.

Based on systems (Equation1), (Equation11) and (Equation19), we give a model for the interaction impulse between Wolbachia-infected male mosquitoes and wild mosquitoes. (20) {A˙(t)=ρϕ(1A(t)θK)Fs(t)(νA+ω+μA)A(t)M˙(t)=(1r)λνAA(t)μMM(t)F˙(t)=νAA(t)(β+μF)F(t)F˙s(t)=βM(t)M(t)+Mi(t)+F(t)F(t)μFsFs(t)M˙i(t)=μiMi(t),}tnTA(nT+)=A(nT),M(nT+)=M(nT),F(nT+)=F(nT),Fs(nT+)=Fs(nT),Mi(nT+)=Mi(nT)+Λ,}t=nTwith n=1,2,,NT.(20) In the formulation of model (Equation20), Wolbachia-infected male mosquitoes are periodically released into the environment (Mi(nT+)=Mi(nT)+Λ), where T is the time elapsed between two successive releases of Wolbachia-infected male mosquitoes and nT+ is the time immediately following the nth release of Wolbachia-infected male mosquitoes. At each instant nT, a constant number of Wolbachia-infected male mosquitoes (Λ) is released into the environment. The items A(nT+),M(nT+),F(nT+),Fs(nT+) denote the mosquito densities of aquatic stages, wild males, young females not yet laying eggs and fertilized females, respectively ones at time nT after pulse releases.

We proceed to the mathematical analysis of model (Equation20). Let X(t)=(A(t),M(t),F(t),Fs(t),Mi(t)) and consider the set, Ω1={X(t)R+5:A(t)θK, M(t)(1r)λνAθKμM, F(t)νAθKβ+μF,Fs(t)νAβθKμFs(β+μF), Mi(t)Λ1eμiT}.Now, we consider the following system (Equation19), which models the periodic release of Wolbachia-infected male mosquitoes. Using two equations of system (Equation19), we get the following result.

Lemma 4.1

The impulsive system (Equation19) has a unique positive periodic solution and Miper(t) is globally asymptotically stable, where (21) Miper(t)=Λeμi(tnT)1eμiT,nT<t(n+1)T.(21)

Proof.

From the first equation of system (Equation19), over the nth impulsive interval it follows that, Mi(t)=Mi(nT+)eμi(tnT),nT<t(n+1)T.The following stroboscopic map can be established using the second equation of system (Equation19), Huang et al. [Citation20], Iboi et al. [Citation48]. Mi[(n+1)T+]=Mi[(n+1)T]+Λ=Mi(nT+).Hence, the unique positive periodic solution of system (Equation19) given by (22) Miper(t)=Λeμi(tnT)1eμiT,nT<t(n+1)T,(22) is globally asymptotically stable.

Taking into account (Equation21), we deduce that solutions of system (Equation20) converge in the sense of L(0,+) norm, to solutions of system (Equation23) [Citation19]. Precisely, only Mi(t) is substituted by Miper(t) in (Equation20) to obtain system (Equation23): (23) {A˙(t)=ρϕ(1A(t)θK)Fs(t)(νA+ω+μA)A(t)M˙(t)=(1r)λνAA(t)μMM(t)F˙(t)=νAA(t)(β+μF)F(t)F˙s(t)=βM(t)M(t)+Miper(t)+F(t)F(t)μFsFs(t),}tnTA(nT+)=A(nT),M(nT+)=M(nT),F(nT+)=F(nT),Fs(nT+)=Fs(nT),}t=nTwith n=1,2,,NT.(23) with initial condition (24) A(0)0,M(0)0,F(0)0,Fs(0)0.(24) System (Equation23) can be written as follows (25) {X˙=A(X,t)X,tnT,X(nT+)=X(nT),t=nT, with n=1,2,,NT.(25) with A(X,t)=((νA+ω+μA)00ρϕ(1A(t)θK)(1r)λνAμM00νA0(β+μF)000βM(t)M(t)+Miper(t)+F(t)μFs)where A(X,t) is a Metzler matrix for all XR+4, and t0.

Remark 4.1

In particular using the fourth equation of system (Equation23), we have F˙s(t)=βM(t)M(t)+Miper(t)+F(t)F(t)μFsFs(t),βM(t)M(t)+F(t)F(t)μFsFs(t),the solutions of system (Equation23) are bounded and above by the solution of system (Equation1), showing that the release of Wolbachia-infected male mosquitoes has a positive effect, i.e. it reduces the number of fertile females, thereby reducing the next generation of new offsprings.

Using the results obtained in (Equation21), we have (26) Mi,minperMiper(t)Mi,maxper,(26) where Mi,minper=Λ1eμiTeμiTandMi,maxper=Λ1eμiT.Using the previous upper and lower bound for Miper(t) , we now consider the following matrix Amin(X,t)=((νA+ω+μA)00ρϕ(1A(t)θK)(1r)λνAμM00νA0(β+μF)000βM(t)M(t)+Mi,minper(t)+F(t)μFs)and Amax(X,t)=((νA+ω+μA)00ρϕ(1A(t)θK)(1r)λνAμM00νA0(β+μF)000βM(t)M(t)+Mi,maxper(t)+F(t)μFs)Then, let us consider the following systems (27) {X_˙=Amin(X_,t)X_,tnT,X_(nT+)=X_(nT),t=nT, with n=1,2,,NT.(27) and (28) {X¯˙=Amax(X¯,t)X¯,tnT,X¯(nT+)=X¯(nT),t=nT, with n=1,2,,NT.(28) with X_(0)=X¯(0)=X0R+4. Amin and Amax being Metzler matrix, we have existence and uniqueness of X_ and X¯ in R+4. We now state and prove the following lemma which will be used in the sequel.

Lemma 4.2

Let systems (Equation23)(Equation27), and (Equation28), have the same initial condition X_(0)=X¯(0)=X0R+4. Then, for all t0, we have X_(t)X(t)X¯(t).

Proof.

Consider Y=XX_. Then, using Equations (Equation25) and (Equation27), we have Y˙=B(t)Y,with B(t)=((νA+ω+μA+ρϕFs(t)θK)00ρϕ(1A_(t)θK)(1r)λνAμM00νA0(β+μF)00b01(t)b02(t)μFs)where b01(t)=βF(t)(F_(t)+Mi,minper(t))(M(t)+Mi,minper(t)+F)(M_(t)+Mi,minper(t)+F_(t)),b02(t)=βM_(t)(M_(t)+Mi,minper(t))(M(t)+Mi,minper(t)+F(t))(M_(t)+Mi,minper(t)+F_(t)).For all t0, B(t) is Metzler, which implies that if Y(0)0, then Y(t)0. Thus, we deduce X(t)X_(t) for all t0. Using the same type of reasoning, we show that X(t)X¯(t) for all t0. This completes the proof.

Theorem 4.3

The solutions of system (Equation20) with initial condition (Equation24) enter eventually into Ω1 and Ω1 is positively invariant for system (Equation20).

Proof.

Let us note that Miper(t) is positive. It should be noted that the right-hand side of each of the four equations of system (Equation23) is continuous and locally-Lipschitz at t = 0. Therefore, a solution of the model with non-negative initial conditions exists and is unique for all time t>0. Moreover, the periodic solution Miper(t) is bounded according to the above. This solution satisfies limt+Miper(t)Λ1eμiT.Using the results obtained in Proposition 2.1, then all solutions of system (Equation20) are bounded for all time t>0.

4.2. The trivial periodic solution and its stability

4.2.1. Existence of trivial periodic solution

Firstly, we discuss the existence of the periodic solution for system (Equation20) when all wild mosquitoes are eliminated. That is to say A(t)=0,M(t)=0,F(t)=0,Fs(t)=0, then system (Equation20) can be simplified to (29) {A˙(t)=0,M˙(t)=0,F˙(t)=0,F˙s(t)=0,M˙i=μiMi,}tnTA(nT+)=A(nT),M(nT+)=M(nT),F(nT+)=F(nT),Fs(nT+)=Fs(nT),Mi(nT+)=Mi(nT)+Λ,}t=nTwith n=1,2,,NT.(29) Thus, system (Equation20) has a trivial periodic solution corresponding to the extinction of wild mosquitoes (30) E5(t)=(A(t),M(t),F(t),Fs(t),Mi(t))=(0,0,0,0,Miper(t))(30) where Miper(t) is given in (Equation21).

4.2.2. The local stability of trivial periodic solutions

Theorem 4.4

The trivial periodic solution E5(t)=(0,0,0,0,Miper(t)) of system (Equation20) is always locally asymptotically stable.

Proof.

In order to determine the local stability of trivial periodic solution E5(t), we define A~(t),M~(t),F~(t), F~s(t) and M~i(t) small amplitude perturbations [Citation20, Citation49–51] of system (Equation20) by {A~(t)=A(t)A(t),M~(t)=M(t)M(t),F~(t)=F(t)F(t),F~s(t)=Fs(t)Fs(t),M~i(t)=Mi(t)Miper(t).The impulse term disappears at this time M~i(nT+)=[Mi(nT)+Λ][M¯i(nT)+Λ]=M~i(nT).The linearization of system (Equation20) is presented by the following system: (31) {A~˙(t)=ρϕ(1A~(t)θK)F~s(t)(νA+ω+μA)A~(t),M~˙(t)=(1r)λνAA~(t)μMM~(t),F~˙(t)=νAA~(t)(β+μF)F~(t),F~˙s(t)=βM~(t)M~(t)+M~i(t)+Miper(t)+F~(t)F~(t)μFsF~s(t),M~˙i=μiM~i,}tnTA~(nT+)=A~(nT),M~(nT+)=M~(nT),F~(nT+)=F~(nT),F~s(nT+)=F~s(nT),M~i(nT+)=M~i(nT),}t=nTwith n=1,2,,NT.(31) The above equation can be written as (32) (A~(t)M~(t)F~(t)F~s(t)M~i(t))=Φ(t)(A~(0)M~(0)F~(0)F~s(0)M~i(0))(32) and the Φ(t) satisfies dΦ(t)dt=U(t)Φ(t)with (33) U(t)=((ρϕF~s(t)θK+νA+ω+μA)00ρϕ(1A~(t)θK)0(1r)λνAμM000νA0(β+μF)000a~1a~2μFsa~30000μi)(33) where (34) {a~1=βF~(t)(F~(t)+M~i(t)+Miper(t))(M~(t)+M~i(t)+Miper(t)+F~(t))2,a~2=βM~(t)(M~(t)+M~i(t)+Miper(t))(M~(t)+M~i(t)+Miper(t)+F~(t))2,a~3=βM~(t)F~(t)(M~(t)+M~i(t)+Miper(t)+F~(t))2,(34) and Φ(0)=I, where I is the identity matrix.

Thus, when A~(t)=0,M~(t)=0,F~(t)=0,F~s(t)=0, it leads to the following linear constant system which is equivalent to Φ˙(t)=U(t)Φ(t)with U(t)=((νA+ω+μA)00ρϕ0(1r)λνAμM000νA0(β+μF)00000μFs00000μi)The resetting impulsive conditions of system (Equation31) become: (35) (A~(nT+)M~(nT+)F~(nT+)F~s(nT+)M~i(nT+))=(1000001000001000001000001)(A~(nT)M~(nT)F~(nT)F~s(nT)M~i(nT))(35) The monodromy matrix [Citation30, Citation52, Citation53] is Γ=(1000001000001000001000001)Φ(T)=Φ(T)It can be seen that Φ(T)=Φ(0)exp(0Ttr(U(t))dt),and Φ(0) is the identity matrix. Let θ1,θ2,θ3,θ4,θ5, the eigenvalues of exp(0Ttr(U(t))dt) then θ1=exp(0T(νA+ω+μA)dt),θ2=exp(0TμMdt),θ3=exp(0T(β+μF)dt),θ4=exp(0TμFsdt),θ5=exp(0Tμidt).Obviously, θ1<1,θ2<1,θ3<1,θ4<1,θ5<1. Using Floquet's theorem [Citation33, Citation51], we deduce that the trivial periodic solution E5(t)=(0,0,0,0,Miper(t)) of system (Equation20) is always locally asymptotically stable. Whatever the value taken by the basic offspring number N, the previous result holds.

4.2.3. The global stability of trivial periodic solutions

In this part, we will further explore the global stability of the trivial periodic solutions for system (Equation20) by using Lyapunov function [Citation51, Citation54].

Theorem 4.5

For system (Equation20), when N<1, then the trivial periodic solution E5(t)=(0,0,0,0,Miper(t)) is globally asymptotically stable.

Proof.

Consider the Lyapunov function defined in (Equation14) by V=μFs[(1r)(β+μF)+rμM](1r)rρϕβA(t)+M(t)+F(t)+(1r)(β+μF)+rμM(1r)rβFs(t)and using the proof of Theorem 3.2, we have V˙(1r)rρϕβμFs[(1r)(β+μF)+rμM](νA+ω+μA)(N1)A(t).Thus V˙(t)0 if N1. Note that V(t)=0 if and only if A(t)=0,M(t)=0,F(t)=0,Fs(t)=0, and we know that the periodic solution Miper(t) is globally asymptotically stable for the dynamics of Wolbachia-infected male mosquitoes. Then, the trivial periodic solution of system (Equation20) is globally attractive provided N<1. We have proved that the trivial periodic solution E5(t) is always locally asymptotically stable in Theorem 4.4. Therefore, the trivial periodic solution E5(t) is globally asymptotically stable when N<1.

4.3. The non-trivial periodic solution and its stability

4.3.1. Existence of non-trivial periodic solution

In practice sometimes we can not completely eradicate mosquito populations, but rather keep them within a certain range. With this aim, we study the non-trivial periodic solution of system (Equation20) corresponding to the partial survival of wild mosquitoes, using similar reasoning in [Citation50, Citation51]. In the following, we continue to use the positive equilibrium of system (Equation31) point to discuss the existence of non-trivial periodic solutions of system (Equation20). System (Equation31) satisfies at the positive equilibrium point (36) {ρϕ(1A~(t)θK)F~s(t)=(νA+ω+μA)A~(t),(1r)λνAA~(t)=μMM~(t),νAA~(t)=(β+μF)F~(t),βM~(t)M~(t)+Miper(t)+F~(t)F~(t)=μFsF~s(t),M~i=0.(36) By solving Equation (Equation36), we obtain the following equation as a function of M~(t), which is given as follows: (37) μMN(1r)λνAθK(M~(t))2(N1)M~(t)+(1r)(β+μF)(1r)(β+μF)+rμMMiper(t)=0(37) From Equation (Equation37), we know μMN(1r)λνAθK>0and(1r)(β+μF)(1r)(β+μF)+rμMMiper(t)>0,then existence of the positive equilibrium of system (Equation31) depends on the sign of N. If N1, then there is no positive equilibrium point for system (Equation31), corresponding to the absence of positive non-trivial periodic solutions for system (Equation20). If N>1 then there exist 0, 1 or 2 positive equilibria for system (Equation31), corresponding to the existence of 0, 1 or 2 positive equilibria for system (Equation20). We assume that N>1, so the number of positive roots in Equation (Equation37) depends on (38) Δ(Miper)=(N1)24μM(1r)(β+μF)NMiper(t)νAθKλ(1r)((1r)(β+μF)+rμM).(38) If Δ(Miper)>0 then there are two different positive roots of Equation (Equation37), and if Δ(Miper)<0 then there are no positive roots of Equation (Equation37). Since the value of Δ(Miper) depends on Miper(t) which varies periodically, the value of Δ(Miper) cannot always be equal to 0. Thus, let us define a release threshold for the number of positive periodic solutions by (39) Micrit=rλ2νA2βρϕθK(1r)4μMμFs(νA+ω+μA)(β+μF)(11N)2.(39) Inegality (Equation26) shows that Miper(t) , has minimum value Mi,minper and maximum value Mi,maxper depend on release period T and release amount Λ. When N>1, therefore, there is a threshold Micrit>0 that makes Δ(Miper)<0 or Δ(Miper)>0 when Mi,minper>Micrit or Mi,maxper<Micrit , respectively. So, system (Equation20) have only no non-trivial periodic solutions when Mi,minper>Micrit and two positive non-trivial periodic solutions when Mi,maxper<Micrit.

Thus, the results for the existence of non-trivial periodic solution for system (Equation20) are summarized as follows.

Theorem 4.6

(i)

When N1, then system (Equation20) has no positive non-trivial periodic solution.

(ii)

When N>1 and Mi,minper>Micrit , then system (Equation20) still has no positive non-trivial periodic solution

(iii)

When N>1 and Mi,maxper<Micrit , then system (Equation20) has two positive non-trivial periodic solutions.

Remark 4.2

Considering system (Equation27), we can define the following release thresholds that ensures the existence of a non-trivial equilibrium for all T>0, (40) Λ=Micrit(eμiT1),(40) where Micrit is the threshold given in (Equation39). From equality (Equation40), we immediately deduce the following result. If T>0 is chosen such that T>T then Λ<Λ, with T=1μilog(1+ΛMicrit)andΛ=Micrit(eμiT1)where T represents the release period threshold and Λ represents the release threshold.

Thus, the conclusion of the non-trivial positive periodic solution of system (Equation20) can be summarized as follows.

Theorem 4.7

(i)

If N1, then system (Equation20) has no positive non-trivial periodic solution.

(ii)

If N>1,T<T and Λ>Λ, then system (Equation20) still has no positive non-trivial periodic solution.

(iii)

If N>1,T>T and Λ<Λ, then system (Equation20) has two positive non-trivial periodic solution.

4.3.2. The local stability of non-trivial periodic solution

In the following, we consider condition (iii) of Theorem 4.7 and discuss the stability of non-trivial periodic solutions using Floquet theory. Under these conditions, the perturbation system (Equation31) has two positive two positive equilibra points, denoted by E6±(t)=(A~±(t),M~±(t),F~±(t),F~s±(t),M~i(t))and satisfies M~i(t)=0,where {A~±(t)=μMM~±(t)(1r)λνA,M~±(t)=νAθKλ(1r)(N1)2μMN(1±Δ(Miper)),F~±(t)=rμMM~±(t)(1r)(β+μF),F~s±(t)=rμMβ(M~±(t))2μFs[((1r)(β+μF)+rμM)M~±(t)+(1r)(β+μF)Miper(t)],corresponding to system (Equation20), there are two positive non-trivial periodic solutions E±(t)=(A±(t),M±(t),F±(t),Fs±(t),Miper(t)).To study the local stability of non-trivial periodic solutions, we use a proof similar to Theorem 4.4. Thus, we consider the matrix U(t) defined at (Equation33) and we evaluate at the point E6±(t)=(A~±(t),M~±(t),F~±(t),F~s±(t),M~i(t)). At this time (41) UE6(t)=((ρϕF~s(t)θK+νA+ω+μA)00ρϕ(1A~(t)θK)0(1r)λνAμM000νA0(β+μF)000a~4a~5μFsa~60000μi)(41) where (42) {a~4=βF~(t)(F~(t)+Miper(t))(M~(t)+Miper(t)+F~(t))2,a~5=βM~(t)(M~(t)+Miper(t))(M~(t)+Miper(t)+F~(t))2,a~6=βM~(t)F~(t)(M~(t)+Miper(t)+F~(t))2,(42) We consider that the eigenvalues of exp(0TUE6(t)dt) are δ1,δ2,δ3,δ4,δ5 and the eigenvalues of UE6(t) are δ1,δ2,δ3,δ4,δ5, then the characteristic polynomial of matrix (Equation41) is (43) (δμi)(δ4+b~1δ3+b~2δ2+b~3δ+b~4)=0,(43) with b~1=ρϕβM~(t)F~(t)μFsθK(M~(t)+Miper(t)+F~(t))+νA+ω+μA+μM+μF+μFs+β,b~2=(μM+μF+μFs+β)[ρϕβM~(t)F~(t)μFsθK(M~(t)+Miper(t)+F~(t))+νA+ω+μA]+μM(μF+μFs+β),b~3=ρϕβ((μM+μFs)(β+μF)+μMμFs)μFsθK(M~(t)+Miper(t)+F~(t))+μM(β+μF)(μM+νA+ω+μA)+μMμFs(νA+ω+μA)(μMM~(t)+(β+μF)F~(t))M~(t)+Miper(t)+F~(t),b~4=rμM2ρϕβ((M~(t))2μFsθK(1r)(νA+ω+μA)(β+μF)Miper(t)rμMρϕβ)θK(1r)(M~(t)+Miper(t)+F~(t)).We recall some of the expressions obtained in Section 3.2.2.2 useful for the rest of our demonstration. f(M~(t))=(M~(t))2μFsθK(1r)(νA+ω+μA)(β+μF)Miper(t)rμMρϕβ,M~1(t)=μFsθK(1r)(νA+ω+μA)(β+μF)Miper(t)rμMρϕβ,g(M~1(t))=2(1r)(β+μF)Miper(t)(1r)(β+μF)+rμM(N1)M~1(t).It should be noted that the eigenvalues of exp(0TUE6(t)dt) are δ1=exp(0Tδ1dt),δ2=exp(0Tδ2dt),δ3=exp(0Tδ3dt),δ4=exp(0Tδ4dt),δ5=exp(0Tδ5dt).We have δ5=μi<0, therefore δ5<1. Then, using the results obtained in Section 3.2.2.2, we can conclude that the signs of δ1,δ2,δ3,δ4 depend on b~4(M~(t)). Consider the two positive equilibria E6(t)=(A~(t),M~(t),F~(t),F~s(t),0) and E6+(t)=(A~+(t),M~+(t),F~+(t),F~s+(t),0) of perturbation system (Equation31). By using the inequality defined in Section 3.2.2.2 (i.e. g(M~1(t)<0)), we have M~(t)<M~1(t)<M~+(t). Since f(M~(t)) is such that f(M~(t))<0 for M~(t)<M~1(t) and f(M~(t))>0 for M~(t)>M~1(t), we deduce that b~4(M~(t))<0 for M~(t)<M~1(t) and b~4(M~+(t))>0 for M~+(t)>M~1(t). Therefore, at E6, we have b~4(M~+(t))>0 the roots of Equation (Equation43) have negative real parts, which leads to δ1,δ2,δ3 and δ4 all being less than 1. Moreover, at E6, we have b~4(M~(t))<0, the roots of the Equation (Equation43) have at least one positive real part, which leads at δ1,δ2,δ3 or δ4 greater than 1. According to the Floquet theory [Citation33, Citation51] and using the fact that Miper(t) is globally asymptotically stable, we deduce that the non-trivial periodic solution E(t)=(A(t),M(t),F(t),Fs(t),Miper(t)) is unstable and E+(t)=(A+(t),M+(t),F+(t),Fs+(t),Miper(t)) is locally asymptotically stable.

We summarize the stability of non-trivial periodic solution of system (Equation20) as follows:

Theorem 4.8

Assumes that N>1, T>T and Λ<Λ. The non-trivial periodic solution E(t)=(A(t),M(t),F(t),Fs(t),Miper(t)) of system (Equation20) is unstable and the non-trivial periodic solution E+(t)=(A+(t),M+(t),F+(t),Fs+(t),Miper(t)) is locally asymptotically stable.

Remark 4.3

The above results show that for wild mosquitoes to persist in the wild, the condition N>1 is necessary. However, when Wolbachia-infected male mosquitoes are released, this condition may not be sufficient, as some of the females may not actually be fertilized. According to Remark 4.2 and Theorem 4.8, when the conditions Λ>Λ and T<T are verified, the wild mosquito population can be successfully suppressed.

5. Numerical simulation and discussion

In this section, we perform some numerical simulations to illustrate the theoretical results established previously. Our numerical simulation has been performed using the MATLAB technical computing software with the fourth order Runge–Kutta method. The values of the parameters are given in Table  and the range of values for the control parameters in Table .

Table 3. Parameters values of the mosquito models (Equation1), (Equation10) and (Equation20).

Table 4. Range of values for the control parameters.

5.1. Numerical simulation of wild mosquito population dynamics

In this section, we illustrate different possible scenarios for reducing the population of wild mosquitoes. More specifically, we considered the following cases:

  1. different values for larvicide treatment,

  2. different values for mechanical control,

  3. the influence of a combination of mechanical control and larvicide treatment on the growth dynamics of wild mosquitoes.

5.1.1. Dynamics of wild mosquitoes when N<1 and N>1

In this subsection, we simulate the wild mosquitoes model (Equation1) with different initial conditions (Figure ). Figure shows the trajectories of model (Equation1) when N<1. These trajectories converge towards the trivial equilibrium E0. This means that the mosquito population disappears in the environment.

Figure 3. Simulation of model (Equation1) using various initial condition when ϕ=8,K=50,000,θ=1,νA=0.003,r=0.5,μA=0.08,ω=0,μM=1/10,β=0.1, μF=1/7, ρ=1, λ=1 and μFs=1/10 this leads to N=0.4217<1. (a) Aquatic. (b) Males. (c) Young females and (d) Fertilized females.

Figure 3. Simulation of model (Equation1(1) {A˙(t)=ρϕ(1−A(t)θK)Fs(t)−(νA+ω+μA)A(t),M˙(t)=(1−r)λνAA(t)−μMM(t),F˙(t)=rλνAA(t)−(β+μF)F(t),F˙s(t)=βM(t)M(t)+F(t)F(t)−μFsFs(t),(1) ) using various initial condition when ϕ=8,K=50,000,θ=1,νA=0.003,r=0.5,μA=0.08,ω=0,μM=1/10,β=0.1, μF=1/7, ρ=1, λ=1 and μFs=1/10 this leads to N=0.4217<1. (a) Aquatic. (b) Males. (c) Young females and (d) Fertilized females.

Figure  shows the trajectories of model (Equation1) when N>1. These trajectories converge towards the unique non-trivial equilibrium E1. Hence, mosquito populations persist in the environment. These numerical verifications illustrate the result stated in the Theorem 2.5.

5.1.2. Dynamics of wild mosquitoes with larvicide control only

The most common larviciding treatment is the use of temephos (Abate) and Bacillus thuringiensis israelensis [Citation34, Citation48]. Here, we have assumed that larvicide spraying takes place at least once every 15 days over a period of 150 days. For this, we fixed θ=1 and consider ϕ=10,K=50,000,νA=0.1,r=0.5,μA=0.0583,μM=1/10,β=0.7, μF=1/7, ρ=1, λ=1 and μFs=1/10. The simulations showing the effectiveness of the larvicide treatment are presented in Figure . We note that the population of mosquitoes in the aquatic stage is more affected than those in the aldute stage as shown in Table . We also observe that the effectiveness of larvicidal treatments has an impact on the basic offspring number, N, more precisely it decreases when the larvicidal control parameter increases.

Figure 4. Simulation of model (Equation1) using various initial condition when ϕ=10,K=50,000,θ=1,νA=0.1,r=0.5,μA=0.0583,ω=0,μM=1/10,β=0.7, μF=1/7, ρ=1, λ=1 and μFs=1/10 this leads to N=23.4499>1. (a) Aquatic. (b) Males. (c) Young females and (d) Fertilized females.

Figure 4. Simulation of model (Equation1(1) {A˙(t)=ρϕ(1−A(t)θK)Fs(t)−(νA+ω+μA)A(t),M˙(t)=(1−r)λνAA(t)−μMM(t),F˙(t)=rλνAA(t)−(β+μF)F(t),F˙s(t)=βM(t)M(t)+F(t)F(t)−μFsFs(t),(1) ) using various initial condition when ϕ=10,K=50,000,θ=1,νA=0.1,r=0.5,μA=0.0583,ω=0,μM=1/10,β=0.7, μF=1/7, ρ=1, λ=1 and μFs=1/10 this leads to N=23.4499>1. (a) Aquatic. (b) Males. (c) Young females and (d) Fertilized females.

Figure 5. Dynamics of mosquito population using different values of larvicide control. (a) Aquatic. (b) Males. (c) Young females and (d) Ferlized females.

Figure 5. Dynamics of mosquito population using different values of larvicide control. (a) Aquatic. (b) Males. (c) Young females and (d) Ferlized females.

Table 5. Non-trivial equilibrium E1=(A,M,F,Fs) of wild mosquitoes with different larvicide control values.

5.1.3. Dynamics of wild mosquitoes with mechanical control only

Mechanical control (mc) is a type of intervention that involves destroying mosquito breeding sites to slow the growth of the immature form of the mosquito population. It has been shown that mechanical control (mc) can be an effective control tool [Citation34, Citation57]. It is an inexpensive form of control, but requires the support of the local population. We assume that mechanical control (mc) leads to a decrease in K, i.e. a decrease in the equilibrium of the wild aquatic stage A. From relationship (Equation8), we deduce that reducing A for mc% corresponds to a decrease in K as follows K=NA(N1)(1mc100),with θ=1mc100.We assume that the effectiveness of mechanical control does not depend on the seasons and should be applied as soon as possible. Thus, θ is defined as being between (0,1], and the numerical simulation is carried out for 150 days. For this, we fixed ω=0 and consider ϕ=10,K=50,000,νA=0.1,r=0.5,μA=0.0583,μM=1/10,β=0.7, μF=1/7,μFs=1/10, ρ=1, λ=1. We use the following initial conditions: A(0)=1000,M(0)=1000,F(0)=1000,Fs(0)=1000. The results of the numerical simulation for the effectiveness of mechanical control is shown in Figure . We observe that the mosquito population decreases when the level of mechanical control parameter increases. These figures show that this method of mechanical control is appropriate to fight against the proliferation of mosquitoes, especially in their immature stage as indicated in Table .

Table 6. Non-trivial equilibrium E1=(A,M,F,Fs) of wild mosquitoes with different mechanic control values.

5.1.4. Dynamics of wild mosquitoes with mechanical and larvicide control

We varied parameter θ,ω and fixed ϕ=10,K=50,000,νA=0.1,r=0.5,μA=0.0583,μM=1/10,β=0.7, μF=1/7,μFs=1/10, ρ=1, λ=1. Figure  shows the impact of the combination of mechanical and larvicide control. Moreover, these figures show interesting results that are better than mechanical control in Figure and larviciding in Figure . We observe that when larval capacity is reduced by 50% (θ=0.5) with low larviciding control (ω=0.05), the mosquito population decreases to a significantly lower level than with high larviciding control (ω=0.25) and low mechanical control 10% (θ=0.90) as show in Table . Thus, this combined control could have a greater impact on reducing the growth dynamics of wild mosquitoes.

Figure 6. Dynamics of mosquito population using different values of mechanical and larvicide controls. (a) Aquatic. (b) Males. (c) Young females and (d) Ferlized females.

Figure 6. Dynamics of mosquito population using different values of mechanical and larvicide controls. (a) Aquatic. (b) Males. (c) Young females and (d) Ferlized females.

Figure 7. Dynamics of mosquito population using different values of mechanical control. (a) Aquatic. (b) Males. (c) Young females and (d) Ferlized females.

Figure 7. Dynamics of mosquito population using different values of mechanical control. (a) Aquatic. (b) Males. (c) Young females and (d) Ferlized females.

Table 7. Non-trivial equilibrium E1=(A,M,F,Fs) of wild mosquitoes with different mechanic and larvicide control values.

5.2. Numerical simulation of mosquito dynamics with constant releases of Wolbachia-infected male mosquitoes

In this section, we test different possible scenarios to reduce or eradicate the population of wild mosquitoes. Specifically, we have considered the following cases:

  1. different release quantities of Wolbachia-infected male mosquitoes,

  2. the influence of a combination of Wolbachia-infected male mosquitoes release, mechanical control and larvicide treatment on the growth dynamics of wild mosquitoes.

5.2.1. Dynamics of mosquito population with constant release without mechanical and larvicide control

We evaluate in this subsection the impact of different release rates of Wolbachia-infected male mosquitoes on wild mosquito populations. We use the following initial conditions: A(0)=1000,M(0)=2000,F(0)=500,Fs(0)=1500,Fst(0)=1. The parameters used are: ϕ=10,K=50,000,θ=1,νA=0.1,r=0.5,μA=0.0583,ω=0,μM=1/10,β=0.7, μF=1/7,μFs=1/10,μi=0.12,μst=1/10, ρ=1, λ=1 and the following release rates are assumed: Λ=10,000,12,500,15,000,17,500 which verifies the condition Mi[0,Micrit) as indicated in Table . Thus, according to Theorem 3.3, system (Equation11) admits two positive equilibra points E3 given in Table  and E3+ given in Table . By verifying the stability conditions given in Table , we conclude that equilibrium E3 is unstable and equilibrium E3+ is locally stable, so the model solutions approach equilibrium E3+ as shown in Figure  which confirms our analytical results stated in Theorem 3.4.

Figure 8. The effect of the release rate of Wolbachia-infected male mosquitoes without mechanical control and larvicide treatment on the population of wild mosquitoes with Mi(0,Micrit) as indicated in Table . (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females and (e) Sterile females.

Figure 8. The effect of the release rate of Wolbachia-infected male mosquitoes without mechanical control and larvicide treatment on the population of wild mosquitoes with Mi∗∈(0,Micrit) as indicated in Table 8. (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females and (e) Sterile females.

Table 8. Equilibrium point E3+ for different values of Wolbachia-infected male releases Mi(0,Micrit) with % mc=0 and ω=0.

Table 9. Equilibrium point E3 for different values ofWolbachia-infected male releases Mi(0,Micrit) with % mc=0 and ω=0.

Table 10. Numerical values verifying the condition of stability of the non-trivial equilibrium for the case of Figure .

Figure  shows the total number of mosquitoes in the aquatic stage, wild males, young females, fertile females and the number of Wolbachia-infected female mosquitoes. When the release rate is increased from 10,000 to 17,500, we observe that the total number of mosquitoes in the aquatic stage, wild male, young female and fertile female stages decreases to a significantly lower level. Thus, when a large number of Wolbachia-infected male mosquitoes are released, this will result in a significant reduction in the wild mosquito population.

5.2.2. Dynamics of mosquito population with constant release, mechanical and larvicide control

In this subsection, we evaluate the constant release of Wolbachia-infected males in combination with mechanical control and larvicide treatment.

Figure  shows the effect of the release rate of Wolbachia-infected male mosquitoes combined mechanical control with larvicide treatment on the population of wild mosquitoes with Mi(0,Micrit) as indicated in Table . By comparing the numerical results in Tables  and , we find that a release of Wolbachia-infected males mosquitoes 10,000 with a mechanical control of 20% and a larvicide treatment of 0.05 gives better results than a release of 17,500 Wolbachia-infected males mosquitoes without control measures.

5.3. Numerical simulation of mosquito dynamics with impulsive releases of Wolbachia-infected male mosquitoes

In this section, we simulate different possible scenarios to reduce or eradicate the population of wild mosquitoes by the impulsive release of Wolbachia-infected males. More specifically, we considered the following cases:

  1. different release periodicities,

  2. different quantities of Wolbachia-infected male mosquitoes for each fixed period,

  3. the influence of a combination of impulsive release of Wolbachia-infected male mosquitoes, mechanical control and larvicide treatment on the growth dynamics of wild mosquitoes.

Figure 9. The effect of the release rate of Wolbachia-infected male mosquitoes combined mechanical control and larvicide treatment on the population of wild mosquitoes with Mi(0,Micrit) as indicated in Table . (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females and (e) Sterile females.

Figure 9. The effect of the release rate of Wolbachia-infected male mosquitoes combined mechanical control and larvicide treatment on the population of wild mosquitoes with Mi∗∈(0,Micrit) as indicated in Table 11. (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females and (e) Sterile females.

Figure 10. Dynamics of impulsive release of Wolbachia-infected male on the mosquito population when T=14,Λ=50,000,150,000 which verifies the conditions given in Table . (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females. (e) Sterile females and (f) Sterile males.

Figure 10. Dynamics of impulsive release of Wolbachia-infected male on the mosquito population when T=14,Λ=50,000,150,000 which verifies the conditions given in Table 12. (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females. (e) Sterile females and (f) Sterile males.

Table 11. Impact of the constant release of Wolbachia-infected males combined with mechanical control and larvicidal treatment on mosquito dynamics with Λ=10,000.

Table 12. Numerical values verifying the condition of existence and stability of the non-trivial periodic solution, (i.e. N>1, T>T and Λ<Λ) for the case of Figure .

5.3.1. Dynamics of mosquito population with impulsive releases without mechanical and larvicide control

In this subsection to illustrate our theoretical results, we use the following values θ=1,ω=0,ϕ=10,K=50,000,νA=0.1,r=0.5,μA=0.0583,μM=1/10,β=0.7, μF=1/7,μFs=1/10,μi=0.12,μst=1/10, ρ=1, λ=1.

In the case of Figure , we choose a pulse period T = 14 days and the number of Wolbachia-infected male mosquitoes release Λ=50,000 and 150,000. Moreover, in the case of Figure , we take a pulse period T = 21 days and the number of Wolbachia-infected male mosquitoes release Λ=50,000 and 150,000. Thus, using these parameters, we have N>1, T>T and Λ<Λ, as indicated by the numerical values obtained in Tables  and . Then, according to Theorem 4.7, system (Equation20) has two positive non-trivial periodic solutions E(t)=(A(t),M(t),F(t),Fs(t),Miper(t)) and E+(t)=(A+(t),M+(t),F+(t),Fs+(t),Miper(t)) . Moreover, E(t)=(A(t),M(t),F(t),Fs(t),Miper(t)) is an unstable non-trivial periodic solution and E+(t)=(A+(t),M+(t),F+(t),Fs+(t),Miper(t)) is a locally asymptotically stable non-trivial periodic solution. Additionally, the range of existence of a stable non-trivial periodic solution is given in Tables , , and . We observe in the Figures and that the number of mosquitoes in the aquatic stage, wild males, young females and fertile females decreases as release pulses increase, while Wolbachia-infected female mosquitoes increase, leading to a reduction in the wild mosquito population.

Remark 5.1

From the numerical results obtained in Tables , , , and Figures , we can deduce that a short pulse period and a higher number of releases of Wolbachia-infected male mosquitoes are the best way to quickly reduce wild mosquito population.

Figure 11. Dynamics of impulsive release of Wolbachia-infected male on the mosquito population when T=21,Λ=50,000,150,000 which verifies the conditions given in Table . (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females. (e) Sterile females and (f) Sterile males.

Figure 11. Dynamics of impulsive release of Wolbachia-infected male on the mosquito population when T=21,Λ=50,000,150,000 which verifies the conditions given in Table 15. (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females. (e) Sterile females and (f) Sterile males.

Table 13. The range of existence of a stable non-trivial periodic solution when T = 14 and Λ=150,000 for Figure (blue colour).

Table 14. The range of existence of a stable non-trivial periodic solution when T = 14 and Λ=50,000 for Figure (red colour).

Table 15. Numerical values verifying the condition of existence and stability of the non-trivial periodic solution (i.e. N>1, T>T and Λ<Λ) for the case of Figure .

Table 16. The range of existence of a stable non-trivial periodic solution when T = 21 and Λ=150,000 for Figure (blue colour).

Table 17. The range of existence of a stable non-trivial periodic solution when T = 21 and Λ=50,000 for Figure (red colour).

5.3.2. Dynamics of mosquito population with impulsive releases, mechanical and larvicide control

Now, we evaluate impulsive release of Wolbachia-infected males combined with mechanical control and larvicide treatment.

In Figure (resp. Figure ), we choose an impulse period T = 14 days (resp. T = 21 days), Λ=50,000, θ=0.70 and ω=0.1 (curves in red), then we consider Λ=150,000, θ=1 and ω=0 (curves in blue). Thus, using the parameters given in Table , we have N>1, T>T and Λ<Λ, as indicated by the numerical values obtained in Table  (resp. Table ). Moreover, using the same interpretation argument from Section 5.3.1 we can deduce that the non-trivial periodic solution E+(t)=(A+(t),M+(t),F+(t),Fs+(t),Miper(t)) is locally asymptotically stable. In Figures and , we observe that the release of 50,000 Wolbachia-infected males combined with a mechanical control of 30% (θ=0.70) and a larvicidal treatment that eliminates the population of aquatic stage mosquitoes by 10% (ω=0.1) is significantly better than the release of 150,000 Wolbachia-infected males without a mechanical control and larvicide treatment.

Remark 5.2

The numerical results obtained in Tables , , , , , and Figures , , , allow us to deduce that a long pulse period with a lower number of releases of Wolbachia-infected male mosquitoes combined with mechanical control and larviciding is a better way of rapidly reducing the wild mosquito population than a huge number of male releases over a short period.

Figure 12. Dynamics of the impulsive release of a Wolbachia-infected male on the mosquito population with mechanical and larvicide control when T=14,Λ=50,000,150,000 verifying the conditions given in Table . (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females. (e) Sterile females and (f) Sterile males.

Figure 12. Dynamics of the impulsive release of a Wolbachia-infected male on the mosquito population with mechanical and larvicide control when T=14,Λ=50,000,150,000 verifying the conditions given in Table 18. (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females. (e) Sterile females and (f) Sterile males.

Figure 13. Dynamics of the impulsive release of a Wolbachia-infected male on the mosquito population with mechanical and larvicide control when T=21,Λ=50,000,150,000 verifying the conditions given in Table . (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females. (e) Sterile females and (f) Sterile males.

Figure 13. Dynamics of the impulsive release of a Wolbachia-infected male on the mosquito population with mechanical and larvicide control when T=21,Λ=50,000,150,000 verifying the conditions given in Table 19. (a) Aquatic. (b) Males. (c) Young females. (d) Ferlized females. (e) Sterile females and (f) Sterile males.

Table 18. Numerical values verifying the condition of existence and stability of the non-trivial periodic solution, (i.e. N>1, T>T and Λ<Λ) for the case of Figure .

Table 19. Numerical values verifying the condition of existence and stability of the non-trivial periodic solution, (i.e. N>1, T>T and Λ<Λ) for the case of Figure .

6. Conclusion

Due to the absence of effective clinical treatments and vaccines, current methods of controlling vector-borne diseases are mainly aimed at eradicating the main mosquito vectors. However, traditional means, including larval source reduction and insecticide application, are not sufficient to keep vector population densities below the epidemic risk threshold. In recent years, a very promising strategy has been to release Wolbachia-infected male mosquitoes into wild areas in order to sterilize wild female mosquitoes through cytoplasmic incompatibility.

In this work, we first presented a mathematical model (Equation1) describing the dynamics of the wild mosquito population on the basis of its life cycle. Measures in terms of mechanical control and larviciding were considered to reduce the growth of wild mosquitoes. We have shown that the problem is well posed mathematically and biologically. We defined the basic offspring number N. Specifically, we have shown that, depending on the basic offspring number N, there is an equilibrium point corresponding to the extinction of the wild mosquito population and another corresponding to the persistence of the mosquito population in the environment. We have also shown that when N1, the trivial equilibrium is globally asymptotically stable and when N>1, there exists a unique non-trivial equilibrium that is globally asymptotically stable. Numerical simulations have been presented to illustrate and support the analytical results. Our numerical result suggests that mechanical control combined with mechanical treatment reduces significantly the wild mosquito population.

Next, we added a compartment of Wolbachia-infected male mosquitoes to this initial model (Equation1). In order to evaluate the effects of constant releases of Wolbachia-infected male mosquitoes on the wild mosquito population. In addition, we considered mechanical control and larviciding to assess its impact on mosquito dynamics. We obtained a critical threshold number of releases of Wolbachia-infected male mosquitoes, Micrit . Then, we established theoretically and numerically that the elimination of wild mosquitoes is possible when the number of Wolbachia-infected male mosquitoes released Mi is greater than this threshold. In addition, below this threshold, we have determined properties of a global nature for system equilibria. We prove that when Mi[0,Micrit) , the model admits exactly two equilibrium points; one of which is asymptotically stable and the other unstable. The trivial equilibrium, corresponding to the elimination of wild mosquitoes, is locally asymptotically stable, and becomes globally asymptotically stable when N1. Our numerical result suggests that releasingWolbachia-infected males in combination with mechanical control and larviciding is more efficient, as shown in Figure , since it reduces the number of Wolbachia-infected male mosquitoes initial to be released, which is economically beneficial.

Finally, we have established an impulse model of the interaction between Wolbachia-infected male mosquitoes and wild mosquito population. Firstly, it is found that there is always a trivial periodic solution corresponding to the extinction of mosquito populations, and the corresponding local and global stability conditions are given using Floquet theory and Lyapunov's stability theorem, respectively. Moreover, the results show that the impulsive release of Wolbachia-infected male mosquitoes can lead to the extinction of wild mosquitoes when N>1, T<T and Λ>Λ. Secondly, we find that there can be two solutions in the system when certain threshold conditions are met. The existence of non-trivial periodic solutions is then obtained, and its local stability is discussed when N>1, T>T and Λ<Λ. The analysis shows that there is only one non-trivial periodic solution which is locally asymptotically stable. The numerical results of the impulsive model show that a long pulse period with a lower number of releases of Wolbachia-infected male mosquitoes combined with mechanical control and larviciding Figures and is a better way of rapidly reducing the wild mosquito population than a huge number of male releases over a short period Figures and .

Our results suggest that the release ofWolbachia-infected male mosquitoes, combined with mechanical control and larviciding, is one of the best ways of significantly reducing the wild mosquito population to a very low level. Indeed, control measures combined with releases reduce the cost of producing Wolbachia-infected male mosquitoes in the laboratory.

In our future work, we intend to use the optimal control approach to minimize the costs associated with controlling the mosquito population through the application of larvicides, by introducing the release of Wolbachia-infected male mosquitoes into the environment. In addition, we aim to couple these models to the transmission of malaria and other vector-borne diseases, in order to assess their impact on transmission.

Disclosure statement

No potential conflict of interest was reported by the author(s).

References

  • Caraballo H, King K. Emergency department management of mosquito-borne illness: malaria, dengue, and West Nile virus. Emerg Med Pract. 2014;16(5):1–23.
  • Tolle M. Mosquito-borne diseases. Curr Probl Pediatr Adolesc Health Care. 2009;39:97–140.
  • Hemingway J, Ranson H. Insecticide resistance in insect vectors of human disease. Annu Rev Entomol. 2000;45:371–91. doi: 10.1146/ento.2000.45.issue-1
  • Saha S, Dutta P, Samanta G. Dynamical behavior of SIRS model incorporating government action and public response in presence of deterministic and fluctuating environments. Chaos Solit Fractals. 2022;164:Article ID 112643. doi: 10.1016/j.chaos.2022.112643
  • Traoré B, Sangaré B, Traoré S. A mathematical model of malaria transmission in a periodic environment. J Biol Dyn. 2018;12:400–432. doi: 10.1080/17513758.2018.1468935
  • Engelstädter J, Telschow A. Cytoplasmic incompatibility and host population structure. J Hered. 2009;103:196–207. doi: 10.1038/hdy.2009.53
  • Zheng B, Liu X, Tang M, et al. Use of age-stage structural models to seek optimal Wolbachia-infected male mosquito releases for mosquito-borne disease control. J Theor Biol. 2019;472:95–109. doi: 10.1016/j.jtbi.2019.04.010
  • Bian G, Joshi D, Dong Y, et al. Wolbachia invades Anopheles stephensi populations and induces refractoriness to plasmodium infection. Science. 2013;340:748–751. doi: 10.1126/science.1236192
  • Xianghong Z, Sanyi T, Robert AC. Models to assess how best to replace dengue virus vectors with Wolbachia-infected mosquito populations. Math Biosci. 2015;269:164–177. doi: 10.1016/j.mbs.2015.09.004
  • Dutta P, Samanta G, Nieto JJ. Periodic transmission and vaccination effects in epidemic dynamics: a study using the SIVIS model. Nonlinear Dyn. 2024;112:2381–2409. doi: 10.1007/s11071-023-09157-4
  • Koutou O, Traoré B, Sangaré B. Mathematical modeling of malaria transmission global dynamics: taking into account the immature stages of the vectors. Adv Differ Equ. 2018;2018(220):1–34.
  • Saha S, Samanta G. Analysis of a host–vector dynamics of a dengue disease model with optimal vector control strategy. Math Comput Simul. 2022;195:31–55. doi: 10.1016/j.matcom.2021.12.021
  • Sangaré B. Contribution à la modélisation mathématique en épidémiologie et en écologie et à la résolution numérique des EDP évolutives. Habilitation à Diriger des Recherches, Université Franche-Comté; 02 Juin 2023.
  • Traoré B, Koutou O, Sangaré B. A global mathematical model of malaria transmission dynamics with structured mosquito population and temperature variations. Nonlinear Anal Real World Appl. 2020;53:Article ID 103081.
  • Luz PM, Codeço CT, Medlock J, et al. Impact of insecticide interventions on the abundance and resistance profile of Aedes aegypti. Epidemiol Infect. 2009;137:1203–1215. doi: 10.1017/S0950268808001799
  • Oki M, Sunahara T, Hashizume M, et al. Optimal timing of insecticide fogging to minimize dengue cases: modeling dengue transmission among various seasonalities and transmission intensities. PLoS Negl Trop Dis. 2011;5:1–7. doi: 10.1371/journal.pntd.0,001,367
  • Alphey N, Alphey L, Bonsall MB. A model framework to estimate impact and cost of genetics-based sterile insect methods for dengue vector control. PLoS One. 2011;6:1–12. doi: 10.1371/journal.pone.0025384
  • Anguelov R, Dumont Y, Lubuma J. Mathematical modeling of sterile insect technology for control of anopheles mosquito. Comput Math Appl. 2012;64:374–389. doi: 10.1016/j.camwa.2012.02.068
  • Dumont Y, Yatat–Djeumen I. Sterile insect technique with accidental releases of sterile females. impact on mosquito-borne diseases control when viruses are circulating. Math Biosci. 2022;343:Article ID 108724. doi: 10.1016/j.mbs.2021.108724
  • Huang M, Song X, Li J. Modelling and analysis of impulsive releases of sterile mosquitoes. J Biol Dyn. 2017;11:147–171. doi: 10.1080/17513758.2016.1254286
  • Roberto CT, Mo YH, Lourdes E. Optimal control of Aedes aegypti mosquitoes by the sterile insect technique and insecticide. Math Biosci. 2010;223(1):12–23. doi: 10.1016/j.mbs.2009.08.009
  • Strugarek M, Bossin H, Dumont Y. On the use of the sterile insect release technique to reduce or eliminate mosquito populations. Appl Math Model. 2018;68:443–470. doi: 10.1016/j.apm.2018.11.026
  • Bliman P-A, Aronna M, Coelho F, et al. Ensuring successful introduction of Wolbachia in natural populations of Aedes aegypti by means of feedback control. J Math Biol. 2018;76:1269–1300. doi: 10.1007/s00285-017-1174-x
  • Himanshu J, Arvind KS. Modeling the effect of Wolbachia to control malaria transmission. Expert Syst Appl. 2023;221:Article ID 119769.
  • Hoffmann A, Montgomery B, Popovici J, et al. Successful establishment of Wolbachia in Aedes populations to suppress dengue transmission. Nature. 2011;476:454–457. doi: 10.1038/nature10356
  • Jair K, Moacyr AHS, Max OS, et al. Aedes, Wolbachia and dengue; 2014.
  • Putri Y, Rozi S, Tasman H, et al. Assessing the effect of extrinsic incubation period (EIP) prolongation in controlling dengue transmission with Wolbachia-infected mosquito intervention.AIP Conference Proceedings, Vol. 1825. Makassar, Indonesia. 2017.
  • Arias J, Martínez H, Sepúlveda-Salcedo L, et al. Predator–prey model for analysis of Aedes aegypti population dynamics in Cali, Colombia. Int J Pure Appl Math. 2015;105(4):561–597. doi: 10.12732/ijpam.v105i4.2
  • Zhongcai Z, Yuanxian H, Linchao H. The impact of predators of mosquito larvae on Wolbachia spreading dynamics. J Biol Dyn. 2023;17(1):Article ID 2249024.
  • Nundloll S, Mailleret L, Grognard F. Two models of interfering predators in impulsive biological control. J Biol Dyn. 2010;4:102–114. doi: 10.1080/17513750902968779
  • Pliego-Pliego E, Vasilieva O, Velázquez-Castro J, et al. Control strategies for a population dynamics model of Aedes aegypti with seasonal variability and their effects on dengue incidence. Appl Math Model. 2020;81:296–319. doi: 10.1016/j.apm.2019.12.025
  • Dumont Y, Chiroleu F, Domerg C. On a temporal model for the Chikungunya disease: modeling, theory and numerics. Math Biosci. 2008;213(1):80–91. doi: 10.1016/j.mbs.2008.02.008
  • Dumont Y, Tchuenche J. Mathematical studies on the sterile insect technique for the Chikungunya disease and Aedes albopictus. J Math Biol. 2012;65:809–854. doi: 10.1007/s00285-011-0477-6
  • Hamdan N, Kilicman A. Mathematical modelling of dengue transmission with intervention strategies using fractional derivatives. Bull Math Biol. 2022;84:138. doi: 10.1007/s11538-022-01096-2
  • Pulecio-Montoya A, López-Montenegro L, Medina-García J. Description and analysis of a mathematical model of population growth of Aedes aegypti. J Appl Math Comput. 2021;65:335–349. doi: 10.1007/s12190-020-01394-9
  • Djeumen V, Couteron P, Jules TJ, et al. An impulsive modelling framework of fire occurrence in a size structured model of tree-grass interactions for savanna ecosystems. J Math Biol. 2017;74:1425–1482. doi: 10.1007/s00285-016-1060-y
  • Miranda GAN, Teboh-Ewungkem I. Infectious diseases and our planet. Cham: Springer; 2021.
  • Ngwa GA, Teboh-Ewungkem MI, Dumont Y, et al. On a three-stage structured model for the dynamics of malaria transmission with human treatment, adult vector demographics and one aquatic stage. J Theor Biol. 2019;481:202–222. doi: 10.1016/j.jtbi.2018.12.043
  • Assefa WW, Rachid O, Jacek B. Mathematical analysis of the impact of transmission-blocking drugs on the population dynamics of malaria. Appl Math Comput. 2021;400:Article ID 126005.
  • Howell P, Knols B. Male mating biology. Malar J. 2009;8(Suppl 2):S8. doi: 10.1186/1475-2875-8-S2-S8
  • Diabaté A, Sangaré B, Koutou O. Mathematical modeling of the dynamics of vector-borne diseases transmitted by mosquitoes: taking into account aquatic stages and gonotrophic cycle. Nonauton Dyn Syst. 2022;9:205–236. doi: 10.1515/msds-2022-0155
  • Diabaté A, Sangaré B, Koutou O. Optimal control analysis of a COVID-19 and tuberculosis (TB) co-infection model with an imperfect vaccine for covid-19. SeMA J. 2023;1–28. doi:10.1007/s40324-023-00330-8.
  • Savadogo A, Sangaré B, Ouédraogo H. A mathematical analysis of Hopf-bifurcation in a prey–predator model with nonlinear functional response. Adv Differ Equ. 2021;2021:275. doi: 10.1186/s13662-021-03437-2
  • Smith HL. Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems. bull new series. Am Math Soc. 1996;33:203–209. doi: 10.1090/S0273-0979-96-00642-8
  • Zhang X, Liu X, Li Y, et al. Modelling the effects of Wolbachia-carrying male augmentation and mating competition on the control of dengue fever. J Dyn Differ Equ. 2023;81:1–41.
  • Koutou O, Diabaté A, Sangaré B. Mathematical analysis of the impact of the media coverage in mitigating the outbreak of covid-19. Math Comput Simul. 2022;205:600–618. doi: 10.1016/j.matcom.2022.10.017
  • Savadogo A, Sangaré B, Ouédraogo H. A mathematical analysis of prey–predator population dynamics in the presence of an SIS infectious disease. RMS Res Math Stat. 2022;9:1–22.
  • Iboi E, Gumel A, Taylor J. Mathematical modeling of the impact of periodic release of sterile male mosquitoes and seasonality on the population abundance of malaria mosquitoes. J Biol Syst. 2020;28:1–34. doi: 10.1142/S0218339020400033
  • Ivanka S, Gani S. Applied impulsive mathematical models. Cham: Springer; 2016.
  • Pang Y, Wang S. Dynamic analysis of the Wolbachia-infected male mosquitoes model by pulsed release. Adv Appl Math. 2021;10:3505–3517. doi: 10.12677/AAM.2021.1010370
  • Pang Y, Wang S, Liu S. Dynamics analysis of stage-structured wild and sterile mosquito interaction impulsive model. J Biol Dyn. 2022;16:464–479. doi: 10.1080/17513758.2022.2079739
  • Drumi DB, Pavel SS. Impulsive differential equations: periodic solutions and applications; 1993.
  • Bainov D, Simeonov PS, Lakshmikantham V. Theory of impulsive differential equations. Plovdiv: Word Scientific Publishing; 1989.
  • Korobeinikov A, Wake G. Lyapunov functions and global stability for SIR, SIRS, and SIS epidemiological models. Appl Math Lett. 2002;15(8):955–960. doi: 10.1016/S0893-9659(02)00069-1
  • Manyombe MLM, Mbang J, Tsanou B, et al. Mathematical analysis of a spatiotemporal model for the population ecology of anopheles mosquito. Math Methods Appl Sci. 2020;43:3524–3555. doi: 10.1002/mma.v43.6
  • Hughes H, Britton N. Modeling the use of Wolbachia to control dengue fever transmission. Bull Math Biol. 2013;75:796–818. doi: 10.1007/s11538-013-9835-4
  • Dumont Y, Chiroleu F. Vector control for the Chikungunya disease. Math Biosci Eng. 2010;7:313–345. doi: 10.3934/mbe.2010.7.313