Controlling Wolbachia Transmission and Invasion Dynamics among Aedes Aegypti Population via Impulsive Control Strategy

Abstract

This work is devoted to analyzing an impulsive control synthesis to maintain the self-sustainability of Wolbachia among Aedes Aegypti mosquitoes. The present paper provides a fractional order Wolbachia invasive model. Through fixed point theory, this work derives the existence and uniqueness results for the proposed model. Also, we performed a global Mittag-Leffler stability analysis via Linear Matrix Inequality theory and Lyapunov theory. As a result of this controller synthesis, the sustainability of Wolbachia is preserved and non-Wolbachia mosquitoes are eradicated. Finally, a numerical simulation is established for the published data to analyze the nature of the proposed Wolbachia invasive model.

1. Introduction

In the 19th century, fractional calculus (FC) theory has been built by some famous mathematicians like Grunwald, Letnikov, Riemann, Liouville, Euler and Caputo [1,2,3]. Fractional order derivatives are the generalization of integer order derivatives. FC is unavoidable due to its extensive applications in the study of real-world problems. The main advantage of FC is that it can provide a path to understand the description of memory and inheritance of various processes [4,5]. The book [6] plays an important role in the area of applied fractional calculus. In recent years, researchers in the field of physics, chemistry, Neural Networks, economic and mathematical modeling, biological problems and engineering have been very much attracted to fractional calculus [7], because FC interprets the whole function geometrically and globalizes its entire function.

Mosquito-borne diseases are primarily spread by female mosquitoes while taking a blood meal from living organisms such as humans, animals and birds. A parasite, virus, or bacteria-infected female mosquito can transmit those foreign agents to humans [8]. For instance, the Dengue virus, Zika virus, Yellow fever virus and Chikungunya are transmitted from infected human to uninfected human via primary vector Aedes Aegypti mosquitoes. Currently, the secondary vector for the above-mentioned diseases is Aedes Albopictus [9,10,11]. In recent years, the death rate due to mosquito-borne diseases has increased dramatically [8]. Gubler et al. [12,13] and Ong et al. [14] explained that dengue and dengue hemorrhagic fever are a more common issue for public health. According to the World Health Organization (WHO) [15], per annum, mosquito-borne diseases cause more than 40,000 deaths and 96 million asymptomatic cases in 129 countries.

Currently, there are several methods to control Aedes Aegypti mosquitoes such as insecticide spraying, sterile insect technique, incompatible insect technique, combined sterile insect technique, and genetic modifications. In [16,17], the authors proposed that the Sterile insect technique is likely to be used in mosquito-borne disease control. The authors of [18], analyzed that the particular transgenic strain can simulate the female-specific flightless phenotype to increase the sterilization in male mosquitoes. In [19,20], the authors discussed that the safe and effective replacement of vector population by genetically modified mosquitoes will play a significant role in mosquito-borne disease control. Furthermore, some other types of mosquito control strategies, such as making changes in feeding behaviors, intervention strategies, using bed nets and mosquito repellents, are also tested [21,22].

A novel Aedes Aegypti suppression technique using the life-shortening bacterium Wolbachia plays an important role [23,24,25]. It is an endosymbiotic bacterium that is reported in nearly 60 percent of insect species by Wolbach (1924) [26]. The World Mosquito Program (WMP) [27] from Australia currently release Wolbachia infected mosquitoes over 10 countries, such as countries in Latin America, India, Sri Lanka, Vietnam, Indonesia and cities in Oceania. In that research, they found that Wolbachia is a self-sustaining bacterium and in the presence of Wolbachia infected mosquitoes there is zero possibility of having Dengue. The Wolbachia releasing strategy is more powerful than that of the above-mentioned control strategies in the sense that it is self-sustaining, affordable, only needs a small amount of release, the area covered is larger than the released area, and the most important thing is it is not harmful to human health. The authors of [28,29,30,31] discussed that Wolbachia can restrict the virus particles of various diseases. We know that the virus is transmitted from infected humans to uninfected humans via female mosquitoes. Meanwhile, if a virus-infected mosquito carries Wolbachia strain, then the virus cannot be transmitted to an uninfected human. Because this Wolbachia strain blocks the virus particles inside the salivary gland of mosquitoes (Ref. Figure 1).

The Wolbachia infection is introduced into wild mosquitoes population through two major processes such as microinjection and Introgression [32].

• Micro injection: In this process, Wolbachia strains are microinjected into aquatic stages such as eggs, larvae and pupae.
• Introgression: In this process, the Wolbachia strains are carried out to next generation through mating. If Wolbachia infected female mated with Wolbachia infected or uninfected male, then the produced offsprings have the Wolbachia strain (Called CI rescue). Suppose the Wolbachia uninfected female mated with a Wolbachia infected male then there is no viable progeny. Finally, if a non-Wolbachia female mated with a non-Wolbachia male then there is no Wolbachia infection in the offspring.

To understand the introgression process, one can refer to Figure 2.

Furthermore, some existing mathematical models consider Wolbachia as a control agent for mosquito-borne diseases. In [33], the author proposed a deterministic model to control mosquito-borne diseases up to 90% via Wolbachia spread, also the author considered both human and mosquito populations to create a mathematical model. In [30,34], the authors proposed a mathematical model depicting the life stages of mosquitoes with Wolbachia and proved that Wolbachia has an excellent quality to control dengue virus spread. In [35], the authors analyzed the integer ordered mathematical model consisting of only four stages (aquatic stage with and without Wolbachia and adult female mosquitoes with and without Wolbachia), which considered the imperfect maternal transmission and Wolbachia invasion. In [36], the two sex mathematical model is discussed to analyze the persistence of Wolbachia. In [37], the age and bite structured mathematical model is proposed and performed the mathematical analysis. In [38], the authors discussed the linear feedback control strategy of a mathematical model containing only three stages such as aquatic, female Wolbachia infected and uninfected mosquitoes. In this, the author analyzed the Wolbachia infected mosquitoes release into the seasonal environment. In [39], the authors presented a mathematical model to depict the mechanism of the virus inside both humans and mosquitoes. In this work, the author utilized two various types of controls like vaccination for humans and Wolbachia infected mosquitoes’ release for mosquitoes. The pontriyagin maximum principle was utilized to analyze the optimal control of the proposed mathematical model. In [40], the authors discussed the Wolbachia infection among Aedes Aegypti mosquitoes via delay differential equations. In that work, the author proposed the delay dependent stability criteria of the proposed model by utilizing the results from spectrum analysis. In [41], the authors proposed an age structured fractional order mathematical model to control the Aedes Aegypti mosquitoes via Wolbachia bacterium using the Linear Matrix Inequality (LMI) approach.

As per the practical results of [27], Wolbachia should be released into every stage to get the optimal control in a short period. Also, by utilizing fractional calculus we can get the memory property and inheritance of this process. In nature, Wolbachia infected mosquitoes may lose the Wolbachia infection. Because of this, invasion in Wolbachia is unavoidable. Motivation by the above discussions, our contributions are listed below:

• A novel mathematical model, which considers the total of ten stages in Aedes Aegypti mosquitoes (combining both Wolbachia infected and Wolbachia uninfected) is proposed and the possible optimal stages to release the Wolbachia are discussed, and the most important concept of Wolbachia invasion and Wolbachia gain are adopted.
• The Wolbachia free equilibrium, Wolbachia present Equilibrium, Zero mosquitoes, and both Wolbachia and Non-Wolbachia mosquitoes co-existence equilibrium are derived. And utilizing fixed point theory results, the Existence and Uniqueness results of the Wolbachia invasive model are proposed. To attain optimal control, we utilized an impulsive control strategy.
• We perform global Mittag-Leffler stability analysis of the proposed model via Linear Matrix Inequality (LMI) theory and Lyapunov theory.
• In the end, by utilizing the data from the published literature, we have presented the numerical simulation of the proposed model using MATLAB software.

The rest of the paper is arranged as follows—in Section 2, we provide some basic Definitions, Lemmas and Theorems. In Section 3, the fractional order complete mathematical model describes the interaction between Wolbachia infected and Non-Wolbachia mosquitoes is presented. In Section 4, the possible equilibrium points are presented. In Section 5, the Wolbachia invasive and gain model with impulsive control is presented. In Section 6, the existence and uniqueness results are analyzed and the global Mittag-Leffler stability results are derived in Section 7. In Section 8, the numerical simulation results are presented. In Section 9, the work is concluded.

Notations.$\mathbb{N}$ denotes the space of all natural numbers, $\mathbb{R}$ denotes the space of all real numbers, $\mathbb{C}$ denotes the space of all complex numbers, ${\mathbb{R}}^n$ denotes the space of n-dimensional Euclidean space, ${\mathbb{Z}}^{+}$ denotes the space of all positive integers. Moreover, $Re(·)$ denotes the real part of a complex number and $[.]$ denotes the integer part of a number. * denotes the corresponding symmetric terms in a symmetric matrix. Also, ${}_k^c{D}_t^{\alpha }(·)$ and ${}_k^c{I}_t^{\alpha }(·)$ denotes the derivative and anti derivative of order $\alpha$ with respect to t respectively, c denotes that its in Caputo sense, k denotes the initial condition and $\Gamma (·)$ denotes the Gamma function.

2. Preliminaries

In this section, we provide some basic Definitions, Lemmas and Theorems, which are used to attain our results.

Definition 1.

Ref. [4] The most important basic function in fractional calculus is the gamma function. It is defined as follows:

$$\begin{array}{ccc}\hfill \Gamma \left(z\right)& =& {\int }_0^{\infty }{e}^{-s}{s}^{z-1}ds,\hfill \end{array}$$

with $Re\left(z\right)>0$.

Definition 2.

Ref. [1] The Caputo fractional derivative of a continuous function $f\left(t\right)$ over $[k,T]$ of order $\alpha \in \mathbb{C}$ (with $Re\left(\alpha \right)>0$, $\alpha \notin \mathbb{N})$ is

$$\begin{array}{ccc}\hfill {}_k^c{D}_t^{\alpha }f\left(t\right)& =& \frac{1}{\Gamma (n-\alpha )}\left[{\int }_k^t{(t-\eta )}^{n-\alpha -1}\frac{d^n}{d{\eta }^n}f\left(\eta \right)d\eta \right],\hfill \end{array}$$

(1)

where, $n=\left[Re\right(\alpha \left)\right]+1$.

If $0<Re\left(\alpha \right)<1$, then the expression (1), can be rewritten as

$$\begin{array}{ccc}\hfill {}_k^c{D}_t^{\alpha }f\left(t\right)& =& \frac{1}{\Gamma (1-\alpha )}\left[{\int }_k^t\frac{f^{{}^{\prime }}\left(\eta \right)d\eta }{{(t-\eta )}^{\alpha }}\right].\hfill \end{array}$$

(2)

Since, $n=1$ for all $0<Re\left(\alpha \right)<1$.

Definition 3.

Ref. [42] The Caputo sense fractional integral of a continuous function f on $L^1([0,T],\mathbb{R})$ over $\alpha \in (0,1]$ with respect to t is defined as

$$\begin{array}{ccc}\hfill {}_0^c{I}_t^{\alpha }f\left(t\right)& =& \frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}f\left(\eta \right)d\eta .\hfill \end{array}$$

(3)

The two parameter Mittag-Leffler function is defined as follows: [4]

$$\begin{array}{ccc}\hfill E_{a,b}\left(z\right)& =& \sum _{l=0}^{\infty }\frac{z^l}{\Gamma (al+b)},\hfill \end{array}$$

where, $z\in \mathbb{C}$, $a>0$, and $b>0$. If $b=1$ the $E_a\left(z\right)={\sum }_{l=0}^{\infty }\frac{z^l}{\Gamma (al+1)}$. If both $a=1$ and $b=1$, the $E_{1,1}\left(z\right)=e^z$.

Lemma 1

(Schur Complement [43]). Let us denote three $n\times n$ matrices as ${\mathsf{\Psi }}_1,{\mathsf{\Psi }}_2,{\mathsf{\Psi }}_3$, where ${\mathsf{\Psi }}_1={\mathsf{\Psi }}_1^{\top }$ and ${\mathsf{\Psi }}_2={\mathsf{\Psi }}_2^{\top }>0$. Then ${\mathsf{\Psi }}_1+{\mathsf{\Psi }}_3^{\top }{\mathsf{\Psi }}_2^{-1}{\mathsf{\Psi }}_3<0$ if and only if $\left[\begin{array}{cc}{\mathsf{\Psi }}_1& {\mathsf{\Psi }}_3^{\top }\\ {\mathsf{\Psi }}_3& -{\mathsf{\Psi }}_2\end{array}\right]<0$ (or) $\left[\begin{array}{cc}-{\mathsf{\Psi }}_2& {\mathsf{\Psi }}_3\\ {\mathsf{\Psi }}_3^{\top }& {\mathsf{\Psi }}_1\end{array}\right]<0$.

Lemma 2.

Ref. [44] For any scalar $ϵ>0$, $A,N\in {\mathbb{R}}^n$ and matrix $P_1$, then

$$\begin{array}{ccc}\hfill A^{\top }{P}_{1}N& \le & \frac{1}{2ϵ}{A}^{\top }{P}_1{P}_1^{\top }A+\frac{ϵ}{2}{N}^{\top }N.\hfill \end{array}$$

Let us consider the fractional order dynamical system with impulse of type,

$$\begin{array}{ccc}\hfill {}_k^c{D}_t^{\alpha }x\left(t\right)& =& -A_{1}x\left(t\right)+A_{2}f\left(x\left(t\right)\right),t\ne t_{\theta },\theta =1,2,\cdots ,m,\hfill \\ \hfill \Delta x\left(t_{\theta }\right)& =& x\left(t_{\theta }^{+}\right)-x\left(t_{\theta }^{-}\right)={\delta }_{\theta }\left(x\left(t_{\theta }\right)\right),t=t_{\theta },\theta =1,2,\cdots ,m,\hfill \end{array}$$

(4)

with initial condition $x\left(t_0\right)=x_0\in {\mathbb{Z}}^{+}$, where the n states is defined by $x\left(t\right)={[x_1\left(t\right),x_2\left(t\right),x_3\left(t\right),\cdots ,x_n\left(t\right)]}^{\top }\in {\mathbb{R}}^n$ and $f\left(x\left(t\right)\right)={[f\left(x_1\left(t\right)\right),f\left(x_2\left(t\right)\right),f\left(x_3\left(t\right)\right),\cdots ,f\left(x_n\left(t\right)\right)]}^{\top }$ be a function, $A_1$ and $A_2$ are constant coefficient matrices with the impulsive operator ${\delta }_{\theta }:{\mathbb{R}}^n\to {\mathbb{R}}^n$.

Definition 4.

Ref. [44] The system (4), is said to be globally Mittag-Leffler stable at its equilibrium points, if the following hold:

$$\begin{array}{ccc}\hfill \left|\right|x\left(t\right)-x^{*}\left|\right|& \le & {\left[h(x_0-x^{*})E_{\alpha }(-\kappa t^a)\right]}^b,\hfill \end{array}$$

where $x^{*}$ is an equilibrium point, $0<\alpha <1$, $\kappa \ge 0$ and $a,b>0$. Moreover, $h\left(0\right)=0$, $h\left(x\right)\ge 0$ and $h\left(x\right)$ is locally Lipschitz with Lipchitz constant $h_0$.

Lemma 3.

Ref. [45] Let us consider the fractional order system with impulsive control of type (4). Suppose $f\left(0\right)=0$, $t>0$ and ${\delta }_{\theta }\left(0\right)=0$, $\theta =1,2,3,\cdots ,m$. If there exists a positive definite function V such that the following hold:

(1) 

There exists positive constants ${\alpha }_1$ and ${\alpha }_2$

$$\begin{array}{ccc}\hfill {\alpha }_1\left|\right|x\left(t\right)\left|\right|& \le & V\left(t\right)\le {\alpha }_2\left|\right|x\left(t\right)\left|\right|,x\left(t\right)\in {\mathbb{R}}^n.\hfill \end{array}$$

(2) 

${}_0^c{D}_t^{\alpha }V\left(t\right)\le -{ϵ}_{1}V\left(t\right)$, $t\ne t_{\theta }$, $\theta =1,2,3,\cdots ,m$ for any scalar ${ϵ}_1$.

(3) 

$V\left(t_{\theta }^{+}\right)\le V\left(t_{\theta }\right)$, $t=t_{\theta }$, $\theta =1,2,3,\cdots ,m$.

then the equilibrium point of the system (4) is globally Mittag-Leffler stable.

Definition 5.

Ref. [46] A map $\nu :H\to H$, H compact Banach space, is said to be a contraction mapping if there exists $h\in (0,1)$ such that

$$\begin{array}{ccc}\hfill \left|\right|\nu \left(m_1\right)-\nu \left(m_2\right)\left|\right|& \le & h\left|\right|m_1-m_2\left|\right|\hfill \end{array}$$

for every $m_1,m_2\in H$.

Theorem 1

(Contraction Mapping Theorem). Ref. [46] Suppose H is a complete metric space and $\nu :H\to H$ is a contraction mapping. Then, ν has a unique fixed point.

3. Model Formulation

In this section, a novel mathematical model is proposed to expose the transmission dynamics of the gram negative bacteria Wolbachia among Aedes Aegypti mosquitoes. While constructing the model we have considered the total of 10 stages such as non-Wolbachia eggs $\left(W_e\right)$, non-Wolbachia larvae $\left(W_l\right)$, non-Wolbachia pupae $\left(W_p\right)$, non-Wolbachia adult female $\left(W_f\right)$, non-Wolbachia adult male $\left(W_a\right)$, Wolbachia infected eggs $\left(I_e\right)$, Wolbachia infected larvae $\left(I_l\right)$, Wolbachia infected pupae $\left(I_p\right)$, Wolbachia infected adult female $\left(I_f\right)$, Wolbachia infected adult male $\left(I_a\right)$. The total population at time t is denoted as $T=W_e\left(t\right)+W_l\left(t\right)+W_p\left(t\right)+W_f\left(t\right)+W_a\left(t\right)+I_e\left(t\right)+I_l\left(t\right)+I_p\left(t\right)+I_f\left(t\right)+I_a\left(t\right)$. The eggs with zero Wolbachia infection are produced at the rate ${\Lambda }_{w_e}$ by the mating process between non-Wolbachia female $\left(W_f\right)$ and non-Wolbachia male $\left(W_a\right)$. There is no other possibilities of having a non-Wolbachia eggs. Therefore, the reproduction rate of non-Wolbachia mosquitoes can be calculated by the term $\frac{{\Lambda }_{w_e}{W}_f{W}_a}{T}$. Along with this, the terms ${\lambda }_{w_e}$ (natural mortality rate of non-Wolbachia eggs) and ${\gamma }_{w_e}$ (maturation rate of non-Wolbachia eggs) denotes the limitations in the growth of wild mosquito eggs. At the same time, after release of Wolbachia infected mosquitoes (in both aquatic and ariel stage) in a common environment, the production of Wolbachia infected mosquito eggs $I_e\left(t\right)$, depends on mating between Wolbachia infected female $I_f\left(t\right)$ and non-Wolbachia male $W_a\left(t\right)$ and from mating between Wolbachia infected female $I_f\left(t\right)$ and Wolbachia infected male $I_a\left(t\right)$. Through this, the birth rate of Wolbachia infected mosquito eggs population $I_e\left(t\right)$ with the reproduction rate ${\Lambda }_{i_e}$ is

$$\begin{array}{ccc}\hfill \frac{{\Lambda }_{i_e}(I_f{W}_a+I_f{I}_a)}{T}& =& \frac{{\Lambda }_{i_e}{I}_f(W_a+I_a)}{T}.\hfill \end{array}$$

Similarly, the increase in the growth of Wolbachia infected eggs is limited by the natural mortality rate ${\lambda }_{i_e}$ and the maturation rate ${\gamma }_{i_e}$ (That is, the rate in which the corresponding compartment moved into the next stage).

Furthermore, the quantity $(1-\alpha ){\gamma }_{i_e}{I}_e$ is added to the wild mosquito larvae population. Because the term $\alpha$ and $(1-\alpha )$ denotes the probability of getting larvae with and without Wolbachia respectively. Similarly, $\beta$ and $(1-\beta )$ denotes the probability of getting pupae with and without Wolbachia respectively, $ϵ$ and $(1-ϵ)$ denotes the probability rate of having Wolbachia infection in adult mosquitoes by introgression. That is, $ϵ$ be the probability of getting Wolbachia infected adults (with ${\rho }_{i_w}$ = probability of getting male and $(1-{\rho }_{i_w})$ = probability of getting female). Because of these reasons, the terms $(1-\alpha ){\gamma }_{i_e}{I}_e$, $(1-\beta ){\gamma }_{i_l}{I}_l$, $(1-ϵ){\gamma }_{i_p}{\rho }_{i_w}{I}_p$ and $(1-ϵ){\gamma }_{i_p}(1-{\rho }_{i_w})I_p$ are added to the corresponding stages and similarly, the terms $\alpha {\gamma }_{i_e}{I}_e$, $\beta {\gamma }_{i_l}{I}_l$ and $ϵ{\gamma }_{i_p}{I}_{i_p}$ are removed from the corresponding stages. The parameter description of the system of Equation (5) is presented in

Table 1. Description of parameters involved in system of Equation (5).

Parameter	Description
${\Lambda }_{w_e}$, ${\Lambda }_{i_e}$	Reproduction rate of non-Wolbachia mosquitoes
	and Wolbachia infected mosquitoes respectively
${\lambda }_{w_e}$	The natural death rate of eggs without Wolbachia infection
${\lambda }_{w_l}$	The natural death of larvae without Wolbachia infection
${\lambda }_{w_p}$	The natural death of pupae without Wolbachia infection
${\lambda }_{w_f}$	The natural death of adult female mosquitoes without Wolbachia infection
${\lambda }_{w_a}$	The natural death of adult male mosquitoes without Wolbachia infection
${\lambda }_{i_e}$	The natural death of eggs with Wolbachia infection
${\lambda }_{i_l}$	The natural death of larvae with Wolbachia infection
${\lambda }_{i_p}$	The natural death of pupae with Wolbachia infection
${\lambda }_{i_f}$	The natural death of adult female mosquitoes with Wolbachia infection
${\lambda }_{i_a}$	The natural death of infected adult male mosquitoes with Wolbachia infection
${\gamma }_{w_e}$	The rate at which the fraction of non-Wolbachia eggs matured into non-Wolbachia larvae
${\gamma }_{w_l}$	The rate at which the fraction of non-Wolbachia larvae matured into non-Wolbachia pupae
${\gamma }_{w_p}$	The rate at which the fraction of non-Wolbachia pupae matured into non-Wolbachia
	immature female or male
${\gamma }_{i_e}$	The rate at which the fraction of the Wolbachia infected mosquito eggs
	matured into Wolbachia infected or uninfected larvae
${\gamma }_{i_l}$	The rate at which the fraction of the Wolbachia infected mosquito larvae
	matured into Wolbachia infected or uninfected pupae
${\gamma }_{i_p}$	The rate at which the fraction of the Wolbachia infected mosquito pupae
	matured into Wolbachia infected or uninfected adults
$\rho$	The probability of having male or female mosquitoes

.

From the above facts, the novel mathematical model that describes the transmission dynamics of Wolbachia among Aedes Aegypti mosquitoes is proposed as follows:

$$\left\{\begin{array}{cc}\hfill {}_0^c{D}_t^{\alpha }{W}_e& =\frac{{\Lambda }_{w_e}{W}_f{W}_a}{T}-{\lambda }_{w_e}{W}_e-{\gamma }_{w_e}{W}_e\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_l& ={\gamma }_{w_e}{W}_e-{\lambda }_{w_l}{W}_l-{\gamma }_{w_l}{W}_l+(1-\alpha ){\gamma }_{i_e}{I}_e\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_p& ={\gamma }_{w_l}{W}_l-{\lambda }_{w_p}{W}_p-{\gamma }_{w_p}{W}_p+(1-\beta ){\gamma }_{i_l}{I}_l\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_f& =\rho {\gamma }_{w_p}{W}_p-{\lambda }_{w_f}{W}_f+(1-ϵ){\gamma }_{i_p}{\rho }_{i_w}{I}_p\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_a& =(1-\rho ){\gamma }_{w_p}{W}_p-{\lambda }_{w_a}{W}_a+(1-ϵ){\gamma }_{i_p}(1-{\rho }_{i_w})I_p\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_e& =\frac{{\Lambda }_{i_e}{I}_f(W_a+I_a)}{T}-{\lambda }_{i_e}{I}_e-\alpha {\gamma }_{i_e}{I}_e\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_l& =\alpha {\gamma }_{i_e}{I}_e-{\lambda }_{i_l}{I}_l-\beta {\gamma }_{i_l}{I}_l\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_p& =\beta {\gamma }_{i_l}{I}_l-{\lambda }_{i_p}{I}_p-ϵ{\gamma }_{i_p}{I}_p\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_f& ={\rho }_iϵ{\gamma }_{i_p}{I}_p-{\lambda }_{i_f}{I}_f\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_a& =(1-{\rho }_i)ϵ{\gamma }_{i_p}{I}_p-{\lambda }_{i_a}{I}_a.\hfill \end{array}\right.$$

(5)

The dynamics of the population can be easily understand by the schematic diagram Figure 3 and the parameters are described in Table 1.

4. Equilibrium Points

In this section, we can find the four cases of possible equilibrium points such as wild mosquitoes only, Wolbachia mosquitoes only, co-existence of both population and zero mosquitoes.

4.1. Zero Mosquitoes

Suppose there is no mosquitoes, then the equilibrium point can be written as $P_1=(0,0,0,0,0,0,0,0,0,0)$. This is trivial but does not exists in nature.

4.2. Wolbachia Infected Mosquitoes Free Equilibrium

Suppose, there is no Wolbachia infected mosquitoes population then the possible equilibrium can be written as

$$\begin{array}{c}\hfill P_2=(W_{e_1}^{*},W_{l_1}^{*},W_{p_1}^{*},W_{f_1}^{*},W_{a_1}^{*},0,0,0,0,0),\end{array}$$

where,

$$\begin{array}{ccc}\hfill W_{e_1}^{*}& =& \frac{T{\lambda }_{w_f}{\lambda }_{w_a}({\lambda }_{w_e}+{\gamma }_{w_e}){({\lambda }_{w_l}+{\gamma }_{w_l})}^2{({\lambda }_{w_p}+{\gamma }_{w_p})}^2}{\rho (1-\rho ){\Lambda }_{w_e}{\gamma }_{w_p}^2{\gamma }_{w_e}^2{\gamma }_{w_l}^2}\hfill \\ \hfill W_{l_1}^{*}& =& \frac{{\gamma }_{w_e}}{{\lambda }_{w_l}+{\gamma }_{w_l}}{W}_{e_1}^{*}\hfill \\ \hfill W_{p_1}^{*}& =& \frac{{\gamma }_{w_l}{\gamma }_{w_e}}{({\lambda }_{w_l}+{\gamma }_{w_l})({\lambda }_{w_p}+{\gamma }_{w_p})}{W}_{e_1}^{*}\hfill \\ \hfill W_{f_1}^{*}& =& \frac{\rho {\gamma }_{w_p}{\gamma }_{w_e}{\gamma }_{w_l}}{{\lambda }_{w_f}({\lambda }_{w_f}+{\gamma }_{w_f})({\lambda }_{w_p}+{\gamma }_{w_p})}{W}_{e_1}^{*}\hfill \\ \hfill W_{a_1}^{*}& =& \frac{(1-\rho ){\gamma }_{w_p}{\gamma }_{w_e}}{{\lambda }_{w_e}({\lambda }_{w_l}+{\gamma }_{w_l})({\lambda }_{w_p}+{\gamma }_{w_p})}{W}_{e_1}^{*}\hfill \end{array}$$

4.3. Wild Mosquitoes Free Equilibrium

After the successful replacement of Wolbachia uninfected mosquitoes by Wolbachia infected mosquitoes the equilibrium point can be represented by

$$\begin{array}{c}\hfill P_3=(0,0,0,0,0,I_{e_2}^{*},I_{l_2}^{*},I_{p_2}^{*},I_{f_2}^{*},I_{a_2}^{*}),\end{array}$$

where,

$$\begin{array}{ccc}\hfill I_{e_2}^{*}& =& \frac{({\lambda }_{i_l}+\beta {\gamma }_{i_l})({\lambda }_{i_p}+ϵ{\gamma }_{i_p})}{\alpha \beta {\gamma }_{i_e}{\gamma }_{i_l}}{I}_{p_2}^{*}\hfill \\ \hfill I_{l_2}^{*}& =& \frac{({\lambda }_{i_p}+ϵ{\gamma }_{i_p})}{\beta {\gamma }_{i_l}}{I}_{p_2}^{*}\hfill \\ \hfill I_{p_2}^{*}& =& \frac{T{\lambda }_{i_f}{\lambda }_{i_a}({\lambda }_{i_e}+\alpha {\gamma }_{i_e})({\lambda }_{i_l}+\beta {\gamma }_{i_l})({\lambda }_{i_p}+ϵ{\gamma }_{i_p})}{{\Lambda }_{i_e}\alpha \beta {\rho }_i(1-{\rho }_i){ϵ}^2{\gamma }_{i_p}^2{\gamma }_{i_e}{\gamma }_{i_l}}\hfill \\ \hfill I_{f_2}^{*}& =& \frac{{\rho }_iϵ{\gamma }_{i_p}}{{\lambda }_{i_f}}{I}_{p_2}^{*}\hfill \\ \hfill I_{a_2}^{*}& =& \frac{(1-{\rho }_i)ϵ{\gamma }_{i_p}}{{\lambda }_{i_a}}{I}_{p_2}^{*}.\hfill \end{array}$$

4.4. Both Wolbachia Infected Mosquitoes and Non-Wolbachia Mosquitoes Co-Existence Equilibrium

If both Wolbachia infected and Wolbachia uninfected mosquitoes present in common environment, then the equilibrium point is

$$\begin{array}{ccc}\hfill S_n& =& \left\{W_{e_n}^{*},W_{l_n}^{*},W_{p_n}^{*},W_{f_n}^{*},W_{a_n}^{*},I_{e_n}^{*},I_{l_n}^{*},I_{p_n}^{*},I_{f_n}^{*},I_{a_n}^{*}\right\},n=3,4.\hfill \end{array}$$

$$\begin{array}{ccc}\hfill W_{e_n}^{*}& =& \left(\frac{{\lambda }_{w_l}+{\gamma }_{w_l}}{{\gamma }_{w_e}}\right)\left(\frac{{\lambda }_{w_p}+{\gamma }_{w_p}}{{\gamma }_{w_l}}\right)\left(\frac{T{B}_1{B}_2{B}_3{\lambda }_{i_f}{\lambda }_{w_a}}{{\Lambda }_{i_e}(1-\rho ){\rho }_i{\gamma }_{w_p}}\right)-\frac{I_{a_n}^{*}}{{\gamma }_{w_e}}\hfill \\ & & \left[B_4({\lambda }_{w_l}+{\gamma }_{w_l})\left(\frac{{\lambda }_{w_p}+{\gamma }_{w_p}}{{\gamma }_{w_l}}\right)+({\lambda }_{w_l}+{\gamma }_{w_l})\left(\frac{(1-\beta ){\gamma }_{i_l}}{{\gamma }_{w_l}}\right)\left(\frac{{\lambda }_{i_a}{B}_1}{\beta {\gamma }_{i_l}(1-{\rho }_i)}\right)+\frac{(1-\alpha ){\gamma }_{i_e}{\lambda }_{i_a}{B}_1{B}_2}{\alpha {\gamma }_{i_e}(1-{\rho }_i)}\right]\hfill \end{array}$$

$$\begin{array}{ccc}\hfill W_{l_n}^{*}& =& \left(\frac{{\lambda }_{w_p}+{\gamma }_{w_p}}{{\gamma }_{w_l}}\right)\left[\frac{T{B}_1{B}_2{B}_3{\lambda }_{i_f}{\lambda }_{w_a}}{{\Lambda }_{i_e}(1-\rho ){\rho }_i{\gamma }_{w_p}}-B_4{I}_{a_n}^{*}\right]-\left(\frac{(1-\beta ){\gamma }_{i_l}}{{\gamma }_{w_l}}\right)\left[\frac{{\lambda }_{i_a}{B}_1}{\beta {\gamma }_{i_l}(1-{\rho }_i)}\right]\hfill \end{array}$$

$$\begin{array}{ccc}\hfill W_{p_n}^{*}& =& \left[\frac{T{B}_1{B}_2{B}_3{\lambda }_{i_f}{\lambda }_{w_a}}{{\Lambda }_{i_e}(1-\rho ){\rho }_i{\gamma }_{w_p}}-B_4{I}_{a_n}^{*}\right]\hfill \end{array}$$

$$\begin{array}{ccc}\hfill W_{f_n}^{*}& =& \frac{\rho {\gamma }_{w_p}}{{\lambda }_{w_f}}\left[\frac{T{B}_1{B}_2{B}_3{\lambda }_{i_f}{\lambda }_{w_a}}{{\Lambda }_{i_e}(1-\rho ){\rho }_i{\gamma }_{w_p}}-B_4{I}_{a_n}^{*}\right]+\frac{(1-ϵ){\rho }_{i_w}{\lambda }_{i_a}}{ϵ{\lambda }_{w_f}(1-{\rho }_i)}{I}_{a_n}^{*}\hfill \end{array}$$

$$\begin{array}{ccc}\hfill W_{a_n}^{*}& =& \frac{T{B}_1{B}_2{B}_3{\lambda }_{i_f}}{{\rho }_i{\Lambda }_{i_e}}-I_{a_n}^{*}\hfill \end{array}$$

$$\begin{array}{ccc}\hfill I_{e_n}^{*}& =& \frac{B_1{B}_2{\lambda }_{i_a}{I}_{a_n}^{*}}{\alpha {\gamma }_{i_e}(1-{\rho }_i)}\hfill \end{array}$$

$$\begin{array}{ccc}\hfill I_{p_n}^{*}& =& \frac{{\lambda }_{i_a}}{(1-{\rho }_i)ϵ{\gamma }_{i_p}}{I}_{a_n}^{*}\hfill \end{array}$$

$$\begin{array}{ccc}\hfill I_{f_n}^{*}& =& \frac{{\rho }_i{\lambda }_{i_a}{I}_{a_n}^{*}}{(1-{\rho }_i){\lambda }_{i_f}}\hfill \end{array}$$

with $I_{a_3}^{*}>I_{a_4}^{*}$, both roots can be found from the quadratic equation

$$\begin{array}{ccc}\hfill a_1{I}_a^{{*}^2}+a_2{I}_a^{*}+a_3& =& 0,\hfill \end{array}$$

where,

$$\begin{array}{ccc}\hfill a_1& =& \frac{{\Lambda }_{w_e}\rho B_4{\gamma }_{w_p}}{T{\lambda }_{w_f}};\hfill \\ \hfill a_2& =& \left(\frac{{\lambda }_{w_e}+{\gamma }_{w_e}}{T{\lambda }_{w_f}}\right)\left(\frac{{\lambda }_{w_e}{\lambda }_{i_f}\rho B_1{B}_2{B}_3}{{\rho }_i{\Lambda }_{i_e}{\lambda }_{w_f}}\right)\left(\frac{{\lambda }_{w_a}}{(1-\rho )}+B_4{\gamma }_{w_p}\right)\hfill \\ & & \left(\frac{({\lambda }_{w_l}+{\gamma }_{w_l})({\lambda }_{w_p}+{\gamma }_{w_p})B_4}{{\gamma }_{w_l}}+\frac{({\lambda }_{w_l}+{\gamma }_{w_l})(1-\beta ){\lambda }_{i_a}{B}_1}{{\gamma }_{w_l}\beta (1-{\rho }_i)}+\frac{(1-\alpha ){\lambda }_{i_a}{B}_1{B}_2}{\alpha (1-{\rho }_i)}\right);\hfill \\ \hfill a_3& =& \frac{{\Lambda }_{w_e}\rho T{B}_1^2{B}_2^2{B}_3^2{\lambda }_{w_a}}{{\Lambda }_{i_e}^2(1-\rho ){\rho }_i^2{\lambda }_{w_f}}.\hfill \end{array}$$

Here,

$$\begin{array}{ccc}\hfill B_1& =& 1+\frac{{\lambda }_{i_p}}{ϵ{\gamma }_{i_p}};\hfill \\ \hfill B_2& =& 1+\frac{{\lambda }_{i_l}}{\beta {\gamma }_{i_l}};\hfill \\ \hfill B_3& =& 1+\frac{{\lambda }_{i_e}}{\alpha {\gamma }_{i_e}};\hfill \\ \hfill B_4& =& 1+\frac{(1-ϵ)(1-{\rho }_{i_w}){\lambda }_{i_a}}{(1-\rho )(1-{\rho }_i)ϵ{\gamma }_{w_p}}.\hfill \end{array}$$

For more details about the calculations of Section 4, kindly refer the Appendix A section.

5. Wolbachia Invasion Model

We considered the possibility of Wolbachia loss in adult mosquitoes and possibility of Wolbachia gain in aquatic stage mosquitoes. Then Equation (5), can be rewritten as

$$\left\{\begin{array}{cc}\hfill {}_0^c{D}_t^{\alpha }{W}_e\left(t\right)& =\frac{{\Lambda }_{w_e}{W}_f{W}_a}{T}-{\lambda }_{w_e}{W}_e-{\gamma }_{w_e}{W}_e-{\eta }_1{I}_e\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_l\left(t\right)& ={\gamma }_{w_e}{W}_e-{\lambda }_{w_l}{W}_l-{\gamma }_{w_l}{W}_l+(1-\alpha ){\gamma }_{i_e}{I}_e-{\eta }_2{I}_l\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_p\left(t\right)& ={\gamma }_{w_l}{W}_l-{\lambda }_{w_p}{W}_p-{\gamma }_{w_p}{W}_p+(1-\beta ){\gamma }_{i_l}{I}_l-{\eta }_3{I}_p\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_f\left(t\right)& =\rho {\gamma }_{w_p}{W}_p-{\lambda }_{w_f}{W}_f+(1-ϵ){\gamma }_{i_p}{\rho }_{i_w}{I}_p+{\eta }_4{W}_f\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_a\left(t\right)& =(1-\rho ){\gamma }_{w_p}{W}_p-{\lambda }_{w_a}{W}_a+(1-ϵ){\gamma }_{i_p}(1-{\rho }_{i_w})I_p+{\eta }_5{W}_a\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_e\left(t\right)& =\frac{{\Lambda }_{i_e}{I}_f(W_a+I_a)}{T}-{\lambda }_{i_e}{I}_e-\alpha {\gamma }_{i_e}{I}_e+{\eta }_1{I}_e\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_l\left(t\right)& =\alpha {\gamma }_{i_e}{I}_e-{\lambda }_{i_l}{I}_l-\beta {\gamma }_{i_l}{I}_l+{\eta }_2{I}_l\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_p\left(t\right)& =\beta {\gamma }_{i_l}{I}_l-{\lambda }_{i_p}{I}_p-ϵ{\gamma }_{i_p}{I}_p+{\eta }_3{I}_p\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_f\left(t\right)& ={\rho }_iϵ{\gamma }_{i_p}{I}_p-{\lambda }_{i_f}{I}_f-{\eta }_4{W}_f\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_a\left(t\right)& =(1-{\rho }_i)ϵ{\gamma }_{i_p}{I}_p-{\lambda }_{i_a}{I}_a-{\eta }_5{W}_a,\hfill \end{array}\right.$$

(6)

where ${\eta }_1,{\eta }_2$ and ${\eta }_3$ all are the rates at which the non-Wolbachia aquatic population gain Wolbachia infected mosquitoes infection and ${\eta }_4$ & ${\eta }_5$ are the rates at which the Wolbachia infected mosquitoes losses their Wolbachia infection.

Impulsive control plays an predominant role in dynamical systems such as Neural Networks [47,48], non–linear delay dynamic systems [49,50,51] and so forth. To optimize the Wolbachia release, we can release the Wolbachia infected eggs, larvae and pupae in the form of ’Zancu kit’ and Wolbachia infected adult female and male mosquitoes (introgression) impulsively. The situation should be monitored weekly once by Biogents trap (BG trap or BG sentinel trap). While monitoring, if there is less number of Wolbachia infected mosquitoes then in that situation we should release Wolbachia infected mosquitoes impulsively.

The mathematical model which describes the transmission dynamics of Wolbachia among Aedes Aegypti mosquitoes along with Wolbachia invasion and impulsive control is defined as follows:

When $t\ne t_{\theta }$ for $\theta =1,2,...m$,

$$\left\{\begin{array}{cc}\hfill {}_0^c{D}_t^{\alpha }{W}_e\left(t\right)& =\frac{{\Lambda }_{w_e}{W}_f{W}_a}{T}-{\lambda }_{w_e}{W}_e-{\gamma }_{w_e}{W}_e-{\eta }_1{I}_e\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_l\left(t\right)& ={\gamma }_{w_e}{W}_e-{\lambda }_{w_l}{W}_l-{\gamma }_{w_l}{W}_l+(1-\alpha ){\gamma }_{i_e}{I}_e-{\eta }_2{I}_l\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_p\left(t\right)& ={\gamma }_{w_l}{W}_l-{\lambda }_{w_p}{W}_p-{\gamma }_{w_p}{W}_p+(1-\beta ){\gamma }_{i_l}{I}_l-{\eta }_3{I}_p\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_f\left(t\right)& =\rho {\gamma }_{w_p}{W}_p-{\lambda }_{w_f}{W}_f+(1-ϵ){\gamma }_{i_p}{\rho }_{i_w}{I}_p+{\eta }_4{W}_f\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_a\left(t\right)& =(1-\rho ){\gamma }_{w_p}{W}_p-{\lambda }_{w_a}{W}_a+(1-ϵ){\gamma }_{i_p}(1-{\rho }_{i_w})I_p+{\eta }_5{W}_a\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_e\left(t\right)& =\frac{{\Lambda }_{i_e}{I}_f(W_a+I_a)}{T}-{\lambda }_{i_e}{I}_e-\alpha {\gamma }_{i_e}{I}_e+{\eta }_1{I}_e\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_l\left(t\right)& =\alpha {\gamma }_{i_e}{I}_e-{\lambda }_{i_l}{I}_l-\beta {\gamma }_{i_l}{I}_l+{\eta }_2{I}_l\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_p\left(t\right)& =\beta {\gamma }_{i_l}{I}_l-{\lambda }_{i_p}{I}_p-ϵ{\gamma }_{i_p}{I}_p+{\eta }_3{I}_p\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_f\left(t\right)& ={\rho }_iϵ{\gamma }_{i_p}{I}_p-{\lambda }_{i_f}{I}_f-{\eta }_4{W}_f\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_a\left(t\right)& =(1-{\rho }_i)ϵ{\gamma }_{i_p}{I}_p-{\lambda }_{i_a}{I}_a-{\eta }_5{W}_a\hfill \end{array}\right.$$

(7)

When $t=t_{\theta }$ for $\theta =1,2,...m$,

$$\left\{\begin{array}{c}\Delta W_e\left(t\right)=0\\ \Delta W_l\left(t\right)=0\\ \Delta W_p\left(t\right)=0\\ \Delta W_f\left(t\right)=0\\ \Delta W_a\left(t\right)=0\\ \Delta I_e\left(t\right)={\delta }_1{I}_e\left(t_{\theta }\right)\\ \Delta I_l\left(t\right)={\delta }_2{I}_l\left(t_{\theta }\right)\\ \Delta I_p\left(t\right)={\delta }_3{I}_p\left(t_{\theta }\right)\\ \Delta I_f\left(t\right)={\delta }_4{I}_f\left(t_{\theta }\right)\\ \Delta I_a\left(t\right)={\delta }_5{I}_a\left(t_{\theta }\right),\end{array}\right.$$

with initial conditions,

$$\begin{array}{ccc}\hfill W_e\left(t_0\right)& =& W_{e_0};W_l\left(t_0\right)=W_{l_0};W_{t_0}\left(0\right)=W_{p_0};W_{t_0}\left(0\right)=W_{f_0};W_{t_0}\left(0\right)=W_{a_0};\hfill \\ \hfill I_e\left(t_0\right)& =& I_{e_0};I_l\left(t_0\right)=I_{l_0};I_p\left(t_0\right)=I_{p_0};I_f\left(t_0\right)=I_{f_0};I_a\left(t_0\right)=I_{a_0};\hfill \end{array}$$

where $W_{e_0}$, $W_{l_0}$, $W_{p_0}$, $W_{f_0}$, $W_{a_0}$, $I_{e_0}$, $I_{l_0}$, $I_{p_0}$, $I_{f_0}$ and $I_{a_0}$ all are positive integers. Moreover,

$$\begin{array}{ccc}\hfill \Delta W_e\left(t_{\theta }\right)& =& W_e\left(t_{\theta }^{+}\right)-W_e\left(t_{\theta }^{-}\right)\hfill \\ \hfill \Delta W_l\left(t_{\theta }\right)& =& W_l\left(t_{\theta }^{+}\right)-W_l\left(t_{\theta }^{-}\right)\hfill \\ \hfill \Delta W_p\left(t_{\theta }\right)& =& W_p\left(t_{\theta }^{+}\right)-W_p\left(t_{\theta }^{-}\right)\hfill \\ \hfill \Delta W_f\left(t_{\theta }\right)& =& W_f\left(t_{\theta }^{+}\right)-W_f\left(t_{\theta }^{-}\right)\hfill \\ \hfill \Delta W_a\left(t_{\theta }\right)& =& W_a\left(t_{\theta }^{+}\right)-W_a\left(t_{\theta }^{-}\right)\hfill \\ \hfill \Delta I_e\left(t_{\theta }\right)& =& I_e\left(t_{\theta }^{+}\right)-I_e\left(t_{\theta }^{-}\right)\hfill \\ \hfill \Delta I_l\left(t_{\theta }\right)& =& I_l\left(t_{\theta }^{+}\right)-I_l\left(t_{\theta }^{-}\right)\hfill \\ \hfill \Delta I_p\left(t_{\theta }\right)& =& I_p\left(t_{\theta }^{+}\right)-I_p\left(t_{\theta }^{-}\right)\hfill \\ \hfill \Delta I_f\left(t_{\theta }\right)& =& I_f\left(t_{\theta }^{+}\right)-I_f\left(t_{\theta }^{-}\right)\hfill \\ \hfill \Delta I_a\left(t_{\theta }\right)& =& I_a\left(t_{\theta }^{+}\right)-I_a\left(t_{\theta }^{-}\right),\hfill \end{array}$$

with $t_1<t_2<t_3\cdots <t_m$. Let us assume that,

$$\begin{array}{ccc}\hfill M\left(t\right)& =& {\left[\begin{array}{cccccccccc}{W}_e\left(t\right)& W_l\left(t\right)& W_p\left(t\right)& W_f\left(t\right)& W_a\left(t\right)& I_e\left(t\right)& I_l\left(t\right)& I_p\left(t\right)& I_f\left(t\right)& I_a\left(t\right)\end{array}\right]}^{\top };\hfill \end{array}$$

$$\begin{array}{ccc}\hfill M^{*}& =& {\left[\begin{array}{cccccccccc}{W}_e^{*}& W_l^{*}& W_p^{*}& W_f^{*}& W_a^{*}& I_e^{*}& I_l^{*}& I_p^{*}& I_f^{*}& I_a^{*}\end{array}\right]}^{\top };\hfill \end{array}$$

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \Delta M\left(t_{\theta }\right)& =\hfill & \hfill \left[\begin{array}{ccccc}\Delta W_e\left(t\right)& \Delta W_l\left(t_{\theta }\right)& \Delta W_p\left(t_{\theta }\right)& \Delta W_f\left(t_{\theta }\right)& \Delta W_a\left(t_{\theta }\right)\end{array}\right.\\ & & \hfill {\left.\begin{array}{ccccc}\Delta I_e\left(t_{\theta }\right)& \Delta I_l\left(t_{\theta }\right)& \Delta I_p\left(t_{\theta }\right)& \Delta I_f\left(t_{\theta }\right)& \Delta I_a\left(t_{\theta }\right)\end{array}\right]}^{\top };\end{array}\end{array}$$

Therefore, (7) can be rewritten as,

$$\left\{\begin{array}{cc}{}_0^c{D}_t^{\alpha }M\left(t\right)\hfill & =-W_{1}M\left(t\right)+g\left(M\left(t\right)\right),t\ne t_{\theta },\theta =1,2,3,\cdots ,m\hfill \\ \Delta M\left(t_{\theta }\right)\hfill & =M\left(t_{\theta }^{+}\right)-M\left(t_{\theta }^{-}\right)={\delta }_{\theta }M\left(t_{\theta }\right),t=t_{\theta },\theta =1,2,3,\cdots ,m\hfill \\ M\left(t_0\right)\hfill & =M_0\in {\mathbb{Z}}^{+},\hfill \end{array}\right.$$

(8)

where,

$$\begin{array}{ccc}\hfill W_1& =\hfill & \hfill -\left[\begin{array}{ccccccc}{\lambda }_{w_e}+{\gamma }_{w_e}& 0& 0& 0& 0& {\eta }_1& 0\\ -{\gamma }_{w_e}& {\lambda }_{w_l}+{\gamma }_{w_l}& 0& 0& 0& -(1-\alpha ){\gamma }_{i_e}& {\eta }_2\\ 0& -{\gamma }_{w_l}& {\lambda }_{w_p}+{\gamma }_{w_p}& 0& 0& 0& -(1-\beta ){\gamma }_{i_l}\\ 0& 0& -\rho {\gamma }_{w_p}& {\lambda }_{w_f}-{\eta }_4& 0& 0& 0\\ 0& 0& -(1-\rho ){\gamma }_{w_p}& 0& {\lambda }_{w_a}-{\eta }_5& 0& 0\\ 0& 0& 0& 0& 0& {\lambda }_{i_e}+\alpha {\gamma }_{i_e}-{\eta }_1& 0\\ 0& 0& 0& 0& 0& -\alpha {\gamma }_{i_e}& {\lambda }_{i_l}+\beta {\gamma }_{i_l}-{\eta }_2\\ 0& 0& 0& 0& 0& 0& -\beta {\gamma }_{i_l}\\ 0& 0& 0& {\eta }_4& 0& 0& 0\\ 0& 0& 0& 0& {\eta }_5& 0& 0\end{array}\right.\\ \hfill & \hfill & \hfill \left.\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ {\eta }_3& 0& 0-(1-ϵ){\gamma }_{i_p}{\rho }_{i_w}& 0& 0\\ -(1-ϵ){\gamma }_{i_p}(1-{\rho }_{i_w})& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ {\lambda }_{i_p}+ϵ{\gamma }_{i_p}-{\eta }_3& 0& 0\\ -{\rho }_iϵ{\gamma }_{i_p}& {\lambda }_{i_f}& 0\\ -(1-{\rho }_iϵ{\gamma }_{i_p})& 0& {\lambda }_{i_a}\end{array}\right];\end{array}$$

$$\begin{array}{cc}\hfill g\left(M\right(t\left)\right)& =\left[\begin{array}{c}\frac{{\Lambda }_{w_e}{m}_4{m}_5}{T}\\ 0\\ 0\\ 0\\ 0\\ \frac{{\Lambda }_{i_e}{m}_9(m_5+m_{10})}{T}\\ 0\\ 0\\ 0\\ 0\end{array}\right];\hfill \end{array}\begin{array}{cc}\hfill {\delta }_{\theta }& =\left[\begin{array}{cccccccccc}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& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& {\delta }_1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& {\delta }_2& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& {\delta }_3& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& {\delta }_4& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& {\delta }_5\end{array}\right];\hfill \end{array}$$

6. Existence and Uniqueness of Solution

By utilizing the results from fixed point theory, the existence and uniqueness results for the system of Equation (7) were derived in this section.

Let $C_{n,m}=H$ be the Banach space of all bounded continuous function defined on $[n,m]\in \mathbb{R}$.

For the sake of simplicity, let

$$\left\{\begin{array}{cc}\hfill {}_0^c{D}_t^{\alpha }{W}_e\left(t\right)& =K_1(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_l\left(t\right)& =K_2(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_p\left(t\right)& =K_3(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_f\left(t\right)& =K_4(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{W}_a\left(t\right)& =K_5(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_e\left(t\right)& =K_6(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_l\left(t\right)& =K_7(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_p\left(t\right)& =K_8(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_f\left(t\right)& =K_9(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right)))\hfill \\ \hfill {}_0^c{D}_t^{\alpha }{I}_a\left(t\right)& =K_{10}(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right)).\hfill \end{array}\right.$$

(9)

where, $m_1\left(t\right)=W_e\left(t\right)$, $m_2\left(t\right)=W_l\left(t\right)$, $m_3\left(t\right)=W_p\left(t\right)$, $m_4\left(t\right)=W_f\left(t\right)$, $m_5\left(t\right)=W_a\left(t\right)$, $m_6\left(t\right)=I_e\left(t\right)$, $m_7\left(t\right)=I_l\left(t\right)$, $m_8\left(t\right)=I_p\left(t\right)$, $m_9\left(t\right)=I_f\left(t\right)$ and $m_{10}\left(t\right)=I_a\left(t\right)$. Moreover, let us assume that,

$$\left\{\begin{array}{cc}\hfill K_1(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& =\frac{{\Lambda }_{w_e}{W}_f{W}_a}{T}-{\lambda }_{w_e}{W}_e-{\gamma }_{w_e}{W}_e-{\eta }_1{I}_e\hfill \\ \hfill K_2(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& ={\gamma }_{w_e}{W}_e-{\lambda }_{w_l}{W}_l-{\gamma }_{w_l}{W}_l+(1-\alpha ){\gamma }_{i_e}{I}_e-{\eta }_2{I}_l\hfill \\ \hfill K_3(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& ={\gamma }_{w_l}{W}_l-{\lambda }_{w_p}{W}_p-{\gamma }_{w_p}{W}_p+(1-\beta ){\gamma }_{i_l}{I}_l-{\eta }_3{I}_p\hfill \\ \hfill K_4(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& =\rho {\gamma }_{w_p}{W}_p-{\lambda }_{w_f}{W}_f+(1-ϵ){\gamma }_{i_p}{\rho }_{i_w}{I}_p+{\eta }_4{W}_f\hfill \\ \hfill K_5(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& =(1-\rho ){\gamma }_{w_p}{W}_p-{\lambda }_{w_a}{W}_a+(1-ϵ){\gamma }_{i_p}(1-{\rho }_{i_w})I_p+{\eta }_5{W}_a\hfill \\ \hfill K_6(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& =\frac{{\Lambda }_{i_e}{I}_f(W_a+I_a)}{T}-{\lambda }_{i_e}{I}_e-\alpha {\gamma }_{i_e}{I}_e+{\eta }_1{I}_e\hfill \\ \hfill K_7(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& =\alpha {\gamma }_{i_e}{I}_e-{\lambda }_{i_l}{I}_l-\beta {\gamma }_{i_l}{I}_l+{\eta }_2{I}_l\hfill \\ \hfill K_8(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& =\beta {\gamma }_{i_l}{I}_l-{\lambda }_{i_p}{I}_p-ϵ{\gamma }_{i_p}{I}_p+{\eta }_3{I}_p\hfill \\ \hfill K_9(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& ={\rho }_iϵ{\gamma }_{i_p}{I}_p-{\lambda }_{i_f}{I}_f-{\eta }_4{W}_f\hfill \\ \hfill K_{10}(t,m_1\left(t\right),m_2\left(t\right),\cdots ,m_{10}\left(t\right))& =(1-{\rho }_i)ϵ{\gamma }_{i_p}{I}_p-{\lambda }_{i_a}{I}_a-{\eta }_5{W}_a.\hfill \end{array}\right.$$

(10)

By the Definition 3 of, fractional order anti derivative in Caputo sense, we have

$$\left\{\begin{array}{cc}\hfill W_e\left(t\right)-W_e\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_1(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill W_l\left(t\right)-W_l\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_2(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill W_p\left(t\right)-W_p\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_3(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill W_f\left(t\right)-W_f\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_4(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill W_a\left(t\right)-W_a\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_5(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \end{array}\right.$$

$$\left\{\begin{array}{cc}\hfill I_e\left(t\right)-I_e\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_6(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill I_l\left(t\right)-I_l\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_7(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill I_p\left(t\right)-I_p\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_8(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill I_f\left(t\right)-I_f\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_9(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill I_a\left(t\right)-I_a\left(0\right)& =\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_{10}(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta .\hfill \end{array}\right.$$

This implies that,

$$\left\{\begin{array}{cc}\hfill W_e\left(t\right)& =W_e\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_1(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill W_l\left(t\right)& =W_l\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_2(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill W_p\left(t\right)& =W_p\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_3(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill W_f\left(t\right)& =W_f\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_4(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill W_a\left(t\right)& =W_a\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_5(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill I_e\left(t\right)& =I_e\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_6(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill I_l\left(t\right)& =I_l\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_7(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill I_p\left(t\right)& =I_p\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_8(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill I_f\left(t\right)& =I_f\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_9(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta \hfill \\ \hfill I_a\left(t\right)& =I_a\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}{K}_{10}(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta .\hfill \end{array}\right.$$

(11)

Now, we define Equation (11) as

$$\begin{array}{c}\hfill M\left(t\right)=M\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}\underset{0}{\overset{t}{\int }}{(t-\eta )}^{\alpha -1}K(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta ,\end{array}$$

where

$$M\left(t\right)=\left(\begin{array}{c}{W}_e\left(t\right)\\ W_l\left(t\right)\\ W_p\left(t\right)\\ W_f\left(t\right)\\ W_a\left(t\right)\\ I_e\left(t\right)\\ I_l\left(t\right)\\ I_p\left(t\right)\\ I_f\left(t\right)\\ I_a\left(t\right)\end{array}\right)andK(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))=\left(\begin{array}{c}{K}_1(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\\ K_2(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\\ K_3(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\\ K_4(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\\ K_5(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\\ K_6(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\\ K_7(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\\ K_8(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\\ K_9(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\\ K_{10}(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\end{array}\right).$$

Let us define $C_{n,m}$ as $C_{n,m}={\mathfrak{F}}_n\left(t_0\right)\times H_m\left(s\right)$ where,

$$\begin{array}{c}\hfill s=min\{W_{e_0},W_{l_0},W_{p_0},W_{f_0},W_{a_0},I_{e_0},I_{l_0},I_{p_0},I_{f_0},I_{a_0}\}\end{array}$$

and

$$\begin{array}{c}\hfill {\mathfrak{F}}_n\left(t_0\right)=[t_0-n,t_0+n]\\ \hfill H_m\left(s\right)=[s-m,s+m].\end{array}$$

Along with this, we assumed that

$$\begin{array}{ccc}\hfill R& =& \underset{C_{n,m}}{max}\{\underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{1}}\left|\right|,\underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{2}}\left|\right|,\underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{3}}\left|\right|,\underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{4}}\left|\right|,\underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{5}}\left|\right|,\underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{6}}\left|\right|,\underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{7}}\left|\right|,\hfill \\ \hfill & & \underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{8}}\left|\right|,\underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{9}}\left|\right|,\underset{C_{n,m}}{sup}\left|\right|{\mathfrak{F}}_{\mathfrak{10}}\left|\right|\}.\hfill \end{array}$$

Let us define the norm at infinity as follows:

$${\left|\right|\mathsf{\Psi }\left|\right|}_{\infty }=\underset{t\in {\mathfrak{F}}_{\mathfrak{n}}}{sup}\left|\mathsf{\Psi }\left(t\right)\right|.$$

Here, the operator $\nu$: $C_{n,m}\to C_{n,m}$ is defined by

$$\nu \left(M\left(t\right)\right)=M\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\eta )}^{\alpha -1}K(\eta ,m_1\left(\eta \right),m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\mathrm{d}\eta .$$

(12)

To prove $\nu$ is well defined operator, we should prove that

$$\left|\right|\nu M\left(t\right)-{M\left(0\right)\left|\right|}_{\infty }<\left(\begin{array}{c}m\\ m\\ m\\ m\\ m\\ m\\ m\\ m\\ m\\ m\end{array}\right)$$

Now, let

$$\begin{array}{ccc}\hfill \left|\right|{\nu }_1{W}_e\left(t\right)-W_e\left(t\right){\left|\right|}_{\infty }& =& \left|\right|\frac{1}{\Gamma \left(\alpha \right)}\underset{0}{\overset{t}{\int }}{(t-\eta )}^{\alpha -1}{K}_1(\eta ,y_1\left(\eta \right),y_2\left(\eta \right),\cdots ,y_{10}\left(\eta \right))\mathrm{d}\eta {\left|\right|}_{\infty }\hfill \\ & \le & \frac{1}{\Gamma \left(\alpha \right)}\underset{0}{\overset{t}{\int }}{(t-\eta )}^{\alpha -1}\left|\right|K_1(\eta ,y_1\left(\eta \right),y_2\left(\eta \right),\cdots ,y_{10}\left(\eta \right)){\left|\right|}_{\infty }\mathrm{d}\eta \hfill \\ & \le & \frac{R}{\Gamma \left(\alpha \right)}\underset{0}{\overset{t}{\int }}{(t-\eta )}^{\alpha -1}\mathrm{d}\eta \hfill \\ & \le & \frac{R{n}^{\alpha }}{\Gamma (\alpha +1)},\hfill \end{array}$$

where,

$$n<{\left(\frac{m\Gamma (\alpha +1)}{R}\right)}^{1/n}.$$

As well as, we can prove that the other equations of (6) can satisfies this inequality.

That is, the operator $\nu$ is well-defined if

$$n<{\left(\frac{m\Gamma (\alpha +1)}{R}\right)}^{1/n}.$$

Now, we should prove that the operator $\nu$ satisfies the Lipschitz condition. That is,

$$\left|\right|{\nu }_{M_1}-{\nu }_{M_2}{\left|\right|}_{\infty }<h\left|\right|M_1-M_2{\left|\right|}_{\infty }$$

To prove this, let

$$\begin{array}{ccc}\hfill \left|\right|{\nu }_1{W}_{e_1}-{\nu }_1{W}_{e_2}\left|\right|& =& \left|\right|\frac{1}{\Gamma \left(\alpha \right)}\underset{0}{\overset{t}{\int }}{(t-\eta )}^{\alpha -1}{K}_1(W_{e_1},m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right)\eta )\mathrm{d}\eta \hfill \\ & & -\frac{1}{\Gamma \left(\alpha \right)}\underset{0}{\overset{t}{\int }}{(t-\eta )}^{\alpha -1}{K}_1(W_{e_2},m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right),\eta )\mathrm{d}\eta {\left|\right|}_{\infty }\hfill \\ & =& \frac{1}{\Gamma \left(\alpha \right)}\left|\right|\underset{0}{\overset{t}{\int }}{K}_1(\eta ,W_{e_1},m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\hfill \\ & & -K_1(\eta ,W_{e_2},m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right)){(t-\eta )}^{\alpha -1}\mathrm{d}\eta \left|\right|\hfill \\ & \le & \frac{1}{\Gamma \left(\alpha \right)}\underset{0}{\overset{t}{\int }}\left|\right|K_1(\eta ,W_{e_1},m_2\left(\eta \right),\cdots ,m_{10}\left(\eta \right))\hfill \\ & & -K_1(\eta ,W_{e_2},y_2,\cdots ,y_{10})\left|\right|{(t-\eta )}^{\alpha -1}\mathrm{d}\eta \hfill \\ & \le & \frac{1}{\Gamma \left(\alpha \right)}\underset{0}{\overset{t}{\int }}\left|\right|\frac{\omega }{N}\left(\mathfrak{F}\left(t\right)(W_{e_1}-W_{e_2})\right)\left|\right|{(t-\eta )}^{\alpha -1}\mathrm{d}\eta \hfill \\ & \le & \frac{n^{\alpha }{\left|\omega \right|\left|\right|\mathfrak{F}\left(t\right)\left|\right|}_{\infty }}{N\Gamma (\alpha +1)}\left|\right|W_{e_1}-W_{e_2}{\left|\right|}_{\infty }\hfill \\ & \le & h_1\left|\right|W_{e_1}-W_{e_2}{\left|\right|}_{\infty },\hfill \end{array}$$

with $h_1=\frac{n^{\alpha }{\left|\omega \right|\left|\right|\mathfrak{F}\left(t\right)\left|\right|}_{\infty }}{N\Gamma (\alpha +1)}$.

Similarly, we can prove that

$$\begin{array}{ccc}\hfill \left|\right|\nu W_{l_1}-\nu W_{l_2}\left|\right|& \le & h_2\left|\right|W_{l_1}-W_{l_2}\left|\right|\hfill \\ \hfill \left|\right|\nu W_{p_1}-\nu W_{p_2}\left|\right|& \le & h_3\left|\right|W_{p_1}-W_{p_2}\left|\right|\hfill \\ \hfill \left|\right|\nu W_{f_1}-\nu W_{f_2}\left|\right|& \le & h_4\left|\right|W_{f_1}-W_{f_2}\left|\right|\hfill \\ \hfill \left|\right|\nu W_{a_1}-\nu W_{a_2}\left|\right|& \le & h_5\left|\right|W_{a_1}-W_{a_2}\left|\right|\hfill \\ \hfill \left|\right|\nu I_{e_1}-\nu I_{e_2}\left|\right|& \le & h_6\left|\right|I_{e_1}-I_{e_2}\left|\right|\hfill \\ \hfill \left|\right|\nu I_{l_1}-\nu I_{l_2}\left|\right|& \le & h_7\left|\right|I_{l_1}-I_{l_2}\left|\right|\hfill \\ \hfill \left|\right|\nu I_{p_1}-\nu I_{p_2}\left|\right|& \le & h_8\left|\right|I_{p_1}-I_{p_2}\left|\right|\hfill \\ \hfill \left|\right|\nu I_{f_1}-\nu I_{f_2}\left|\right|& \le & h_9\left|\right|I_{f_1}-I_{f_2}\left|\right|\hfill \\ \hfill \left|\right|\nu I_{a_1}-\nu I_{a_2}\left|\right|& \le & h_{10}\left|\right|I_{a_1}-I_{a_2}\left|\right|.\hfill \end{array}$$

By the definition of Contraction mapping Definition 5, the map $\nu$ is a contraction map if $0<h_i<1$ for all $i=1,2,3,\cdots ,10$. Therefore, $\nu$ is a contraction mapping on a compact Banach space H. Then by Contraction mapping Theorem 1, $\nu$ has a solution and it is unique.

This implies that, the system of Equation (7) has a solution and its unique.

7. Stability Analysis

In the present section, the global Mittag-Leffler stability results were derived via LMI (Linear Matrix Inequality) approach and Lyapunov method.

Assumption (A1): Assume that the function $g\left(M\right(t\left)\right)$ satisfies the following:

For any $e_1,e_2\in {\mathbb{R}}^n$ there exists $S_1\in {\mathbb{R}}^{n\times n}$, such that $\left|\right|g\left(e_1\right)-g\left(e_2\right)\left|\right|\le \left|\right|S_1(e_1-e_2)\left|\right|$.

Theorem 2.

Assume that the system (8) satisfies the assumption (A1) and the impulsive operator satisfies that

$$\begin{array}{ccc}\hfill {\delta }_{\theta }\left(M\left(t_{\theta }\right)\right)& =& -\overline{\delta }(M\left(t_{\theta }\right)-M^{*}),\theta =1,2,\cdots ,m,\hfill \end{array}$$

where $M^{*}$ is an equilibrium point of system (8).

The system (8) is said to be globally Mittag-Leffler stable if there exists a positive definite matrix Q and positive scalars ξ and ${\gamma }_1$ such that the following inequalities hold:

$$\begin{array}{ccc}\hfill Q^{\frac{-1}{2}}{[{\gamma }_1+\overline{\delta }]}^{\top }Q[{\gamma }_1+\overline{\delta }]Q^{\frac{-1}{2}}& \le & {\gamma }_1\hfill \end{array}$$

(13)

and

$$\begin{array}{ccc}\hfill \tilde{\mathsf{\Omega }}& =& \left[\begin{array}{ccc}-2Q{W}_1& Q& \xi S_1\\ *& -\xi & 0\\ *& *& -\xi .\end{array}\right]<0.\hfill \end{array}$$

(14)

Proof. 

Let us consider the system (8) with the initial condition $M\left(t_0\right)=M_0\in {\mathbb{Z}}^{+}$ and an equilibrium point $M^{*}$. By using the transformation, $\mathcal{N}\left(t\right)=M\left(t\right)-M^{*}$, then the system (8) is transformed into

$$\begin{array}{ccc}\hfill {}_0^C{D}_t^{\alpha }\mathcal{N}\left(t\right)& =& -W_1\mathcal{N}\left(t\right)+\overline{g}\left(\mathcal{N}\left(t\right)\right),t\ne t_{\theta },\theta =1,2,3,...m\hfill \\ \hfill \Delta \mathcal{N}\left(t_{\theta }\right)& =& \mathcal{N}\left(t_{\theta }^{+}\right)-\mathcal{N}\left(t_{\theta }^{-}\right)=-\overline{\delta }\mathcal{N}\left(t_{\theta }\right),t=t_{\theta },\theta =1,2,3,...m\hfill \\ \hfill \mathcal{N}\left(t_0\right)& =& {\mathcal{N}}_0\in {\mathbb{Z}}^{+}.\hfill \end{array}$$

(15)

where, $\mathcal{N}\left(t\right)={({\mathcal{N}}_1,{\mathcal{N}}_2,{\mathcal{N}}_3,...,{\mathcal{N}}_{10})}^{\top }$ and $\overline{g}\left(\mathcal{N}\left(t\right)\right)={(\overline{g}\left({\mathcal{N}}_1\right),\overline{g}\left({\mathcal{N}}_2\right),..,\overline{g}\left({\mathcal{N}}_{10}\right))}^{\top }$ and ${\mathcal{N}}_0=M_0-M^{*}$. Let us consider a Lyapunov function as:

$$\begin{array}{ccc}\hfill V\left(t\right)& =& {\mathcal{N}}^{\top }\left(t\right)Q\mathcal{N}\left(t\right),\hfill \end{array}$$

(16)

where Q is a positive definite matrix. Now, the time derivative of $V\left(t\right)$ along with the trajectories of the system (16) is

$$\begin{array}{ccc}\hfill {}_0^C{D}_t^{\alpha }V\left(t\right)& \le & 2{\mathcal{N}}^{\top }\left(t\right)Q{}_0^C{D}_t^{\alpha }\mathcal{N}\left(t\right)\hfill \\ & =& {\mathcal{N}}^{\top }\left(t\right)2Q[-W_1\mathcal{N}\left(t\right)+\overline{g}\left(\mathcal{N}\left(t\right)\right)]\hfill \\ & =& {\mathcal{N}}^{\top }\left(t\right)(-2Q{W}_1)\mathcal{N}\left(t\right)+{\mathcal{N}}^{\top }\left(t\right)\left(2Q\right)\overline{g}\left(\mathcal{N}\left(t\right)\right)\hfill \end{array}$$

(17)

By Lemma 2,

$$\begin{array}{ccc}\hfill {\mathcal{N}}^{\top }\left(t\right)\left(2Q\right)\overline{g}\left(\mathcal{N}\left(t\right)\right)& \le & \frac{1}{\xi }{\mathcal{N}}^{\top }\left(t\right)\left(Q{Q}^{\top }\right)\mathcal{N}\left(t\right)+\xi {\overline{g}}^{\top }\left(\mathcal{N}\left(t\right)\right)\overline{g}\left(\mathcal{N}\left(t\right)\right).\hfill \end{array}$$

(18)

By assumption (A1),

$$\begin{array}{ccc}\hfill {\overline{g}}^{\top }\left(\mathcal{N}\left(t\right)\right)\overline{g}\left(\mathcal{N}\left(t\right)\right)& =& \langle \overline{g}\left(M\left(t\right)\right)-\overline{g}\left(M^{*}\right),\overline{g}\left(M\left(t\right)\right)-\overline{g}\left(M^{*}\right)\rangle \hfill \\ & =& \langle \left(\overline{g}(\mathcal{N}\left(t\right)+M^{*})\right)-\overline{g}\left(M^{*}\right),\overline{g}(\mathcal{N}\left(t\right)+M^{*})-\overline{g}\left(M^{*}\right)\rangle \hfill \\ & \le & {\mathcal{N}}^{\top }\left(t\right)S_1^{\top }{S}_1\mathcal{N}\left(t\right).\hfill \end{array}$$

(19)

Combine (18) and (19) and substitute in (17) we have,

$$\begin{array}{ccc}\hfill {}_0^C{D}_t^{\alpha }V\left(t\right)& \le & {\mathcal{N}}^{\top }\left(t\right)(-2Q{W}_1)\mathcal{N}\left(t\right)+{\mathcal{N}}^{\top }\left(t\right)\left({\xi }^{-1}Q{Q}^{\top }\right)\mathcal{N}\left(t\right)+{\mathcal{N}}^{\top }\left(t\right)\left(\xi S_1^{\top }{S}_1\right)\mathcal{N}\left(t\right)\hfill \\ & =& {\mathcal{N}}^{\top }\left(t\right)[-2Q{W}_1+{\xi }^{-1}Q{Q}^{\top }+\xi S_1^{\top }{S}_1]\mathcal{N}\left(t\right).\hfill \end{array}$$

(20)

Let, $\mathsf{\Omega }=-2Q{W}_1+{\xi }^{-1}Q{Q}^{\top }+\xi S_1^{\top }{S}_1$ and $\mathsf{\Omega }$ can be rewritten as

$$\begin{array}{ccc}\hfill \mathsf{\Omega }& =& \left[\begin{array}{ccc}-2Q{W}_1& Q& S_1\\ *& -\xi & 0\\ *& *& -{\xi }^{-1}\end{array}\right].\hfill \end{array}$$

Now, pre and post multiply $\mathsf{\Omega }$ by diag$\{I,I,\xi \}$, we get

$$\begin{array}{ccc}\hfill \tilde{\mathsf{\Omega }}& =& \left[\begin{array}{ccc}-2Q{W}_1& Q& \xi S_1\\ *& -\xi & 0\\ *& *& -\xi \end{array}\right].\hfill \end{array}$$

(21)

By Schur compliment Lemma 1, $\tilde{\mathsf{\Omega }}<0$.

Furthermore, the Equation (20), can be modified as

$$\begin{array}{ccc}\hfill {}_0^C{D}_t^{\alpha }V\left(t\right)& \le & {\mathcal{N}}^{\top }\left(t\right)\tilde{\mathsf{\Omega }}\mathcal{N}\left(t\right)\hfill \\ & =& -{\mathcal{N}}^{\top }\left(t\right)Q^{\frac{1}{2}}[-Q^{-\frac{1}{2}}\tilde{\mathsf{\Omega }}{Q}^{-\frac{1}{2}}]Q^{\frac{1}{2}}\mathcal{N}\left(t\right))\hfill \end{array}$$

let, ${ϵ}_1={\lambda }_{min}(-Q^{\frac{-1}{2}}\tilde{\mathsf{\Omega }}{Q}^{\frac{-1}{2}})$ and we know that $V\left(t\right)={\mathcal{N}}^{\top }\left(t\right)Q\mathcal{N}\left(t\right)$. This implies that,

$$\begin{array}{ccc}\hfill {}_0^C{D}_t^{\alpha }V\left(t\right)& \le & -{ϵ}_{1}V\left(t\right).\hfill \end{array}$$

(22)

For, $t_{\theta }=t$, $\theta =1,2,3,\cdots m$

$$\begin{array}{ccc}\hfill V\left(t_{\theta }^{+}\right)& =& {\mathcal{N}}^{\top }\left(t_{\theta }^{+}\right)Q\mathcal{N}\left(t_{\theta }^{+}\right)\hfill \\ & =& {[\mathcal{N}\left(t_{\theta }^{-}\right)+\overline{\delta }\mathcal{N}\left(t_{\theta }^{-}\right)]}^{\top }Q[\mathcal{N}\left(t_{\theta }^{-}\right)+\overline{\delta }\mathcal{N}\left(t_{\theta }^{-}\right)]\hfill \\ & =& {\mathcal{N}}^{\top }\left(t_{\theta }^{-}\right){[{\gamma }_1+\overline{\delta }]}^{\top }Q[{\gamma }_1+\overline{\delta }]\mathcal{N}\left(t_{\theta }^{-}\right)\hfill \\ & =& {\mathcal{N}}^{\top }\left(t_{\theta }^{-}\right)Q^{\frac{1}{2}}{[Q^{\frac{-1}{2}}{\gamma }_1+\overline{\delta }]}^{\top }Q[{\gamma }_1+\overline{\delta }{Q}^{\frac{-1}{2}}]Q^{\frac{1}{2}}\mathcal{N}\left(t_{\theta }^{-}\right)\hfill \\ & \le & {\mathcal{N}}^{\top }\left(t_{\theta }^{-}\right)Q\mathcal{N}\left(t_{\theta }^{-}\right)=V\left(\mathcal{N}\left(t_{\theta }^{-}\right)\right)\hfill \\ \hfill V\left(t_{\theta }^{+}\right)& \le & V\left(t_{\theta }^{-}\right)\hfill \end{array}$$

(23)

Therefore, we can easily prove that,

$$\begin{array}{ccc}\hfill {\lambda }_{min}{\left(Q\right)\left|\right|M\left(t\right)\left|\right|}^2& \le & V\left(t\right)\le {\lambda }_{max}{\left(Q\right)\left|\right|M\left(t\right)\left|\right|}^2.\hfill \end{array}$$

(24)

Conditions (22)–(24) satisfies the conditions of Lemma 3. Therefore by Lemma 3, our system (8) is globally Mittag-Leffler stable at its equilibrium point. □

8. Numerical Simulation

In this section, we provide an example to show the benefits of the proposed models (5)–(7). In this, we have analyzed three cases by published data mentioned in

Table 2. Data from published literature.

Parameters	Description	Data
${\Lambda }_{w_e}$	Reproduction rate of Wolbachia uninfected mosquitoes	1.25/day [52]
${\lambda }_{w_e}$, ${\lambda }_{w_l}$, ${\lambda }_{w_p}$	The death rate of aquatic Wolbachia uninfected mosquitoes	$\frac{1}{7.78}$/day [53]
${\gamma }_{w_e},{\gamma }_{w_l},{\gamma }_{w_p}$	The Maturation rate of Wolbachia uninfected mosquitoes	$\frac{1}{6.67}$/day [54]
${\lambda }_{w_f},{\lambda }_{w_a}$	The death rate of adult Wolbachia uninfected mosquitoes	$\frac{1}{14}$/day [53]
${\lambda }_{i_e},{\lambda }_{i_l},{\lambda }_{i_p}$	The death rate of aquatic Wolbachia infected mosquitoes	$\frac{1}{7.78}$ /day [53]
${\lambda }_{i_f},{\lambda }_{i_a}$	The death rate of adult Wolbachia infected mosquitoes	$\frac{1}{7}$/day [24]
${\Lambda }_{i_e}$	Reproduction rate of Wolbachia infected mosquitoes	$0.95\ast {\Lambda }_{w_e}/day$ [52]
${\gamma }_{i_e},{\gamma }_{i_l},{\gamma }_{i_p}$	The maturation rate of Wolbachia infected mosquitoes	$\frac{1}{6.67}$/day [24]

.

Case 1.

In this case, we have analyzed the transmission dynamics of Wolbachia among Aedes Aegypti mosquitoes via substituting the values mentioned in Table 2.

For this consider the system (5), with initial conditions $W_{e_0}=0.9$, $W_{l_0}=0.9$, $W_{p_0}=0.9$, $W_{f_0}=0.3$, $W_{a_0}=0.3$, $I_{e_0}=0.9$, $I_{l_0}=0.9$, $I_{p_0}=0.9$, $I_{f_0}=0.3$, $I_{a_0}=0.3$, total population $T=3000$, and the positive scalar used in Theorem 2 as $\xi =0.8513$

The Figure 4, Figure 5, Figure 6 and Figure 7 are depicts the dynamics of Equation (5) along with the parameters in Table 2 at various orders of $\alpha$ such as $\alpha =0.28,0.68,0.98$ and 1. We can observe by simulation results that, there is a notable decrease in non-Wolbachia mosquitoes and increase in Wolbachia infected mosquitoes.

Case 2.

In this case, we have analyzed the merits and demerits of considering the Wolbachia invasion. For this consider the system of Equation (6) with parameters mentioned in Table 2. We have plotted (6) with initial conditions and total population as considered in Case 1. Along with this, the other parameters ${\eta }_1=0.03$, ${\eta }_2=0.03$, ${\eta }_3=0.03$, ${\eta }_4=0.5$ and ${\eta }_1=0.5$ are fitted.

Figure 8, Figure 9, Figure 10 and Figure 11 are analyzed the dynamics of the system of Equation (6), with Wolbachia invasion and natural Wolbachia gain at various orders $\alpha =0.28,0.68,0.98$ and 1. From this we can observe that, Wolbachia infected mosquitoes tends to annihilation before the eradication of non-Wolbachia mosquitoes. It will lead to the decay in natural CI rescue.

Case 3.

In this case, the decay due to the natural Wolbachia invasion is managed by releasing Wolbachia infected mosquitoes impulsively. For this case, along with the parameters mentioned in Table 2, we have fitted the values of impulsive control as ${\delta }_1=0.4$, ${\delta }_2=0.4$, ${\delta }_3=0.3$, ${\delta }_4=0.5$ and ${\delta }_5=0.5$, invasion rates are ${\eta }_1=0.03$, ${\eta }_2=0.03$, ${\eta }_3=0.03$, and gain rates are ${\eta }_4=0.5$ and ${\eta }_1=0.5$.

Figure 12, Figure 13, Figure 14 and Figure 15 explicitly shows the dynamics of the systems of Equation (7) with impulsive control at orders $\alpha =0.28,0.68,0.98$ and 1. From this we get that, at order $\alpha =0.28$ the system leads to instability, when $\alpha =0.68$ the system started to posses stable state and at $\alpha =1$ the both population are annihilated at initial stage compared with Figure 7 and Figure 11.

By observing all the three cases, we can conclude that an impulsive control is an effective control strategy at Wolbachia invasion environment.

9. Conclusions

The effect of Wolbachia invasion and gain in vector population can lead to non-negligible in disease prevalence. Our impulsive control strategy shows that it is possible to control the transmission and invasion dynamics of Wolbachia bacterium. Our results shows that this method will increase the self-sustainability of Wolbachia bacterium among Aedes Aegypti mosquitoes. Another key result of the proposed fractional order model is, both mosquitoes population tends to annihilation after an impulsive controller synthesis. Further works on this model such as linearization, Lyapunov construction depicts that the created mathematical model is global Mittag-Leffler stable. In simulation performed here, depicts the effectiveness of the proposed model. In thus, we incorporated the real-world data from existing literature to compare the dynamical simulation of the 3 cases of model such as in the absence of Wolbachia invasion, the presence of Wolbachia invasion and the presence of Wolbachia invasion along with the impulsive control.