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

**Joseph Dianavinnarasi <sup>1</sup> , Ramachandran Raja <sup>2</sup> , Jehad Alzabut 3,\* , Michał Niezabitowski <sup>4</sup> and Ovidiu Bagdasar <sup>5</sup>**


**Abstract:** This work is devoted to analyzing an impulsive control synthesis to maintain the selfsustainability 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.

**Keywords:** sustainability; mosquito borne diseases; Aedes Aegypti; Wolbachia invasion; impulsive control

#### **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–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–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

**Citation:** Dianavinnarasi, J.; Raja, R.; Alzabut, J.; Niezabitowski, M.; Bagdasar, O. Controlling Wolbachia Transmission and Invasion Dynamics among Aedes Aegypti Population via Impulsive Control Strategy. *Symmetry* **2021**, *13*, 434. https://doi.org/ 10.3390/sym13030434

Academic Editors: Marin Marin and Jan Awrejcewicz

Received: 15 February 2021 Accepted: 3 March 2021 Published: 8 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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 femalespecific 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–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–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).

**Figure 1.** Mechanism of Wolbachia among mosquitoes and human.

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.

**Figure 2.** Block diagram representing the mechanism of Wolbachia infection in mosquitoes.

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:


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.** N denotes the space of all natural numbers, R denotes the space of all real numbers, C denotes the space of all complex numbers, R*<sup>n</sup>* denotes the space of *n*dimensional Euclidean space, Z <sup>+</sup> 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. <sup>∗</sup> denotes the corresponding symmetric terms in a symmetric matrix. Also, *<sup>c</sup> kD<sup>α</sup> t* (·) and *c k I α t* (·) denotes the derivative and anti derivative of order *α* with respect to *t* respectively, *c* denotes that its in Caputo sense, *k* denotes the initial condition and Γ(·) 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:*

$$\Gamma(z) \quad = \quad \int\_0^\infty e^{-s} s^{z-1} \,\mathbf{d} \, s \, \boldsymbol{z}$$

*with Re*(*z*) > 0*.*

**Definition 2.** *Ref. [1] The Caputo fractional derivative of a continuous function f*(*t*) *over* [*k*, *T*] *of order α* ∈ C *(with Re*(*α*) > 0*, α* ∈/ N) *is*

$$\begin{aligned} \, \_k^cD\_t^a f(t) \; &= \quad \_l\frac{1}{\Gamma(n-a)} \Big[ \int\_k^t (t-\eta)^{n-a-1} \frac{d^n}{d\eta^n} f(\eta) d\eta \Big] \, & \tag{1} \end{aligned} \tag{1}$$

*where, n* = [*Re*(*α*)] + 1*. If* 0 < *Re*(*α*) < 1*, then the expression* (1)*, can be rewritten as*

$$\, \_k^c D\_t^\alpha f(t) \quad = \quad \frac{1}{\Gamma(1-\alpha)} \left[ \int\_k^t \frac{f'(\eta)d\eta}{(t-\eta)^\alpha} \right]. \tag{2}$$

Since, *n* = 1 for all 0 < *Re*(*α*) < 1.

**Definition 3.** *Ref. [42] The Caputo sense fractional integral of a continuous function f on L* 1 ([0, *T*], R) *over α* ∈ (0, 1] *with respect to t is defined as*

$$\Gamma\_0^{\epsilon} I\_t^{\mu} f(t) \quad = \quad \frac{1}{\Gamma(\alpha)} \int\_0^t (t - \eta)^{\alpha - 1} f(\eta) d\eta. \tag{3}$$

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

$$E\_{a,b}(z) \quad = \quad \sum\_{l=0}^{\infty} \frac{z^l}{\Gamma(al+b)},$$

where, *z* ∈ C, *a* > 0, and *b* > 0. If *b* = 1 the *Ea*(*z*) = ∑ ∞ *l*=0 *z l* Γ(*al*+1) . If both *a* = 1 and *b* = 1, the *E*1,1(*z*) = *e z* .

**Lemma 1** (Schur Complement [43])**.** *Let us denote three n* × *n matrices as* Ψ1, Ψ2, Ψ3*, where* Ψ<sup>1</sup> = Ψ<sup>&</sup>gt; 1 *and* Ψ<sup>2</sup> = Ψ<sup>&</sup>gt; <sup>2</sup> > 0*. Then* Ψ<sup>1</sup> + Ψ<sup>&</sup>gt; <sup>3</sup> Ψ −1 <sup>2</sup> <sup>Ψ</sup><sup>3</sup> <sup>&</sup>lt; <sup>0</sup> *if and only if* Ψ<sup>1</sup> Ψ<sup>&</sup>gt; 3 Ψ<sup>3</sup> −Ψ<sup>2</sup> < 0 *(or)* −Ψ<sup>2</sup> Ψ<sup>3</sup> Ψ> <sup>3</sup> Ψ<sup>1</sup> < 0*.*

**Lemma 2.** *Ref. [44] For any scalar <sup>e</sup>* <sup>&</sup>gt; <sup>0</sup>*, A*, *<sup>N</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup> and matrix P*1*, then*

$$A^\top P\_1 N \quad \le \quad \frac{1}{2\epsilon} A^\top P\_1 P\_1^\top A + \frac{\epsilon}{2} N^\top N.$$

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

$$\begin{aligned} \mathbf{x}\_k^c D\_t^\mathbf{u} \mathbf{x}(t) &= -A\_1 \mathbf{x}(t) + A\_2 f(\mathbf{x}(t)), t \neq t\_\theta, \theta = 1, 2, \dots, m, \\ \Delta \mathbf{x}(t\_\theta) &= \mathbf{x}(t\_\theta^+) - \mathbf{x}(t\_\theta^-) = \delta\_\theta(\mathbf{x}(t\_\theta)), t = t\_\theta, \theta = 1, 2, \dots, m, \end{aligned} \tag{4}$$

with initial condition *x*(*t*0) = *x*<sup>0</sup> ∈ Z <sup>+</sup>, where the *n* states is defined by *<sup>x</sup>*(*t*) = [*x*1(*t*), *<sup>x</sup>*2(*t*), *<sup>x</sup>*3(*t*), · · · , *<sup>x</sup>n*(*t*)]<sup>&</sup>gt; <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* and *<sup>f</sup>*(*x*(*t*)) = [ *<sup>f</sup>*(*x*1(*t*)), *<sup>f</sup>*(*x*2(*t*)), *<sup>f</sup>*(*x*3(*t*)), · · · , *f*(*xn*(*t*))]<sup>&</sup>gt; be a function, *A*<sup>1</sup> and *A*<sup>2</sup> are constant coefficient matrices with the impulsive operator *δ<sup>θ</sup>* : <sup>R</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup>*<sup>n</sup>* .

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

$$||\mathfrak{x}(t) - \mathfrak{x}^\*|| \quad \leq \quad [h(\mathfrak{x}\_0 - \mathfrak{x}^\*)E\_\mathfrak{a}(-\mathfrak{x}t^a)]^b.$$

*where x* ∗ *is an equilibrium point,* 0 < *α* < 1*, κ* ≥ 0 *and a*, *b* > 0*. Moreover, h*(0) = 0*, h*(*x*) ≥ 0 *and h*(*x*) *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*(0) = 0*, t* > 0 *and δ<sup>θ</sup>* (0) = 0*, θ* = 1, 2, 3, · · · , *m. If there exists a positive definite function V such that the following hold:*

*(1) There exists positive constants α*<sup>1</sup> *and α*<sup>2</sup>

$$
\alpha\_1 ||\mathfrak{x}(t)|| \quad \le \quad V(t) \le \alpha\_2 ||\mathfrak{x}(t)|| / \varkappa(t) \in \mathbb{R}^n.
$$


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

**Definition 5.** *Ref. [46] A map ν* : *H* → *H, H compact Banach space, is said to be a contraction mapping if there exists h* ∈ (0, 1) *such that*

$$||\nu(m\_1) - \nu(m\_2)|| \quad \le \quad h||m\_1 - m\_2||$$

*for every m*1, *m*<sup>2</sup> ∈ *H.*

**Theorem 1** (Contraction Mapping Theorem)**.** *Ref. [46] Suppose H is a complete metric space and ν* : *H* → *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(*We*), non- Wolbachia larvae (*W<sup>l</sup>* ), non-Wolbachia pupae (*Wp*), non-Wolbachia adult female (*W<sup>f</sup>* ), non- Wolbachia adult male (*Wa*), Wolbachia infected eggs (*Ie*), Wolbachia infected larvae (*I<sup>l</sup>* ), Wolbachia infected pupae (*Ip*), Wolbachia infected adult female (*If* ), Wolbachia infected adult male (*Ia*). The total population at time *t* is denoted as *T* = *We*(*t*) + *W<sup>l</sup>* (*t*) + *Wp*(*t*) + *W<sup>f</sup>* (*t*) + *Wa*(*t*) + *Ie*(*t*) + *I<sup>l</sup>* (*t*) + *Ip*(*t*) + *I<sup>f</sup>* (*t*) + *Ia*(*t*). The eggs with zero Wolbachia infection are produced at the rate Λ*w<sup>e</sup>* by the mating process between non-Wolbachia female (*W<sup>f</sup>* ) and non-Wolbachia male (*Wa*). 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 <sup>Λ</sup>*weW<sup>f</sup> <sup>W</sup><sup>a</sup> T* . Along with this, the terms *λw<sup>e</sup>* (natural mortality rate of non-Wolbachia eggs) and *γw<sup>e</sup>* (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 *Ie*(*t*), depends on mating between Wolbachia infected female *I<sup>f</sup>* (*t*) and non-Wolbachia male *Wa*(*t*) and from mating between Wolbachia infected female *I<sup>f</sup>* (*t*) and Wolbachia infected male *Ia*(*t*). Through this, the birth rate of Wolbachia infected mosquito eggs population *Ie*(*t*) with the reproduction rate Λ*i<sup>e</sup>* is

$$\frac{\Lambda\_{i\_{\mathcal{C}}}(I\_f \mathcal{W}\_a + I\_f I\_a)}{T} \quad = \quad \frac{\Lambda\_{i\_{\mathcal{C}}} I\_f \left(\mathcal{W}\_a + I\_a\right)}{T}.$$

Similarly, the increase in the growth of Wolbachia infected eggs is limited by the natural mortality rate *λi<sup>e</sup>* and the maturation rate *γi<sup>e</sup>* (That is, the rate in which the corresponding compartment moved into the next stage).

Furthermore, the quantity (1 − *α*)*γi<sup>e</sup> Ie* is added to the wild mosquito larvae population. Because the term *α* and (1 − *α*) denotes the probability of getting larvae with and without Wolbachia respectively. Similarly, *β* and (1 − *β*) denotes the probability of getting pupae with and without Wolbachia respectively, *e* and (1 − *e*) denotes the probability rate of having Wolbachia infection in adult mosquitoes by introgression. That is, *e* be the probability of getting Wolbachia infected adults (with *ρi<sup>w</sup>* = probability of getting male and (1 − *ρi<sup>w</sup>* ) = probability of getting female). Because of these reasons, the terms (1 − *α*)*γi<sup>e</sup> Ie* , (1 − *β*)*γ<sup>i</sup> l Il* , (1 − *e*)*γi<sup>p</sup> ρiw I<sup>p</sup>* and (1 − *e*)*γi<sup>p</sup>* (1 − *ρi<sup>w</sup>* )*I<sup>p</sup>* are added to the corresponding stages and similarly, the terms *αγi<sup>e</sup> Ie* , *βγ<sup>i</sup> l I<sup>l</sup>* and *eγi<sup>p</sup> Iip* are removed from the

corresponding stages. The parameter description of the system of Equation (5) is presented in Table 1.

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


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

 *c* <sup>0</sup>*D α <sup>t</sup> W<sup>e</sup>* = Λ*weWfW<sup>a</sup> T* − *λweW<sup>e</sup>* − *γweW<sup>e</sup> c* <sup>0</sup>*D α <sup>t</sup> W<sup>l</sup>* = *γweW<sup>e</sup>* − *λwlW<sup>l</sup>* − *γwlW<sup>l</sup>* + (1 − *α*)*γi<sup>e</sup> Ie c* <sup>0</sup>*D α <sup>t</sup> W<sup>p</sup>* = *γwlW<sup>l</sup>* − *λwpW<sup>p</sup>* − *γwpW<sup>p</sup>* + (1 − *β*)*γ<sup>i</sup> l Il c* <sup>0</sup>*D α <sup>t</sup> W<sup>f</sup>* = *ργwpW<sup>p</sup>* − *λw<sup>f</sup> W<sup>f</sup>* + (1 − *e*)*γi<sup>p</sup> ρiw Ip c* <sup>0</sup>*D α <sup>t</sup> W<sup>a</sup>* = (1 − *ρ*)*γwpW<sup>p</sup>* − *λwaW<sup>a</sup>* + (1 − *e*)*γi<sup>p</sup>* (1 − *ρi<sup>w</sup>* )*I<sup>p</sup> c* <sup>0</sup>*D α t I<sup>e</sup>* = Λ*i<sup>e</sup> If* (*W<sup>a</sup>* + *Ia*) *T* − *λi<sup>e</sup> I<sup>e</sup>* − *αγi<sup>e</sup> Ie c* <sup>0</sup>*D α t I<sup>l</sup>* = *αγi<sup>e</sup> I<sup>e</sup>* − *λ<sup>i</sup> l I<sup>l</sup>* − *βγ<sup>i</sup> l Il c* <sup>0</sup>*D α t I<sup>p</sup>* = *βγ<sup>i</sup> l I<sup>l</sup>* − *λi<sup>p</sup> I<sup>p</sup>* − *eγi<sup>p</sup> Ip c* <sup>0</sup>*D α t I<sup>f</sup>* = *ρieγi<sup>p</sup> I<sup>p</sup>* − *λ<sup>i</sup> f If c* <sup>0</sup>*D α t I<sup>a</sup>* = (1 − *ρi*)*eγi<sup>p</sup> I<sup>p</sup>* − *λi<sup>a</sup> Ia*. (5)

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

**Figure 3.** Schematic representation Wolbachia spread dynamics among Aedes Aegypti mosquitoes.

#### **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*<sup>1</sup> = (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

2

$$P\_2 = (\mathcal{W}^\*\_{\varepsilon\_1 \prime} \mathcal{W}^\*\_{l\_1 \prime} \mathcal{W}^\*\_{p\_1 \prime} \mathcal{W}^\*\_{f\_1 \prime} \mathcal{W}^\*\_{a\_1 \prime}, 0, 0, 0, 0, 0)\_{\prime}$$

where,

$$\begin{array}{rcl} W\_{\varepsilon\_{1}}^{\*} &=& \frac{T\lambda\_{w\_{f}}\lambda\_{w\_{\varepsilon}}(\lambda\_{w\_{\varepsilon}}+\gamma\_{w\_{\varepsilon}})(\lambda\_{w\_{l}}+\gamma\_{w\_{l}})^{2}(\lambda\_{w\_{p}}+\gamma\_{w\_{p}})}{\rho(1-\rho)\Lambda\_{w\_{\varepsilon}}\gamma\_{w\_{p}}^{2}\gamma\_{w\_{\varepsilon}}^{2}\gamma\_{w\_{l}}^{2}}\\ W\_{l\_{1}}^{\*} &=& \frac{\gamma\_{w\_{\varepsilon}}}{\lambda\_{w\_{l}}+\gamma\_{w\_{l}}}W\_{\varepsilon\_{1}}^{\*}\\ W\_{p\_{1}}^{\*} &=& \frac{\gamma\_{w\_{l}}\gamma\_{w\_{\varepsilon}}}{(\lambda\_{w\_{l}}+\gamma\_{w\_{l}})(\lambda\_{w\_{p}}+\gamma\_{w\_{p}})}W\_{\varepsilon\_{1}}^{\*}\\ W\_{f\_{1}}^{\*} &=& \frac{\rho\gamma\_{w\_{p}}\gamma\_{w\_{\varepsilon}}\gamma\_{w\_{l}}}{\lambda\_{w\_{f}}(\lambda\_{w\_{f}}+\gamma\_{w\_{f}})(\lambda\_{w\_{p}}+\gamma\_{w\_{p}})}W\_{\varepsilon\_{1}}^{\*}\\ W\_{\tilde{a}\_{1}}^{\*} &=& \frac{(1-\rho)\gamma\_{w\_{p}}\gamma\_{w\_{\varepsilon}}}{\lambda\_{w\_{\varepsilon}}(\lambda\_{w\_{l}}+\gamma\_{w\_{l}})(\lambda\_{w\_{p}}+\gamma\_{w\_{p}})}W\_{\varepsilon\_{1}}^{\*}\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

$$P\_3 = (0,0,0,0,0,I\_{\mathfrak{e}\_2 \prime}^\* I\_{l\_2 \prime}^\* I\_{p\_2 \prime}^\* I\_{f\_2 \prime}^\* I\_{a\_2}^\*)\_{\prime\prime}$$

where,

$$\begin{split} I\_{e\_2}^\* &= \quad \frac{(\lambda\_{i\_l} + \beta \gamma\_{i\_l})(\lambda\_{i\_p} + \epsilon \gamma\_{i\_p})}{\alpha \beta \gamma\_{i\_l} \gamma\_{i\_l}} I\_{p\_2}^\* \\ I\_{l\_2}^\* &= \quad \frac{(\lambda\_{i\_p} + \epsilon \gamma\_{i\_p})}{\beta \gamma\_{i\_l}} I\_{p\_2}^\* \\ I\_{p\_2}^\* &= \quad \frac{T \lambda\_{i\_f} \lambda\_{i\_d} (\lambda\_{i\_c} + \alpha \gamma\_{i\_c})(\lambda\_{i\_l} + \beta \gamma\_{i\_l})(\lambda\_{i\_p} + \epsilon \gamma\_{i\_p})}{\lambda\_{i\_c} \alpha \beta \rho (1 - \rho\_l)\epsilon^2 \gamma\_{i\_p}^2 \gamma\_{i\_c} \gamma\_{i\_l}} \\ I\_{f\_2}^\* &= \quad \frac{\rho\_l \epsilon \gamma\_{i\_p}}{\lambda\_{i\_f}} I\_{p\_2}^\* \\ I\_{q\_2}^\* &= \quad \frac{(1 - \rho\_l)\epsilon \gamma\_{i\_p}}{\lambda\_{i\_q}} I\_{p\_2}^\*. \end{split}$$

*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

$$\mathcal{S}\_{\mathfrak{n}} \quad = \left\{ \mathcal{W}\_{\mathfrak{e}\_{n'}}^{\*} \mathcal{W}\_{l\_{n'}}^{\*} \mathcal{W}\_{p\_{n'}}^{\*} \mathcal{W}\_{f\_{n'}}^{\*} \mathcal{W}\_{a\_{n'}}^{\*} I\_{\mathfrak{e}\_{n'}}^{\*} I\_{l\_{n'}}^{\*} I\_{p\_{n'}}^{\*} I\_{f\_{n'}}^{\*} I\_{a\_{n}}^{\*} \right\}, n = 3, 4.$$

*W*∗ *<sup>e</sup><sup>n</sup>* = *λw<sup>l</sup>* + *γw<sup>l</sup> γw<sup>e</sup> λw<sup>p</sup>* + *γw<sup>p</sup> γw<sup>l</sup> TB*1*B*2*B*3*λ<sup>i</sup> f λw<sup>a</sup>* Λ*i<sup>e</sup>* (1 − *ρ*)*ρiγw<sup>p</sup>* ! − *I* ∗ *an γw<sup>e</sup> B*4(*λw<sup>l</sup>* + *γw<sup>l</sup>* ) *λw<sup>p</sup>* + *γw<sup>p</sup> γw<sup>l</sup>* + (*λw<sup>l</sup>* + *γw<sup>l</sup>* ) (1 − *β*)*γ<sup>i</sup> l γw<sup>l</sup> λi<sup>a</sup> B*1 *βγ<sup>i</sup> l* (1 − *ρi*) + (1 − *α*)*γi<sup>e</sup> λia B*1*B*<sup>2</sup> *αγi<sup>e</sup>* (1 − *ρi*) *W*∗ *ln* = *λw<sup>p</sup>* + *γw<sup>p</sup> γw<sup>l</sup>* " *TB*1*B*2*B*3*λ<sup>i</sup> f λw<sup>a</sup>* Λ*i<sup>e</sup>* (1 − *ρ*)*ρiγw<sup>p</sup>* − *B*<sup>4</sup> *I* ∗ *an* # − (1 − *β*)*γ<sup>i</sup> l γw<sup>l</sup> λi<sup>a</sup> B*1 *βγ<sup>i</sup> l* (1 − *ρi*) *W*∗ *<sup>p</sup><sup>n</sup>* = " *TB*1*B*2*B*3*λ<sup>i</sup> f λw<sup>a</sup>* Λ*i<sup>e</sup>* (1 − *ρ*)*ρiγw<sup>p</sup>* − *B*<sup>4</sup> *I* ∗ *an* # *W*∗ *fn* = *ργw<sup>p</sup> λw<sup>f</sup>* " *TB*1*B*2*B*3*λ<sup>i</sup> f λw<sup>a</sup>* Λ*i<sup>e</sup>* (1 − *ρ*)*ρiγw<sup>p</sup>* − *B*<sup>4</sup> *I* ∗ *an* # + (1 − *e*)*ρi<sup>w</sup> λia eλw<sup>f</sup>* (1 − *ρi*) *I* ∗ *an W*∗ *<sup>a</sup><sup>n</sup>* = *TB*1*B*2*B*3*λ<sup>i</sup> f ρi*Λ*i<sup>e</sup>* − *I* ∗ *an I* ∗ *<sup>e</sup><sup>n</sup>* = *B*1*B*2*λi<sup>a</sup> I* ∗ *an αγi<sup>e</sup>* (1 − *ρi*) *I* ∗ *<sup>p</sup><sup>n</sup>* = *λia* (1 − *ρi*)*eγi<sup>p</sup> I* ∗ *an I* ∗ *fn* = *ρiλi<sup>a</sup> I* ∗ *an* (1 − *ρi*)*λ<sup>i</sup> f*

with *I* ∗ *<sup>a</sup>*<sup>3</sup> > *I* ∗ *a*4 , both roots can be found from the quadratic equation

$$a\_1 I\_a^{\*^2} + a\_2 I\_a^\* + a\_3 \quad = \quad 0, 1$$

where,

$$\begin{array}{rcl} a\_{1} &=& \frac{\Lambda\_{\overline{w}\_{l}\theta}\theta\_{4}\gamma\_{w\_{p}}}{T\lambda\_{w\_{f}}};\\ a\_{2} &=& \left(\frac{\lambda\_{\overline{w}\_{l}}+\gamma\_{w\_{l}}}{T\lambda\_{w\_{f}}}\right)\left(\frac{\lambda\_{\overline{w}\_{l}\lambda\_{i}}\rho B\_{1}B\_{2}B\_{3}}{\rho\_{i}\Lambda\_{i}\lambda\_{w\_{f}}}\right)\left(\frac{\lambda\_{\overline{w}\_{l}}}{(1-\rho)}+B\_{4}\gamma\_{w\_{p}}\right)\\ & & \left(\frac{(\lambda\_{\overline{w}\_{l}}+\gamma\_{w\_{l}})(\lambda\_{\overline{w}\_{p}}+\gamma\_{w\_{p}})B\_{4}}{\gamma\_{w\_{l}}}+\frac{(\lambda\_{\overline{w}\_{l}}+\gamma\_{w\_{l}})(1-\beta)\lambda\_{i\_{d}}B\_{1}}{\gamma\_{w\_{l}}\beta(1-\rho\_{i})}+\frac{(1-\alpha)\lambda\_{i\_{d}}B\_{1}B\_{2}}{\alpha(1-\rho\_{i})}\right);\\ a\_{3} &=& \frac{\Lambda\_{\overline{w}\_{l}}\rhoTB\_{1}^{2}B\_{2}^{2}B\_{3}^{2}\lambda\_{w\_{l}}}{\Lambda\_{i\_{l}}^{2}(1-\rho)\rho\_{i}^{2}\lambda\_{w\_{f}}}.\end{array}$$

Here,

$$\begin{aligned} B\_1 &=& 1 + \frac{\lambda\_{i\_p}}{\varepsilon \gamma\_{i\_p}}; \\ B\_2 &=& 1 + \frac{\lambda\_{i\_l}}{\beta \gamma\_{i\_l}}; \\ B\_3 &=& 1 + \frac{\lambda\_{i\_c}}{\alpha \gamma\_{i\_c}}; \\ B\_4 &=& 1 + \frac{(1 - \varepsilon)(1 - \rho\_{i\_w})\lambda\_{i\_d}}{(1 - \rho)(1 - \rho\_i)\epsilon \gamma\_{w\_p}} \end{aligned}$$

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

 *c* <sup>0</sup>*D α <sup>t</sup> We*(*t*) = Λ*weWfW<sup>a</sup> T* − *λweW<sup>e</sup>* − *γweW<sup>e</sup>* − *η*<sup>1</sup> *I<sup>e</sup> c* <sup>0</sup>*D α <sup>t</sup> W<sup>l</sup>* (*t*) = *γweW<sup>e</sup>* − *λwlW<sup>l</sup>* − *γwlW<sup>l</sup>* + (1 − *α*)*γi<sup>e</sup> I<sup>e</sup>* − *η*<sup>2</sup> *I<sup>l</sup> c* <sup>0</sup>*D α <sup>t</sup> Wp*(*t*) = *γwlW<sup>l</sup>* − *λwpW<sup>p</sup>* − *γwpW<sup>p</sup>* + (1 − *β*)*γ<sup>i</sup> l I<sup>l</sup>* − *η*<sup>3</sup> *I<sup>p</sup> c* <sup>0</sup>*D α <sup>t</sup> W<sup>f</sup>* (*t*) = *ργwpW<sup>p</sup>* − *λw<sup>f</sup> W<sup>f</sup>* + (1 − *e*)*γi<sup>p</sup> ρiw I<sup>p</sup>* + *η*4*W<sup>f</sup> c* <sup>0</sup>*D α <sup>t</sup> Wa*(*t*) = (1 − *ρ*)*γwpW<sup>p</sup>* − *λwaW<sup>a</sup>* + (1 − *e*)*γi<sup>p</sup>* (1 − *ρi<sup>w</sup>* )*I<sup>p</sup>* + *η*5*W<sup>a</sup> c* <sup>0</sup>*D α t Ie*(*t*) = Λ*i<sup>e</sup> If* (*W<sup>a</sup>* + *Ia*) *T* − *λi<sup>e</sup> I<sup>e</sup>* − *αγi<sup>e</sup> I<sup>e</sup>* + *η*<sup>1</sup> *I<sup>e</sup> c* <sup>0</sup>*D α t Il* (*t*) = *αγi<sup>e</sup> I<sup>e</sup>* − *λ<sup>i</sup> l I<sup>l</sup>* − *βγ<sup>i</sup> l I<sup>l</sup>* + *η*<sup>2</sup> *I<sup>l</sup> c* <sup>0</sup>*D α t Ip*(*t*) = *βγ<sup>i</sup> l I<sup>l</sup>* − *λi<sup>p</sup> I<sup>p</sup>* − *eγi<sup>p</sup> I<sup>p</sup>* + *η*<sup>3</sup> *I<sup>p</sup> c* <sup>0</sup>*D α t If* (*t*) = *ρieγi<sup>p</sup> I<sup>p</sup>* − *λ<sup>i</sup> f I<sup>f</sup>* − *η*4*W<sup>f</sup> c* <sup>0</sup>*D α t Ia*(*t*) = (1 − *ρi*)*eγi<sup>p</sup> I<sup>p</sup>* − *λi<sup>a</sup> I<sup>a</sup>* − *η*5*Wa*, (6)

where *η*1, *η*<sup>2</sup> and *η*<sup>3</sup> all are the rates at which the non-Wolbachia aquatic population gain Wolbachia infected mosquitoes infection and *η*<sup>4</sup> & *η*<sup>5</sup> 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–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* 6= *t<sup>θ</sup>* for *θ* = 1, 2, ...*m*,

 *c* <sup>0</sup>*D α <sup>t</sup> We*(*t*) = Λ*weWfW<sup>a</sup> T* − *λweW<sup>e</sup>* − *γweW<sup>e</sup>* − *η*<sup>1</sup> *I<sup>e</sup> c* <sup>0</sup>*D α <sup>t</sup> W<sup>l</sup>* (*t*) = *γweW<sup>e</sup>* − *λwlW<sup>l</sup>* − *γwlW<sup>l</sup>* + (1 − *α*)*γi<sup>e</sup> I<sup>e</sup>* − *η*<sup>2</sup> *I<sup>l</sup> c* <sup>0</sup>*D α <sup>t</sup> Wp*(*t*) = *γwlW<sup>l</sup>* − *λwpW<sup>p</sup>* − *γwpW<sup>p</sup>* + (1 − *β*)*γ<sup>i</sup> l I<sup>l</sup>* − *η*<sup>3</sup> *I<sup>p</sup> c* <sup>0</sup>*D α <sup>t</sup> W<sup>f</sup>* (*t*) = *ργwpW<sup>p</sup>* − *λw<sup>f</sup> W<sup>f</sup>* + (1 − *e*)*γi<sup>p</sup> ρiw I<sup>p</sup>* + *η*4*W<sup>f</sup> c* <sup>0</sup>*D α <sup>t</sup> Wa*(*t*) = (1 − *ρ*)*γwpW<sup>p</sup>* − *λwaW<sup>a</sup>* + (1 − *e*)*γi<sup>p</sup>* (1 − *ρi<sup>w</sup>* )*I<sup>p</sup>* + *η*5*W<sup>a</sup> c* <sup>0</sup>*D α t Ie*(*t*) = Λ*i<sup>e</sup> If* (*W<sup>a</sup>* + *Ia*) *T* − *λi<sup>e</sup> I<sup>e</sup>* − *αγi<sup>e</sup> I<sup>e</sup>* + *η*<sup>1</sup> *I<sup>e</sup> c* <sup>0</sup>*D α t Il* (*t*) = *αγi<sup>e</sup> I<sup>e</sup>* − *λ<sup>i</sup> l I<sup>l</sup>* − *βγ<sup>i</sup> l I<sup>l</sup>* + *η*<sup>2</sup> *I<sup>l</sup> c* <sup>0</sup>*D α t Ip*(*t*) = *βγ<sup>i</sup> l I<sup>l</sup>* − *λi<sup>p</sup> I<sup>p</sup>* − *eγi<sup>p</sup> I<sup>p</sup>* + *η*<sup>3</sup> *I<sup>p</sup> c* <sup>0</sup>*D α t If* (*t*) = *ρieγi<sup>p</sup> I<sup>p</sup>* − *λ<sup>i</sup> f I<sup>f</sup>* − *η*4*W<sup>f</sup> c* <sup>0</sup>*D α t Ia*(*t*) = (1 − *ρi*)*eγi<sup>p</sup> I<sup>p</sup>* − *λi<sup>a</sup> I<sup>a</sup>* − *η*5*W<sup>a</sup>* (7)

When *t* = *t<sup>θ</sup>* for *θ* = 1, 2, ...*m*,

$$\begin{cases} \Delta W\_{\mathcal{C}}(t) = 0 \\ \Delta W\_{\mathcal{I}}(t) = 0 \\ \Delta W\_{\mathcal{P}}(t) = 0 \\ \Delta W\_{\mathcal{I}}(t) = 0 \\ \Delta W\_{\mathcal{I}}(t) = 0 \\ \Delta I\_{\mathcal{C}}(t) = \delta\_1 I\_{\mathcal{C}}(t\_\theta) \\ \Delta I\_{\mathcal{I}}(t) = \delta\_2 I\_{\mathcal{I}}(t\_\theta) \\ \Delta I\_{\mathcal{P}}(t) = \delta\_3 I\_{\mathcal{P}}(t\_\theta) \\ \Delta I\_{\mathcal{I}}(t) = \delta\_4 I\_{\mathcal{I}}(t\_\theta) \\ \Delta I\_{\mathcal{A}}(t) = \delta\_5 I\_{\mathcal{A}}(t\_\theta) \end{cases}$$

with initial conditions,

$$\begin{aligned} \mathcal{W}\_{\varepsilon}(t\_0) &= \mathcal{W}\_{\varepsilon\_0}; \mathcal{W}\_{l}(t\_0) = \mathcal{W}\_{l\_0}; \mathcal{W}\_{t\_0}(0) = \mathcal{W}\_{p\_0}; \mathcal{W}\_{t\_0}(0) = \mathcal{W}\_{f\_0}; \mathcal{W}\_{t\_0}(0) = \mathcal{W}\_{d\_0}; \\ I\_{\varepsilon}(t\_0) &= \ I\_{\varepsilon\_0}; I\_{l}(t\_0) = I\_{l\_0}; I\_{p}(t\_0) = I\_{p\_0}; I\_{f}(t\_0) = I\_{f\_0}; I\_{a}(t\_0) = I\_{a\_0}; \end{aligned}$$

where *We*<sup>0</sup> , *Wl*<sup>0</sup> , *Wp*<sup>0</sup> , *Wf*<sup>0</sup> , *Wa*<sup>0</sup> , *Ie*<sup>0</sup> , *Il*<sup>0</sup> , *Ip*<sup>0</sup> , *If*<sup>0</sup> and *Ia*<sup>0</sup> all are positive integers. Moreover,

$$\begin{array}{rcl}\Delta W\_{\mathfrak{c}}(t\_{\theta})&=&W\_{\mathfrak{c}}(t\_{\theta}^{+})-W\_{\mathfrak{c}}(t\_{\theta}^{-})\\\Delta W\_{\mathfrak{l}}(t\_{\theta})&=&W\_{\mathfrak{l}}(t\_{\theta}^{+})-W\_{\mathfrak{l}}(t\_{\theta}^{-})\\\Delta W\_{p}(t\_{\theta})&=&W\_{p}(t\_{\theta}^{+})-W\_{p}(t\_{\theta}^{-})\\\Delta W\_{f}(t\_{\theta})&=&W\_{f}(t\_{\theta}^{+})-W\_{f}(t\_{\theta}^{-})\\\Delta W\_{\mathfrak{a}}(t\_{\theta})&=&W\_{\mathfrak{a}}(t\_{\theta}^{+})-W\_{\mathfrak{a}}(t\_{\theta}^{-})\\\Delta I\_{\mathfrak{c}}(t\_{\theta})&=&I\_{\mathfrak{c}}(t\_{\theta}^{+})-I\_{\mathfrak{c}}(t\_{\theta}^{-})\\\Delta I\_{\mathfrak{l}}(t\_{\theta})&=&I\_{\mathfrak{l}}(t\_{\theta}^{+})-I\_{\mathfrak{l}}(t\_{\theta}^{-})\\\Delta I\_{\mathfrak{p}}(t\_{\theta})&=&I\_{\mathfrak{p}}(t\_{\theta}^{+})-I\_{\mathfrak{p}}(t\_{\theta}^{-})\\\Delta I\_{f}(t\_{\theta})&=&I\_{f}(t\_{\theta}^{+})-I\_{f}(t\_{\theta}^{-})\\\Delta I\_{\mathfrak{d}}(t\_{\theta})&=&I\_{\mathfrak{d}}(t\_{\theta}^{+})-I\_{\mathfrak{d}}(t\_{\theta}^{-})\\\end{array}$$

with *t*<sup>1</sup> < *t*<sup>2</sup> < *t*<sup>3</sup> · · · < *tm*. Let us assume that,

$$\begin{array}{rclclcl}\Delta M(t) &=& \begin{bmatrix} \mathcal{W}\_{\varepsilon}(t) & \mathcal{W}\_{l}(t) & \mathcal{W}\_{p}(t) & \mathcal{W}\_{f}(t) & \mathcal{W}\_{a}(t) & I\_{\varepsilon}(t) & I\_{l}(t) & I\_{p}(t) & I\_{l}(t) & I\_{a}(t) \end{bmatrix}^{\intercal}; \\\\ \\ &M^{\*} &=& \begin{bmatrix} \mathcal{W}\_{\varepsilon}^{\*} & \mathcal{W}\_{l}^{\*} & \mathcal{W}\_{p}^{\*} & \mathcal{W}\_{f}^{\*} & \mathcal{W}\_{a}^{\*} & I\_{\varepsilon}^{\*} & I\_{l}^{\*} & I\_{p}^{\*} & I\_{f}^{\*} & I\_{a}^{\*} \end{bmatrix}^{\intercal}; \\\\ \\ \Delta M(t\_{\theta}) &=& \begin{bmatrix} \Delta \mathcal{W}\_{\varepsilon}(t) & \Delta \mathcal{W}\_{l}(t\_{\theta}) & \Delta \mathcal{W}\_{p}(t\_{\theta}) & \Delta \mathcal{W}\_{f}(t\_{\theta}) & \Delta \mathcal{W}\_{a}(t\_{\theta}) \\\\ \Delta I\_{\varepsilon}(t\_{\theta}) & \Delta I\_{l}(t\_{\theta}) & \Delta I\_{p}(t\_{\theta}) & \Delta I\_{f}(t\_{\theta}) & \Delta I\_{d}(t\_{\theta}) \end{bmatrix}^{\intercal}; \end{array}$$

Therefore, (7) can be rewritten as,

$$\begin{cases} ^c\_0 D\_t^\mu M(t) &= -\mathcal{W}\_1 M(t) + \mathcal{g}(M(t)), t \neq t\_\theta, \theta = 1, 2, 3, \cdots \text{, } m \\ ^c\_1 \Delta M(t\_\theta) &= M(t\_\theta^+) - M(t\_\theta^-) = \delta\_\theta M(t\_\theta), t = t\_\theta, \theta = 1, 2, 3, \cdots \text{, } m \\ M(t\_0) &= M\_0 \in \mathbb{Z}^+ \end{cases} \tag{8}$$

where,

*W*<sup>1</sup> = − *λw<sup>e</sup>* + *γw<sup>e</sup>* 0 0 0 0 *η*<sup>1</sup> 0 −*γw<sup>e</sup> λw<sup>l</sup>* + *γw<sup>l</sup>* 0 0 0 −(1 − *α*)*γi<sup>e</sup> η*2 0 −*γw<sup>l</sup> λw<sup>p</sup>* + *γw<sup>p</sup>* 0 0 0 −(1 − *β*)*γ<sup>i</sup> l* 0 0 −*ργw<sup>p</sup> λw<sup>f</sup>* − *η*<sup>4</sup> 0 0 0 0 0 −(1 − *ρ*)*γw<sup>p</sup>* 0 *λw<sup>a</sup>* − *η*<sup>5</sup> 0 0 0 0 0 0 0 *λi<sup>e</sup>* + *αγi<sup>e</sup>* − *η*<sup>1</sup> 0 0 0 0 0 0 −*αγi<sup>e</sup> λi <sup>l</sup>* + *βγ<sup>i</sup> <sup>l</sup>* − *η*<sup>2</sup> 0 0 0 0 0 0 −*βγ<sup>i</sup> l* 0 0 0 *η*<sup>4</sup> 0 0 0 0 0 0 0 *η*<sup>5</sup> 0 0 0 0 0 0 0 0 *η*<sup>3</sup> 0 0 − (1 − *e*)*γi<sup>p</sup> ρiw* 0 0 −(1 − *e*)*γi<sup>p</sup>* (1 − *ρi<sup>w</sup>* ) 0 0 0 0 0 0 0 0 *λi<sup>p</sup>* + *eγi<sup>p</sup>* − *η*<sup>3</sup> 0 0 −*ρieγi<sup>p</sup> λi f* 0 −(1 − *ρieγi<sup>p</sup>* ) 0 *λi<sup>a</sup>* ; Λ*wem*4*m*<sup>5</sup> 

$$\begin{aligned} \mathcal{g}(M(t)) = \begin{bmatrix} \frac{\Delta\mathcal{E}\_{\theta}\circ\mathcal{E}\_{\theta}\circ\mathcal{E}\_{\theta}}{T} \\ 0 \\ 0 \\ 0 \\ 0 \\ \frac{\Delta\_{\mathcal{E}}\circ\mathcal{M}(\mathfrak{m}\_{5}+\mathfrak{m}\_{10})}{T} \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}; \quad \delta\_{\theta} = \begin{bmatrix} 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{bmatrix}; \end{aligned}$$

#### **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 *Cn*,*<sup>m</sup>* = *H* be the Banach space of all bounded continuous function defined on [*n*, *m*] ∈ R.

For the sake of simplicity, let


where, *m*1(*t*) = *We*(*t*), *m*2(*t*) = *W<sup>l</sup>* (*t*), *m*3(*t*) = *Wp*(*t*), *m*4(*t*) = *W<sup>f</sup>* (*t*), *m*5(*t*) = *Wa*(*t*), *m*6(*t*) = *Ie*(*t*), *m*7(*t*) = *I<sup>l</sup>* (*t*), *m*8(*t*) = *Ip*(*t*), *m*9(*t*) = *I<sup>f</sup>* (*t*) and *m*10(*t*) = *Ia*(*t*). Moreover, let us assume that,

$$\begin{cases} \mathcal{K}\_{1}(t,m\_{1}(t),m\_{2}(t),\ldots,m\_{10}(t)) &= \frac{\Lambda\_{W\_{f}}W\_{f}\mathcal{W}\_{e}}{T} - \lambda\_{w\_{i}}\mathcal{W}\_{e} - \gamma\_{w\_{i}}\mathcal{W}\_{e} - \eta\_{1}I\_{e} \\ \mathcal{K}\_{2}(t,m\_{1}(t),m\_{2}(t),\ldots,m\_{10}(t)) &= \gamma\_{w\_{i}}W\_{e} - \lambda\_{W\_{i}}W\_{\mathcal{W}} - \gamma\_{w\_{i}}W\_{\mathcal{W}} + (1-a)\gamma\_{\bar{\chi}\_{i}}I\_{\mathcal{L}} - \eta\_{2}I\_{i} \\ \mathcal{K}\_{3}(t,m\_{1}(t),m\_{2}(t),\ldots,m\_{10}(t)) &= \gamma\_{\bar{\chi}\_{i}}W\_{\mathcal{I}} - \lambda\_{\bar{\mathcal{W}}\_{e}}W\_{\mathcal{W}} - \gamma\_{\bar{\mathcal{W}}\_{e}}W\_{\mathcal{V}} + (1-\bar{\rho})\gamma\_{\bar{\mathcal{W}}\_{i}}I\_{\mathcal{U}} - \eta\_{3}I\_{f} \\ \mathcal{K}\_{4}(t,m\_{1}(t),m\_{2}(t),\ldots,m\_{10}(t)) &= \rho\gamma\_{\bar{\mathcal{W}}\_{e}}W\_{\mathcal{V}} - \lambda\_{\bar{\mathcal{W}}\_{e}}W\_{\mathcal{V}} + (1-\epsilon)\gamma\_{\bar{\mathcal{W}}\_{i}}I\_{\mathcal{U}} + \eta\_{1}W\_{f} \\ \mathcal{K}\_{5}(t,m\_{1}(t),m\_{2}(t),\ldots,m\_{10}(t)) &= (1-\rho)\gamma\_{\bar{\mathcal{W}}\_{e}}W\_{\mathcal{V}} - \lambda\_{\bar{\mathcal{W}}\_{e}}W\_{\mathcal{$$

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

$$\begin{cases} \begin{aligned} \mathcal{W}\_{\varepsilon}(t) - \mathcal{W}\_{\varepsilon}(0) &= \frac{1}{\Gamma(\alpha)} \int\_{0}^{t} (t-\eta)^{\alpha-1} \mathcal{K}\_{1}(\eta, m\_{1}(\eta), m\_{2}(\eta), \dots, m\_{10}(\eta)) d\eta \\ \mathcal{W}\_{l}(t) - \mathcal{W}\_{l}(0) &= \frac{1}{\Gamma(\alpha)} \int\_{0}^{t} (t-\eta)^{\alpha-1} \mathcal{K}\_{2}(\eta, m\_{1}(\eta), m\_{2}(\eta), \dots, m\_{10}(\eta)) d\eta \\ \mathcal{W}\_{p}(t) - \mathcal{W}\_{p}(0) &= \frac{1}{\Gamma(\alpha)} \int\_{0}^{t} (t-\eta)^{\alpha-1} \mathcal{K}\_{3}(\eta, m\_{1}(\eta), m\_{2}(\eta), \dots, m\_{10}(\eta)) d\eta \\ \mathcal{W}\_{f}(t) - \mathcal{W}\_{f}(0) &= \frac{1}{\Gamma(\alpha)} \int\_{0}^{t} (t-\eta)^{\alpha-1} \mathcal{K}\_{4}(\eta, m\_{1}(\eta), m\_{2}(\eta), \dots, m\_{10}(\eta)) d\eta \\ \mathcal{W}\_{a}(t) - \mathcal{W}\_{a}(0) &= \frac{1}{\Gamma(\alpha)} \int\_{0}^{t} (t-\eta)^{\alpha-1} \mathcal{K}\_{5}(\eta, m\_{1}(\eta), m\_{2}(\eta), \dots, m\_{10}(\eta)) d\eta \end{aligned} \end{cases}$$

$$\begin{cases} I\_{\varepsilon}(t) - I\_{\varepsilon}(0) &= \frac{1}{\Gamma(\mathfrak{a})} \int\_{0}^{t} (t-\eta)^{\mathfrak{a}-1} \mathcal{K}\_{\mathbb{G}}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdot, \cdot, m\_{10}(\eta)) \mathrm{d}\eta \\ I\_{l}(t) - I\_{l}(0) &= \frac{1}{\Gamma(\mathfrak{a})} \int\_{0}^{t} (t-\eta)^{\mathfrak{a}-1} \mathcal{K}\_{\mathbb{G}}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdot, \cdot, m\_{10}(\eta)) \mathrm{d}\eta \\ I\_{p}(t) - I\_{p}(0) &= \frac{1}{\Gamma(\mathfrak{a})} \int\_{0}^{t} (t-\eta)^{\mathfrak{a}-1} \mathcal{K}\_{\mathbb{G}}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdot, \cdot, m\_{10}(\eta)) \mathrm{d}\eta \\ I\_{f}(t) - I\_{f}(0) &= \frac{1}{\Gamma(\mathfrak{a})} \int\_{0}^{t} (t-\eta)^{\mathfrak{a}-1} \mathcal{K}\_{\mathbb{G}}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdot, \cdot, m\_{10}(\eta)) \mathrm{d}\eta \\ I\_{a}(t) - I\_{a}(0) &= \frac{1}{\Gamma(\mathfrak{a})} \int\_{0}^{t} (t-\eta)^{\mathfrak{a}-1} \mathcal{K}\_{\mathbb{G}}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdot, \cdot, m\_{10}(\eta)) \mathrm{d}\eta. \end{cases}$$

This implies that,

 *<sup>W</sup>e*(*t*) = *<sup>W</sup>e*(0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*1(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*<sup>η</sup> Wl* (*t*) = *W<sup>l</sup>* (0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*2(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*<sup>η</sup> <sup>W</sup>p*(*t*) = *<sup>W</sup>p*(0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*3(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*<sup>η</sup> W<sup>f</sup>* (*t*) = *W<sup>f</sup>* (0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*4(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*<sup>η</sup> <sup>W</sup>a*(*t*) = *<sup>W</sup>a*(0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*5(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*<sup>η</sup> <sup>I</sup>e*(*t*) = *<sup>I</sup>e*(0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*6(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*<sup>η</sup> Il* (*t*) = *I<sup>l</sup>* (0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*7(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*<sup>η</sup> <sup>I</sup>p*(*t*) = *<sup>I</sup>p*(0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*8(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*<sup>η</sup> If* (*t*) = *I<sup>f</sup>* (0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*9(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*<sup>η</sup> <sup>I</sup>a*(*t*) = *<sup>I</sup>a*(0) + <sup>1</sup> Γ(*α*) Z *t* 0 (*t* − *η*) *<sup>α</sup>*−1*K*10(*η*, *<sup>m</sup>*1(*η*), *<sup>m</sup>*2(*η*), · · · , *<sup>m</sup>*10(*η*))d*η*.

Now, we define Equation (11) as

$$M(t) = M(0) + \frac{1}{\Gamma(\mathfrak{a})} \int\_0^t (t - \eta)^{a-1} K(\eta, m\_1(\eta), m\_2(\eta), \dots, m\_{10}(\eta)) \mathrm{d}\eta \,\omega$$

$$\mathcal{M}(t) = \begin{pmatrix} W\_{\varepsilon}(t) \\ W\_{\eta}(t) \\ W\_{p}(t) \\ W\_{p}(t) \\ W\_{a}(t) \\ I\_{\varepsilon}(t) \\ I\_{\varrho}(t) \\ I\_{\varrho}(t) \\ I\_{\varrho}(t) \end{pmatrix} and \ K(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) = \begin{pmatrix} K\_{1}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) \\ K\_{2}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) \\ K\_{3}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) \\ K\_{4}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) \\ K\_{5}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) \\ K\_{6}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) \\ K\_{7}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) \\ K\_{8}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) \\ K\_{9}(\eta, m\_{1}(\eta), m\_{2}(\eta), \cdots, m\_{10}(\eta)) \\ \end{pmatrix}$$

Let us define *Cn*,*<sup>m</sup>* as *Cn*,*<sup>m</sup>* = F*n*(*t*0) × *Hm*(*s*) where,

$$s = \min\{\mathcal{W}\_{\mathfrak{e}\_{0'}}, \mathcal{W}\_{l\_{0'}}, \mathcal{W}\_{p\_{0'}}, \mathcal{W}\_{f\_{0'}}, \mathcal{W}\_{a\_{0'}}, I\_{\mathfrak{e}\_{0'}}, I\_{p\_{0'}}, I\_{f\_{0'}}, I\_{a\_{0}}\}$$

and

$$\begin{aligned} \mathfrak{F}\_n(t\_0) &= [t\_0 - n, t\_0 + n], \\ H\_m(s) &= [s - m, s + m]. \end{aligned}$$

Along with this, we assumed that

$$\begin{array}{ccl} R & = & \max\{\sup ||\mathfrak{F}\_{1}||, \sup ||\mathfrak{F}\_{2}||, \sup ||\mathfrak{F}\_{3}||, \sup ||\mathfrak{F}\_{3}||, \sup ||\mathfrak{F}\_{4}||, \sup ||\mathfrak{F}\_{5}||, \sup ||\mathfrak{F}\_{6}||, \sup ||\mathfrak{F}\_{7}||\right\} \\ & \sup ||\mathfrak{F}\_{8}|| / \sup ||\mathfrak{F}\_{9}||, \sup ||\mathfrak{F}\_{9}||, \sup ||\mathfrak{F}\_{10}||, \\ & \mathbb{C}\_{\mathfrak{n},\mathfrak{u}} \quad \mathbb{C}\_{\mathfrak{n},\mathfrak{u}} \quad \mathbb{C}\_{\mathfrak{n},\mathfrak{u}} \end{array}$$

Let us define the norm at infinity as follows:

$$||\Psi||\_{\infty} = \sup\_{t \in \mathfrak{F}\_{\mathbf{n}}} |\Psi(t)|.$$

Here, the operator *ν*: *Cn*,*<sup>m</sup>* → *Cn*,*<sup>m</sup>* is defined by

$$\nu(M(t)) = M(0) + \frac{1}{\Gamma(a)} \int\_0^t (t - \eta)^{a-1} K(\eta, m\_1(\eta), m\_2(\eta), \dots, m\_{10}(\eta)) d\eta. \tag{12}$$

To prove *ν* is well defined operator, we should prove that

$$||\iota M(t) - M(0)||\_{\infty} < \begin{pmatrix} m \\ m \\ m \\ m \\ m \\ m \\ m \\ m \\ m \\ m \end{pmatrix}$$

Now, let

$$\begin{split} ||\boldsymbol{\nu}\_{1}\mathbf{W}\_{\varepsilon}(t) - \mathbf{W}\_{\varepsilon}(t)||\_{\infty} &= \mathop{\mathrm{||}\limits\_{\boldsymbol{\Gamma}(\alpha)} \int\_{0}^{t} (t-\eta)^{\alpha-1} \mathbf{K}\_{1}(\boldsymbol{\eta}, \mathbf{y}\_{1}(\eta), \mathbf{y}\_{2}(\eta), \dots, \mathbf{y}\_{10}(\eta)) \mathrm{d}\boldsymbol{\eta}||\_{\infty}}{\displaystyle} \\ &\leq \mathop{\mathrm{||}\limits\_{\boldsymbol{\Gamma}(\alpha)} \int\_{0}^{t} (t-\eta)^{\alpha-1} ||\mathbf{K}\_{1}(\boldsymbol{\eta}, \mathbf{y}\_{1}(\eta), \mathbf{y}\_{2}(\eta), \dots, \mathbf{y}\_{10}(\eta))||\_{\infty} \mathrm{d}\boldsymbol{\eta} \\ &\leq \mathop{\mathrm{||}\limits\_{\boldsymbol{\Gamma}(\alpha)} \int\_{0}^{t} (t-\eta)^{\alpha-1} \mathbf{d}\boldsymbol{\eta} \\ &\leq \mathop{\mathrm{||}\limits\_{\boldsymbol{\Gamma}(\alpha+1)} \int\_{0}^{t} \mathbf{K}\mathbf{r}^{\alpha}} \,. \end{split}$$

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 *ν* is well-defined if

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

Now, we should prove that the operator *ν* satisfies the Lipschitz condition. That is,

.

$$||\nu\_{M\_1} - \nu\_{M\_2}||\_{\infty} < h||M\_1 - M\_2||\_{\infty}$$

To prove this, let


with *h*<sup>1</sup> = *N*Γ(*α*+1) .

Similarly, we can prove that

$$\begin{array}{rcl} ||\nu W\_{l\_1} - \nu W\_{l\_2}|| &\leq & h\_2 ||W\_{l\_1} - W\_{l\_2}|| \\ ||\nu W\_{p\_1} - \nu W\_{p\_2}|| &\leq & h\_3 ||W\_{p\_1} - W\_{p\_2}|| \\ ||\nu W\_{f\_1} - \nu W\_{f\_2}|| &\leq & h\_4 ||W\_{f\_1} - W\_{f\_2}|| \\ ||\nu W\_{d\_1} - \nu W\_{d\_2}|| &\leq & h\_5 ||W\_{d\_1} - W\_{d\_2}|| \\ ||\nu I\_{\varepsilon\_1} - \nu I\_{\varepsilon\_2}|| &\leq & h\_6 ||I\_{\varepsilon\_1} - I\_{\varepsilon\_2}|| \\ ||\nu I\_{l\_1} - \nu I\_{l\_2}|| &\leq & h\_7 ||I\_{l\_1} - I\_{l\_2}|| \\ ||\nu I\_{p\_1} - \nu I\_{p\_2}|| &\leq & h\_8 ||I\_{p\_1} - I\_{p\_2}|| \\ ||\nu I\_{f\_1} - \nu I\_{f\_2}|| &\leq & h\_9 ||I\_{f\_1} - I\_{f\_2}|| \\ ||\nu I\_{d\_1} - \nu I\_{d\_2}|| &\leq & h\_{10} ||I\_{d\_1} - I\_{d\_2}||. \end{array}$$

By the definition of Contraction mapping Definition 5, the map *ν* is a contraction map if 0 < *h<sup>i</sup>* < 1 for all *i* = 1, 2, 3, · · · , 10. Therefore, *ν* is a contraction mapping on a compact Banach space *H*. Then by Contraction mapping Theorem 1, *ν* 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*(*M*(*t*)) satisfies the following:

For any *<sup>e</sup>*1,*e*<sup>2</sup> <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* there exists *<sup>S</sup>*<sup>1</sup> <sup>∈</sup> <sup>R</sup>*n*×*<sup>n</sup>* , such that ||*g*(*e*1) − *g*(*e*2)|| ≤ ||*S*1(*e*<sup>1</sup> − *e*2)||.

**Theorem 2.** *Assume that the system* (8) *satisfies the assumption (A1) and the impulsive operator satisfies that*

$$\delta\_{\theta}(M(t\_{\theta})) \quad = \quad - \check{\delta}(M(t\_{\theta}) - M^\*) , \theta = 1, 2, \cdots \text{ } , m\_{\theta}$$

*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 γ*<sup>1</sup> *such that the following inequalities hold:*

$$Q^{\frac{-1}{2}}[\gamma\_1 + \bar{\delta}]^\top Q[\gamma\_1 + \bar{\delta}]Q^{\frac{-1}{2}} \le \quad \gamma\_1 \tag{13}$$

*and*

$$
\boldsymbol{\Omega} \quad = \begin{bmatrix}
\ast & -\boldsymbol{\mathfrak{J}} & \boldsymbol{0} \\
\ast & \ast & -\boldsymbol{\mathfrak{J}} \\
\end{bmatrix} < 0. \tag{14}
$$

**Proof.** Let us consider the system (8) with the initial condition *M*(*t*0) = *M*<sup>0</sup> ∈ Z <sup>+</sup> and an equilibrium point *M*∗ . By using the transformation, N (*t*) = *M*(*t*) − *M*<sup>∗</sup> , then the system (8) is transformed into

$$\begin{aligned} \, \_0^\mathbb{C}D\_t^\mathbb{A} \mathcal{N}(t) &= \, \_0\mathcal{N}(t) + \mathfrak{g}(\mathcal{N}(t)), t \neq t\_\theta, \theta = 1, 2, 3, \dots & \text{in } \Omega\\ \Delta \mathcal{N}(t\_\theta) &= \, \_0\mathcal{N}(t\_\theta^+) - \mathcal{N}(t\_\theta^-) = -\delta \mathcal{N}(t\_\theta), t = t\_\theta, \theta = 1, 2, 3, \dots & \text{in } \Omega\\ \mathcal{N}(t\_0) &= \quad \mathcal{N}\_0 \in \mathbb{Z}^+. \end{aligned} \tag{15}$$

where, N (*t*) = (N1, N2, N3, ..., N10) <sup>&</sup>gt; and *g*¯(N (*t*)) = (*g*¯(N1), *g*¯(N2), .., *g*¯(N10))<sup>&</sup>gt; and N<sup>0</sup> = *M*<sup>0</sup> − *M*<sup>∗</sup> . Let us consider a Lyapunov function as:

*V*(*t*) = N <sup>&</sup>gt;(*t*)*Q*N (*t*), (16)

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

$$\begin{aligned} \, \_0^\mathsf{C} D\_t^\alpha V(t) &\quad \le \, \_2\mathcal{N}^\top(t) Q\_0^\mathsf{C} D\_t^\alpha \mathcal{N}(t) \\ &= \, \_0\mathcal{N}^\top(t) 2\mathcal{Q}[-\mathcal{W}\_1 \mathcal{N}(t) + \mathfrak{g}(\mathcal{N}(t))] \\ &= \, \_0\mathcal{N}^\top(t)(-2\mathcal{Q}\mathcal{W}\_1)\mathcal{N}(t) + \mathcal{N}^\top(t)(2\mathcal{Q})\bar{\mathfrak{g}}(\mathcal{N}(t)) \end{aligned} \tag{17}$$

By Lemma 2,

$$\mathcal{N}^{\top}(t)(\mathcal{Q}\mathcal{Q})\bar{\mathcal{g}}(\mathcal{N}(t)) \quad \leq \ \frac{1}{\xi}\mathcal{N}^{\top}(t)(\mathcal{Q}\mathcal{Q}^{\top})\mathcal{N}(t) + \xi\bar{\mathcal{g}}^{\top}(\mathcal{N}(t))\bar{\mathcal{g}}(\mathcal{N}(t)).\tag{18}$$

By assumption (A1),

$$\begin{split} \left| \mathfrak{g}^{\top}(\mathcal{N}(t))\mathfrak{g}(\mathcal{N}(t)) \right| &= \left| \mathfrak{g}(M(t)) - \mathfrak{g}(M^{\*}) , \mathfrak{g}(M(t)) - \mathfrak{g}(M^{\*}) \right\rangle \\ &= \left| \left( \mathfrak{g}(\mathcal{N}(t) + M^{\*}) \right) - \mathfrak{g}(M^{\*}) , \mathfrak{g}(\mathcal{N}(t) + M^{\*}) - \mathfrak{g}(M^{\*}) \right\rangle \\ &\leq \left| \begin{array}{c} \mathcal{N}^{\top}(t)\mathcal{S}\_{1}^{\top}\mathcal{S}\_{1}\mathcal{N}(t) . \end{array} \right. \tag{19} \end{split} \tag{10}$$

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

$$\begin{split} \, ^C\_0 D\_t^\mathbf{u} V(t) &\leq \quad \mathcal{N}^\top(t)(-2Q\mathcal{W}\_1)\mathcal{N}(t) + \mathcal{N}^\top(t)(\xi^{-1}QQ^\top)\mathcal{N}(t) + \mathcal{N}^\top(t)(\xi \mathbf{S}\_1^\top \mathbf{S}\_1)\mathcal{N}(t) \\ &= \quad \mathcal{N}^\top(t)[-2Q\mathcal{W}\_1 + \xi^{-1}QQ^\top + \xi \mathbf{S}\_1^\top \mathbf{S}\_1]\mathcal{N}(t). \end{split} \tag{20}$$

Let, Ω = −2*QW*<sup>1</sup> + *ξ* <sup>−</sup>1*QQ*<sup>&</sup>gt; + *ξS* > 1 *S*<sup>1</sup> and Ω can be rewritten as

$$
\Omega\_- = \begin{bmatrix}
\* & -\xi & 0 \\
\* & \* & -\xi^{-1}
\end{bmatrix}.
$$

Now, pre and post multiply Ω by diag{*I*, *I*, *ξ*}, we get

$$
\tilde{\Omega} \quad = \begin{bmatrix}
\* & -\tilde{\xi} & 0 \\
\* & \* & -\tilde{\xi}
\end{bmatrix} . \tag{21}
$$

By Schur compliment Lemma 1, Ω˜ < 0. Furthermore, the Equation (20), can be modified as

$$\begin{aligned} \ \_0^C D\_t^\alpha V(t) &\quad \le \ \_0^\top (t) \tilde{\Omega} \mathcal{N}(t) \\ &= \ \_0^\top (t) \mathcal{Q}^{\frac{1}{2}} [-\mathcal{Q}^{-\frac{1}{2}} \tilde{\Omega} \mathcal{Q}^{-\frac{1}{2}}] \mathcal{Q}^{\frac{1}{2}} \mathcal{N}(t) \end{aligned}$$

let, *e*<sup>1</sup> = *λmin*(−*Q* −1 <sup>2</sup> Ω˜ *Q* −1 <sup>2</sup> ) and we know that *V*(*t*) = N <sup>&</sup>gt;(*t*)*Q*N (*t*). This implies that,

$$\begin{array}{rcl} \, \_0^\mathbb{C} D\_t^\alpha V(t) & \le & -\varepsilon\_1 V(t). \end{array} \tag{22}$$

For, *t<sup>θ</sup>* = *t*, *θ* = 1, 2, 3, · · · *m*

*C*

$$\begin{split} V(t\_{\theta}^{+}) &= \mathcal{N}^{\top}(t\_{\theta}^{+}) Q \mathcal{N}(t\_{\theta}^{+}) \\ &= \left[ \mathcal{N}(t\_{\theta}^{-}) + \bar{\delta} \mathcal{N}(t\_{\theta}^{-}) \right]^{\top} Q [\mathcal{N}(t\_{\theta}^{-}) + \bar{\delta} \mathcal{N}(t\_{\theta}^{-})] \\ &= \mathcal{N}^{\top}(t\_{\theta}^{-}) [\gamma\_{1} + \bar{\delta}]^{\top} Q [\gamma\_{1} + \bar{\delta}] \mathcal{N}(t\_{\theta}^{-}) \\ &= \mathcal{N}^{\top}(t\_{\theta}^{-}) Q^{\frac{1}{2}} [\mathcal{Q}^{\frac{-1}{2}} \gamma\_{1} + \bar{\delta}]^{\top} Q [\gamma\_{1} + \bar{\delta} Q^{\frac{-1}{2}}] Q^{\frac{1}{2}} \mathcal{N}(t\_{\theta}^{-}) \\ &\leq \quad \mathcal{N}^{\top}(t\_{\theta}^{-}) Q \mathcal{N}(t\_{\theta}^{-}) = V(\mathcal{N}(t\_{\theta}^{-})) \\ V(t\_{\theta}^{+}) &\leq \quad V(t\_{\theta}^{-}) \end{split} \tag{23}$$

Therefore, we can easily prove that,

$$
\lambda\_{\min}(Q)||M(t)||^2 \quad \le \quad V(t) \le \lambda\_{\max}(Q)||M(t)||^2. \tag{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.


**Table 2.** Data from published literature.

	- 0.9, *Wf*<sup>0</sup> = 0.3, *Wa*<sup>0</sup> = 0.3, *Ie*<sup>0</sup> = 0.9, *Il*<sup>0</sup> = 0.9, *Ip*<sup>0</sup> = 0.9, *If*<sup>0</sup> = 0.3, *Ia*<sup>0</sup> = 0.3, total population *T* = 3000, and the positive scalar used in Theorem 2 as *ξ* = 0.8513 The Figures 4–7 are depicts the dynamics of Equation (5) along with the parameters in Table 2 at various orders of *α* such as *α* = 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.

**Figure 4.** Population dynamics of both WU and WI mosquitoes at *α* = 0.28.

**Figure 5.** Population dynamics of both WU and WI mosquitoes at *α* = 0.68.

**Figure 6.** Population dynamics of both WU and WI mosquitoes at *α* = 0.98.

**Figure 7.** Population dynamics of both WU and WI mosquitoes at *α* = 1.

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 *η*<sup>1</sup> = 0.03, *η*<sup>2</sup> = 0.03, *η*<sup>3</sup> = 0.03, *η*<sup>4</sup> = 0.5 and *η*<sup>1</sup> = 0.5 are fitted.

Figures 8–11 are analyzed the dynamics of the system of Equation (6), with Wolbachia invasion and natural Wolbachia gain at various orders *α* = 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.

**Figure 8.** Population dynamics of Wolbachia invasive model at *α* = 0.28.

**Figure 9.** Population dynamics of Wolbachia invasive model at *α* = 0.68.

**Figure 10.** Population dynamics of Wolbachia invasive model at *α* = 0.98.

**Figure 11.** Population dynamics of Wolbachia invasive model at *α* = 0.28.

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 *δ*<sup>1</sup> = 0.4, *δ*<sup>2</sup> = 0.4, *δ*<sup>3</sup> = 0.3, *δ*<sup>4</sup> = 0.5 and *δ*<sup>5</sup> = 0.5, invasion rates are *η*<sup>1</sup> = 0.03, *η*<sup>2</sup> = 0.03, *η*<sup>3</sup> = 0.03, and gain rates are *η*<sup>4</sup> = 0.5 and *η*<sup>1</sup> = 0.5.

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

**Figure 12.** Population dynamics of Wolbachia invasive model after impulsive control at *α* = 0.28.

**Figure 13.** Population dynamics of Wolbachia invasive model after impulsive control at *α* = 0.68.

**Figure 14.** Population dynamics of Wolbachia invasive model after impulsive control at *α* = 0.98.

**Figure 15.** Population dynamics of Wolbachia invasive model after impulsive control at *α* = 1.

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 nonnegligible 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.

**Author Contributions:** Conceptualization, J.D. and R.R.; methodology, J.D. and R.R.; software, J.D. and R.R.; validation, J.D., R.R., J.A., M.N. and O.B.; formal analysis, J.D. and R.R.; investigation, J.D. and R.R.; resources, J.D. and R.R.; data curation, J.D. and R.R.; writing—original draft preparation, J.D. and R.R.; writing—review and editing, J.D. and R.R.; visualization, J.D., R.R., J.A., M.N. and O.B.; supervision, R.R., J.A., M.N. and O.B.; project administration, J.D. and R.R.; funding acquisition, J.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** J. Alzabut would like to thank Prince Sultan University for supporting and funding this work through research group Nonlinear Analysis Methods in Applied Mathematics (NAMAM) group number RG-DES-2017-01.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This article has been written with the joint partial financial support of SERB-EEQ/2019/000365, the National Science Centre in Poland Grant DEC-2017/25/ B/ST7/02888, RUSA Phase 2.0 Grant No. F 24–51/2014-U, Policy (TN Multi-Gen), Dept.of Edn. Govt. of India, UGC-SAP (DRS-I) Grant No. F.510/8/DRS-I/ 2016(SAP-I), DST-PURSE 2nd Phase programme vide letter No. SR/ PURSE Phase 2/38 (G), DST (FIST - level I) 657876570 Grant No.SR/FIST/MS-I/ 2018/17.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

*Appendix A.1. Wolbachia Infected Mosquitoes Free Equilibrium*

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

$$P\_2 = (\mathcal{W}\_{e\_1 \prime}^\* \mathcal{W}\_{l\_1 \prime}^\* \mathcal{W}\_{p\_1 \prime}^\* \mathcal{W}\_{f\_1 \prime}^\* \mathcal{W}\_{a\_1 \prime}^\* 0, 0, 0, 0, 0)$$

where,

$$\begin{array}{rcl} W\_{\ell\_{1}}^{\*} &=& \frac{T\lambda\_{w\_{f}}\lambda\_{w\_{\mathfrak{z}}}(\lambda\_{w\_{\mathfrak{e}}}+\gamma\_{w\_{\mathfrak{e}}})(\lambda\_{w\_{\mathfrak{l}}}+\gamma\_{w\_{\mathfrak{l}}})^{2}(\lambda\_{w\_{p}}+\gamma\_{w\_{p}})^{2}}{\rho(1-\rho)\Lambda\_{w\_{\mathfrak{e}}}\gamma\_{w\_{p}}^{2}\gamma\_{w\_{\mathfrak{z}}}^{2}\gamma\_{w\_{\mathfrak{l}}}^{2}}\\ W\_{l\_{1}}^{\*} &=& \frac{\gamma\_{w\_{\mathfrak{e}}}}{\lambda\_{w\_{l}}+\gamma\_{w\_{l}}}W\_{\mathfrak{c}\_{1}}^{\*}\\ W\_{p\_{1}}^{\*} &=& \frac{\gamma\_{w\_{l}}\gamma\_{w\_{\mathfrak{e}}}}{(\lambda\_{w\_{l}}+\gamma\_{w\_{l}})(\lambda\_{w\_{p}}+\gamma\_{w\_{p}})}W\_{\mathfrak{c}\_{1}}^{\*}\\ W\_{f\_{1}}^{\*} &=& \frac{\rho\gamma\_{w\_{\mathfrak{z}}}\gamma\_{w\_{\mathfrak{z}}}\gamma\_{w\_{\mathfrak{z}}}}{\lambda\_{w\_{\mathfrak{f}}}(\lambda\_{w\_{\mathfrak{f}}}+\gamma\_{w\_{\mathfrak{f}}})(\lambda\_{w\_{p}}+\gamma\_{w\_{\mathfrak{g}}})}W\_{\mathfrak{c}\_{1}}^{\*}\\ W\_{a\_{1}}^{\*} &=& \frac{(1-\rho)\gamma\_{w\_{\mathfrak{z}}}\gamma\_{w\_{\mathfrak{z}}}}{\lambda\_{w\_{\mathfrak{e}}}(\lambda\_{w\_{\mathfrak{l}}}+\gamma\_{w\_{\mathfrak{l}}})(\lambda\_{w\_{p}}+\gamma\_{w\_{\mathfrak{l}}})}W\_{\mathfrak{c}\_{1}}^{\*} \end{array}$$

These equilibrium points were derived by the following system of equations by putting *I* ∗ *<sup>e</sup>*<sup>1</sup> = 0, *I* ∗ *l*1 = 0, *I* ∗ *<sup>p</sup>*<sup>1</sup> = 0, *I* ∗ *f*1 = 0, *I* ∗ *<sup>a</sup>*<sup>1</sup> = 0. That is,

$$\begin{cases} \frac{\Lambda\_{\overline{w}\epsilon}\mathcal{W}\_{f\_1}^{\ast}\mathcal{W}\_{a\_1}^{\ast}}{T} - \lambda\_{\overline{w}\_{\varepsilon}}\mathcal{W}\_{\varepsilon\_1}^{\ast} - \gamma\_{\overline{w}\_{\varepsilon}}\mathcal{W}\_{\varepsilon\_1}^{\ast} &= 0 \\ \gamma\_{\overline{w}\epsilon}\mathcal{W}\_{\varepsilon\_1}^{\ast} - \lambda\_{\overline{w}\_{\varepsilon}}\mathcal{W}\_{I\_1}^{\ast} - \gamma\_{\overline{w}\_{\mathbb{T}}}\mathcal{W}\_{I\_1}^{\ast} + (1 - \mathfrak{a})\gamma\_{i\_{\varepsilon}}I\_{\varepsilon\_1}^{\ast} &= \mathbf{0} \\ \gamma\_{\overline{w}\iota}\mathcal{W}\_{I\_1}^{\ast} - \lambda\_{\overline{w}\_{\mathbb{T}}}\mathcal{W}\_{p\_1}^{\ast} - \gamma\_{\overline{w}\_{\mathbb{T}}}\mathcal{W}\_{p\_1}^{\ast} + (1 - \beta)\gamma\_{i\_{\mathbb{T}}}I\_{1\_1}^{\ast} &= \mathbf{0} \\ \rho\gamma\_{\overline{w}\_{\mathbb{T}}}\mathcal{W}\_{p\_1}^{\ast} - \lambda\_{\overline{w}\_{\mathbb{T}}}\mathcal{W}\_{f\_1}^{\ast} + (1 - \epsilon)\gamma\_{i\_{\mathbb{T}}}\rho\_{i\_{\mathbb{T}}}I\_{p\_1}^{\ast} &= \mathbf{0} \\ (1 - \rho)\gamma\_{\overline{w}\_{\mathbb{T}}}\mathcal{W}\_{p\_1}^{\ast} - \lambda\_{\overline{w}\_{\mathbb{T}}}\mathcal{W}\_{q\_1}^{\ast} + (1 - \epsilon)\gamma\_{i\_{\mathbb{T}}}(1 - \rho\_{i\_{\mathbb{T}}})I\_{p\_1}^{\ast} &= \mathbf{0} \end{cases}$$

That is,

(i). By solving,

$$
\gamma\_{w\_\varepsilon} \mathcal{W}\_{\varepsilon\_1}^\* - \lambda\_{\overline{w}\_l} \mathcal{W}\_{l\_1}^\* - \gamma\_{\overline{w}\_l} \mathcal{W}\_{l\_1}^\* + (1 - \alpha) \gamma\_{i\_\varepsilon} I\_{\varepsilon\_1}^\* \ = 0
$$

We get the value of *W*∗ *l* as,

$$\mathcal{W}\_l^\* \quad = \begin{array}{c} \gamma\_{w\_\varepsilon} \\ \hline (\lambda\_{w\_l} + \gamma\_{w\_l}) \end{array} \mathcal{W}\_\varepsilon^\* $$

(ii). By solving

$$
\gamma\_{w\_l} \mathcal{W}\_{l\_1}^\* - \lambda\_{w\_p} \mathcal{W}\_{p\_1}^\* - \gamma\_{w\_p} \mathcal{W}\_{p\_1}^\* + (1 - \beta) \gamma\_{l\_l} I\_{l\_1}^\* = 0
$$

We get the value of *W*∗ *<sup>p</sup>* as,

$$\mathbf{W}\_p^\* \quad = \quad \frac{\gamma\_{w\_l}}{\lambda\_{w\_p} + \gamma\_{w\_p}} \mathbf{W}\_l^\*.$$

Substitute the value of *W*∗ *l* from (i),

$$\mathcal{W}\_p^\* \quad = \frac{\gamma\_{w\_l}\gamma\_{w\_\varepsilon}}{(\lambda\_{w\_p} + \gamma\_{w\_p})(\lambda\_{w\_l} + \gamma\_{w\_l})} \mathcal{W}\_\varepsilon^\*$$

(iii). By solving

$$
\rho \gamma\_{w\_p} \mathcal{W}\_{p\_1}^\* - \lambda\_{w\_f} \mathcal{W}\_{f\_1}^\* + (1 - \epsilon) \gamma\_{i\_p} \rho\_{i\_w} I\_{p\_1}^\* = 0
$$

We get the value of *W*∗ *f* as,

$$\begin{array}{rcl} \mathcal{W}\_f^\* &=& \frac{\rho \gamma\_{w\_p}}{\lambda\_{w\_f}} \mathcal{W}\_p^\* \end{array}$$

Substitute the value of *W*∗ *p* from (ii),

$$\mathcal{W}\_f^\* = \begin{array}{c} \rho \gamma\_{w\_p} \gamma\_{w\_\varepsilon} \gamma\_{w\_l} \\ \overline{\lambda\_{w\_f} (\lambda\_{w\_p} + \gamma\_{w\_p}) (\lambda\_{w\_l} + \gamma\_{w\_l})} \mathcal{W}\_\varepsilon^\* \end{array}$$

(iv). By solving

$$(1 - \rho)\gamma\_{w\_p} \mathcal{W}\_{p\_1}^\* - \lambda\_{w\_a} \mathcal{W}\_{a\_1}^\* + (1 - \epsilon)\gamma\_{i\_p}(1 - \rho\_{i\_w})I\_{p\_1}^\* \ = 0$$

We get the value of *W*∗ *<sup>a</sup>* as,

$$\mathcal{W}\_{a}^{\*} \quad = \quad \frac{(1 - \rho)\gamma\_{w\_{p}}}{\lambda\_{w\_{a}}} \mathcal{W}\_{p}^{\*}$$

Substitute the value of *W*∗ *p* from (ii),

$$\mathcal{W}\_a^\* \;= \; \frac{(1 - \rho)\gamma\_{w\_p}\gamma\_{w\_l}\gamma\_{w\_c}}{\lambda\_{w\_a}(\lambda\_{w\_l} + \gamma\_{w\_l})(\lambda\_{w\_p} + \gamma\_{w\_p})} \mathcal{W}\_c^\*$$

(v). By solving

$$\frac{\Lambda\_{\overline{w}\_{\mathfrak{c}}} \mathcal{W}\_{f\_1}^\* \mathcal{W}\_{a\_1}^\*}{T} - \lambda\_{\overline{w}\_{\mathfrak{c}}} \mathcal{W}\_{\mathfrak{c}\_1}^\* - \gamma\_{\overline{w}\_{\mathfrak{c}}} \mathcal{W}\_{\mathfrak{c}\_1}^\* \quad = \mathbf{0}$$

We get the value of *W*∗ *<sup>e</sup>* as,

$$\begin{array}{rcl} W\_{\mathfrak{e}}^{\*} & = & \frac{\Lambda\_{w\_{\mathfrak{e}}}}{T(\lambda\_{w\_{\mathfrak{e}}} + \gamma\_{w\_{\mathfrak{e}}})} W\_{f}^{\*} W\_{a}^{\*} \end{array}$$

Substitute the value of *W*∗ *f* and *W*∗ *a* from (iii) and (iv),

$$\mathcal{W}\_{\varepsilon}^{\*} = \frac{T\lambda\_{w\_f}\lambda\_{w\_\mathfrak{u}}(\lambda\_{w\_\mathfrak{e}} + \gamma\_{w\_\mathfrak{e}})(\lambda\_{w\_l} + \gamma\_{w\_l})^2(\lambda\_{w\_p} + \gamma\_{w\_p})^2}{\rho(1-\rho)\Lambda\_{w\_\mathfrak{e}}\gamma\_{w\_p}^2\gamma\_{w\_\mathfrak{e}}^2\gamma\_{w\_l}^2}$$

#### *Appendix A.2. Wild Mosquitoes Free Equilibrium*

Suppose a successful release of Wolbachia infected mosquitoes replaces the wild mosquitoes by Wolbachia infected mosquitoes. Then the possible equilibrium points can be found by substituting *W*∗ *<sup>e</sup>*<sup>2</sup> = 0, *W*<sup>∗</sup> *l*2 = 0, *W*∗ *<sup>p</sup>*<sup>2</sup> = 0, *W*<sup>∗</sup> *f*2 = 0 and *W*∗ *<sup>a</sup>*<sup>2</sup> = 0 in the following system of equations

$$\begin{array}{rcl} 0 &=& \frac{\Lambda\_{\dot{\iota}\_{\varepsilon}}I\_{f\_{2}}^{\*}\left(W\_{a\_{2}}^{\*}+I\_{a\_{2}}^{\*}\right)}{T} - \lambda\_{\dot{\iota}\_{\varepsilon}}I\_{e\_{2}}^{\*} - \alpha\gamma\_{\dot{\iota}\_{\varepsilon}}I\_{e\_{2}}^{\*} \\ 0 &=& \alpha\gamma\_{\dot{\iota}\_{\varepsilon}}I\_{e\_{2}}^{\*} - \lambda\_{\dot{\iota}\_{\dot{\iota}}}I\_{l\_{2}}^{\*} - \beta\gamma\_{\dot{\iota}\_{\dot{\iota}}}I\_{l\_{2}}^{\*} \\ 0 &=& \beta\gamma\_{\dot{\iota}\_{\dot{\iota}}}I\_{l\_{2}}^{\*} - \lambda\_{\dot{\iota}\_{\dot{\iota}}}I\_{p\_{2}}^{\*} - \varepsilon\gamma\_{\dot{\iota}\_{\dot{\iota}}}I\_{p\_{2}}^{\*} \\ 0 &=& \rho\_{\dot{\iota}}\varepsilon\gamma\_{\dot{\iota}\_{\dot{\iota}}}I\_{p\_{2}}^{\*} - \lambda\_{\dot{\iota}\_{\dot{\iota}}}I\_{f\_{2}}^{\*} \\ 0 &=& (1-\rho\_{\dot{\iota}})\varepsilon\gamma\_{\dot{\iota}\_{\dot{p}}}I\_{p\_{2}}^{\*} - \lambda\_{\dot{\iota}\_{\dot{a}}}I\_{a\_{2}}^{\*}. \end{array}$$

(i) By solving

$$0 \quad = \ (1 - \rho\_i)\epsilon\gamma\_{i\_p}I\_{p\_2}^\* - \lambda\_{i\_a}I\_{a\_2}^\*$$

We get,

$$I\_{a\_2}^\* \quad = \quad \frac{(1 - \rho\_i)\varepsilon \gamma\_{i\_p}}{\lambda\_{i\_a}} I\_{p\_2}^\*$$

(ii) By solving

We get,

*ρieγi<sup>p</sup>*

*I* ∗ *<sup>p</sup>*<sup>2</sup> − *λ<sup>i</sup> f I* ∗ *f*2

$$I\_{f\_2}^\* = -\frac{\rho\_T \gamma\_{n\_p}}{\lambda\_{i\_f}} I\_{p\_2}^\*$$

(iii) By solving

We get,

$$
\beta \gamma\_{i\_l} I\_{l\_2}^\* - \lambda\_{i\_p} I\_{p\_2}^\* - \epsilon \gamma\_{i\_p}
$$

0 = *ρieγi<sup>p</sup>*

$$I\_{l\_2}^\* = -\frac{(\lambda\_{i\_p} + \epsilon \gamma\_{i\_p})}{\beta \gamma\_{i\_l}} I\_{p\_2}^\*$$

*I* ∗ *p*2

(iv) By solving

$$0 \quad = \ \mathfrak{a}\gamma\_{\dot{l}\_{\mathfrak{c}}}I\_{\mathfrak{c}\_2}^\* - \lambda\_{\dot{l}\_{\mathfrak{l}}}I\_{l\_2}^\* - \beta\gamma\_{\dot{l}\_{\mathfrak{l}}}I\_{l\_2}^\*$$

We get,

$$I\_{\mathfrak{e}\_2}^\* \quad = \quad \frac{(\lambda\_{i\_l} + \beta \gamma\_{i\_l})}{\alpha \gamma\_{i\_\varepsilon}} I\_{I\_2}^\*$$

Substitute the value of *I* ∗ *l*2 from (iii),

$$I\_{\mathfrak{e}\_2}^\* = -\frac{(\lambda\_{i\_l} + \beta \gamma\_{i\_l})(\lambda\_{i\_p} + \varepsilon \gamma\_{i\_p})}{\alpha \beta \gamma\_{i\_k} \gamma\_{i\_l}} I\_{p\_2}^\*$$

(v) By solving,

$$0 \quad = \frac{\Lambda\_{\dot{i}\_{\varepsilon}} I\_{f\_{2}}^{\*} (W\_{a\_{2}}^{\*} + I\_{d\_{2}}^{\*})}{T} - \lambda\_{\dot{i}\_{\varepsilon}} I\_{\varepsilon\_{2}}^{\*} - \mathfrak{a} \gamma\_{\dot{i}\_{\varepsilon}} I\_{\varepsilon\_{2}}^{\*}$$

Put *W*∗ *<sup>a</sup>*<sup>2</sup> = 0,

$$\begin{pmatrix} \frac{\Lambda\_{i\_{\ell}}I\_{f\_{2}}^{\*}I\_{a\_{2}}^{\*}}{T} - \lambda\_{i\_{\ell}}I\_{e\_{2}}^{\*} - a\gamma\_{i\_{\ell}}I\_{e\_{2}}^{\*} & = & 0\\ \left(\frac{\Lambda\_{i\_{\ell}}}{T}\right) \left(\frac{\rho\_{i}\varepsilon\gamma\_{i\_{p}}}{\lambda\_{i\_{f}}}I\_{p\_{2}}^{\*}\right) \left(\frac{(1-\rho\_{i})\varepsilon\gamma\_{i\_{p}}}{\lambda\_{i\_{\ell}}}I\_{p\_{2}}^{\*}\right) & = & (\lambda\_{i\_{\ell}} + a\gamma\_{i\_{\ell}})I\_{e\_{2}}^{\*}\\ I\_{p\_{2}}^{\*} & = & \frac{T\lambda\_{i\_{\tilde{p}}}\lambda\_{i\_{\tilde{\ell}}}(\lambda\_{i\_{\tilde{\ell}}} + a\gamma\_{\tilde{\ell}i\_{\tilde{\ell}}})(\lambda\_{i\_{\tilde{\ell}}} + \beta\gamma\_{\tilde{\ell}i\_{\tilde{\ell}}})(\lambda\_{i\_{\tilde{p}}} + \epsilon\gamma\_{i\_{\tilde{\ell}}})}{\Lambda\_{i\_{\ell}}a\beta\rho\_{\ell}(1-\rho\_{i})\epsilon^{2}\gamma\_{i\_{p}}^{2}\gamma\_{i\_{\tilde{\ell}}}\gamma\_{i\_{\tilde{\ell}}}} \end{pmatrix}$$

From (i)–(v) we have the following equilibrium point,

$$P\_3 = (0,0,0,0,0,I\_{\mathfrak{e}\_2 \prime}^\* I\_{l\_2 \prime}^\* I\_{p\_2 \prime}^\* I\_{f\_2 \prime}^\* I\_{a\_2}^\*)$$

where,

$$\begin{aligned} I\_{2\_2}^\* &= \quad \frac{(\lambda\_{i\_l} + \beta \gamma\_{i\_l})(\lambda\_{i\_p} + \epsilon \gamma\_{i\_p})}{\alpha \beta \gamma \gamma\_{i\_r} \gamma\_{i\_l}} I\_{p\_2}^\* \\\ I\_{l\_2}^\* &= \quad \frac{(\lambda\_{i\_p} + \epsilon \gamma \gamma\_{i\_p})}{\beta \gamma\_{i\_l}} I\_{p\_2}^\* \\\ I\_{p\_2}^\* &= \quad \frac{T \lambda\_{i\_f} \lambda\_{i\_a} (\lambda\_{i\_c} + \epsilon \gamma \gamma\_{i\_c})(\lambda\_{i\_l} + \beta \gamma\_{i\_l})(\lambda\_{i\_p} + \epsilon \gamma \gamma\_{i\_p})}{\lambda\_{i\_c} \alpha \beta \rho\_i (1 - \rho\_i) \epsilon^2 \gamma\_{i\_p}^2 \gamma\_{i\_c} \gamma\_{i\_l}} \\\ I\_{f\_2}^\* &= \quad \frac{\rho\_i \epsilon \gamma\_{i\_p}}{\lambda\_{i\_f}} I\_{p\_2}^\* \\\ I\_{q\_2}^\* &= \quad \frac{(1 - \rho\_i) \epsilon \gamma\_{i\_p}}{\lambda\_{i\_q}} I\_{p\_2}^\* \end{aligned}$$

*Appendix A.3. Both Wolbachia and Non-Wolbachia Mosquitoes Co-Existence Equilibrium*

The equilibrium point for the co-existence state can be found by solving the following systems of equations

 Λ*wenW*<sup>∗</sup> *<sup>f</sup>nW*<sup>∗</sup> *an T* − *λweW*<sup>∗</sup> *<sup>e</sup><sup>n</sup>* − *γweW*<sup>∗</sup> *<sup>e</sup><sup>n</sup>* = 0 *γweW*<sup>∗</sup> *<sup>e</sup><sup>n</sup>* − *λwlW*<sup>∗</sup> *ln* − *γwlW*<sup>∗</sup> *ln* + (1 − *α*)*γi<sup>e</sup> I* ∗ *<sup>e</sup><sup>n</sup>* = 0 *γwlW*<sup>∗</sup> *ln* − *λwpW*<sup>∗</sup> *<sup>p</sup><sup>n</sup>* − *γwpW*<sup>∗</sup> *<sup>p</sup><sup>n</sup>* + (1 − *β*)*γ<sup>i</sup> l I* ∗ *ln* = 0 *ργwpW*<sup>∗</sup> *<sup>p</sup><sup>n</sup>* − *λw<sup>f</sup> W*<sup>∗</sup> *fn* + (1 − *e*)*γi<sup>p</sup> ρiw I* ∗ *<sup>p</sup><sup>n</sup>* = 0 (1 − *ρ*)*γwpW*<sup>∗</sup> *<sup>p</sup><sup>n</sup>* − *λwaW*<sup>∗</sup> *<sup>a</sup><sup>n</sup>* + (1 − *e*)*γi<sup>p</sup>* (1 − *ρi<sup>w</sup>* )*I* ∗ *<sup>p</sup><sup>n</sup>* = 0 Λ*i<sup>e</sup> I* ∗ *fn* (*W*∗ *<sup>a</sup><sup>n</sup>* + *I* ∗ *an* ) *T* − *λi<sup>e</sup> I* ∗ *<sup>e</sup><sup>n</sup>* − *αγi<sup>e</sup> I* ∗ *<sup>e</sup><sup>n</sup>* = 0 *αγi<sup>e</sup> I* ∗ *<sup>e</sup><sup>n</sup>* − *λ<sup>i</sup> l I* ∗ *ln* − *βγ<sup>i</sup> l I* ∗ *ln* = 0 *βγ<sup>i</sup> l I* ∗ *ln* − *λi<sup>p</sup> I* ∗ *<sup>p</sup><sup>n</sup>* − *eγi<sup>p</sup> I* ∗ *<sup>p</sup><sup>n</sup>* = 0 *ρieγi<sup>p</sup> I* ∗ *<sup>p</sup><sup>n</sup>* − *λ<sup>i</sup> f I* ∗ *fn* = 0 (1 − *ρi*)*eγi<sup>p</sup> I* ∗ *<sup>p</sup><sup>n</sup>* − *λi<sup>a</sup> I* ∗ *<sup>a</sup><sup>n</sup>* = 0.

$$\text{(i)}$$

$$(1 - \rho\_i)\varepsilon \gamma\_{i\_p} I\_{p\_n}^\* - \lambda\_{i\_d} I\_{a\_n}^\* = \begin{array}{c} 0 \\ I\_{p\_n}^\* \end{array} = \begin{array}{c} 0 \end{array}$$

$$I\_{p\_n}^\* \quad = \begin{array}{c} \lambda\_{i\_d} \\ (1 - \rho\_i)\varepsilon \gamma\_{i\_p} \end{array} I\_{a\_n}^\*$$

(ii)

$$\begin{aligned} \rho\_i \varepsilon \gamma\_{i\_p} I\_{p\_u}^\* - \lambda\_{i\_f} I\_{f\_u}^\* &=& 0\\ I\_{f\_u}^\* &=& \frac{\rho\_i \varepsilon \gamma\_{i\_p}}{\lambda\_{i\_f}} I\_{p\_u}^\*\\ I\_{f\_u}^\* &=& \frac{\rho\_i \lambda\_{i\_d}}{(1 - \rho\_i)\lambda\_{i\_f}} I\_{d\_u}^\* \end{aligned}$$

(iii)

$$\begin{aligned} \beta \gamma\_{\dot{i}\_l} I\_{l\_n}^\* - \lambda\_{\dot{i}\_p} I\_{p\_n}^\* - \varepsilon \gamma\_{\dot{i}\_p} I\_{p\_n}^\* &= \quad 0 \\ I\_{l\_n}^\* &= \quad \frac{\lambda\_{\dot{i}\_d}}{\beta \gamma\_{\dot{i}\_l} (1 - \rho\_i)} \left[ 1 + \frac{\lambda\_{\dot{i}\_p}}{\varepsilon \gamma\_{\dot{i}\_p}} \right] I\_{a\_n}^\* \end{aligned}$$

$$\text{Let } B\_1 = 1 + \frac{\lambda\_{ip}}{\varepsilon \gamma\_{ip}}$$

$$I\_{l\_n}^\* = -\frac{\lambda\_{i\_a} B\_1}{\beta \gamma\_{i\_l} (1 - \rho\_i)} I\_{a\_n}^\*$$

(iv)

(v)

$$\begin{array}{rcl} \alpha \gamma\_{i\_{\varepsilon}} I\_{\varepsilon\_{n}}^{\*} - \lambda\_{i\_{l}} I\_{l\_{n}}^{\*} - \beta \gamma\_{l\_{1}} I\_{l\_{n}}^{\*} &=& 0\\ I\_{\varepsilon\_{n}}^{\*} &=& \frac{(\lambda\_{i\_{l}} + \beta \gamma\_{i\_{l}})}{\alpha \gamma\_{i\_{\varepsilon}}} I\_{l\_{n}}^{\*}\\ I\_{\varepsilon\_{n}}^{\*} &=& \frac{B\_{1} B\_{2} \lambda\_{i\_{l}}}{\alpha \gamma\_{i\_{\varepsilon}} (1 - \rho\_{i})} I\_{a\_{l}}^{\*}\\ \text{Where, } B\_{1} = \left[1 + \frac{\lambda\_{i\_{l}}}{\varepsilon \gamma\_{i\_{l}}}\right]; B\_{2} = \left[1 + \frac{\lambda\_{i\_{l}}}{\beta \gamma\_{i\_{l}}}\right] \end{array}$$

$$\begin{aligned} \frac{\Lambda\_{I\_{n}}I\_{f\_{n}}^{\*}W\_{d\_{n}}^{\*} + \Lambda\_{I\_{n}}I\_{f\_{n}}^{\*}I\_{n}^{\*}}{T} - \lambda\_{i\_{\ell}}I\_{\ell\_{n}}^{\*} - a\gamma\_{i\_{\ell}}I\_{\ell\_{n}}^{\*} &= \quad 0 \\ W\_{d\_{n}}^{\*} &= \quad \frac{T(\lambda\_{i\_{\ell}} + a\gamma\_{i\_{\ell}})}{\Lambda\_{i\_{\ell}}I\_{f\_{n}}^{\*}}I\_{\ell\_{n}}^{\*} - I\_{d\_{n}}^{\*} \\ W\_{d\_{n}}^{\*} &= \quad \frac{TB\_{1}B\_{2}\lambda\_{i\_{\ell}}(\lambda\_{i\_{\ell}} + a\gamma\_{i\_{\ell}})}{a\Lambda\_{i\_{\ell}}\rho\_{i}\gamma\_{i\_{\ell}}} - I\_{d\_{n}}^{\*} \\ W\_{d\_{n}}^{\*} &= \quad \frac{TB\_{1}B\_{2}B\_{3}\lambda\_{i\_{\ell}}}{\rho\_{i}\Lambda\_{i\_{\ell}}} - I\_{d\_{n}}^{\*} \end{aligned}$$

 $\text{Where, } B\_3 = 1 + \frac{\lambda\_{i\_\ell}}{\alpha \gamma\_{i\_\ell}}$  $\text{(vi)}$ 

$$\begin{aligned} (1 - \rho)\gamma\_{w\_p} W\_{p\_u}^\* - \lambda\_{w\_t} W\_{a\_u}^\* + (1 - \varepsilon)\gamma\_{l\_p} (1 - \rho\_{i\_D}) I\_{p\_u}^\* &= 0\\ W\_{p\_u}^\* &= \frac{\lambda\_{w\_d}}{(1 - \rho)\gamma\_{w\_p}} W\_{a\_u}^\* - \frac{(1 - \varepsilon)\gamma\_{l\_p} (1 - \rho\_{i\_D})}{(1 - \rho)\gamma\_{w\_p}} I\_{p\_u}^\*\\ W\_{p\_u}^\* &= \frac{TB\_1 B\_2 B\_3 \lambda\_{i\_f} \lambda\_{w\_d}}{\Lambda\_{i\_\ell} (1 - \rho)\rho\_{i} \gamma\_{w\_p}} - B\_4 I\_{a\_u}^\* \end{aligned}$$

where,  $B\_{4} = 1 + \frac{(1 - \varepsilon)(1 - \rho\_{iw})\lambda\_{ia}}{(1 - \rho)(1 - \rho\_{i})\varepsilon\gamma\_{wp}}$  (viii)

$$\begin{aligned} \rho \gamma\_{w\_p} \mathcal{W}\_{p\_n}^\* - \lambda\_{w\_f} \mathcal{W}\_{f\_n}^\* + (1 - \varepsilon) \gamma\_{i\_p} \rho\_{i\_w} I\_{p\_n}^\* &= 0 \\ \mathcal{W}\_f^\* &= \frac{\rho \gamma\_{w\_p}}{\lambda\_{w\_f}} \left[ \frac{TB\_1 B\_2 B\_3 \lambda\_{i\_f} \lambda\_{w\_d}}{\Lambda\_{i\_l} (1 - \rho) \rho\_i \gamma\_{w\_p}} - B\_4 I\_{a\_l}^\* \right] + \frac{(1 - \varepsilon) \rho\_{i\_w} \lambda\_{i\_u}}{\varepsilon \lambda\_{w\_f} (1 - \rho\_i)} I\_{a\_l}^\* \end{aligned}$$

(viii) *γwlW*<sup>∗</sup> *ln* − *λwpW*<sup>∗</sup> *<sup>p</sup><sup>n</sup>* − *γwpW*<sup>∗</sup> *<sup>p</sup><sup>n</sup>* + (1 − *β*)*γ<sup>i</sup> l I* ∗ *ln* = 0 *W*∗ *ln* = *λw<sup>p</sup>* + *γw<sup>p</sup> γw<sup>l</sup> W*∗ *p<sup>n</sup>* − (1 − *β*)*γ<sup>i</sup> l γw<sup>l</sup> I* ∗ *ln W*∗ *ln* = *λw<sup>p</sup>* + *γw<sup>p</sup> γw<sup>l</sup>* " *TB*1*B*2*B*3*λ<sup>i</sup> f λw<sup>a</sup>* Λ*i<sup>e</sup>* (1 − *ρ*)*ρiγw<sup>p</sup>* − *B*<sup>4</sup> *I* ∗ *an* # − (1 − *β*)*γ<sup>i</sup> l γw<sup>l</sup> λia B*1 *βγ<sup>i</sup> l* (1 − *ρi*) *I* ∗ *an* 

$$\text{(i\*\*)}$$

*γweW*<sup>∗</sup> *<sup>e</sup><sup>n</sup>* − *λwlW*<sup>∗</sup> *ln* − *γwlW*<sup>∗</sup> *ln* + (1 − *α*)*γi<sup>e</sup> I* ∗ *<sup>e</sup><sup>n</sup>* = 0 *W*∗ *<sup>e</sup><sup>n</sup>* = *λw<sup>l</sup>* + *γw<sup>l</sup> γw<sup>e</sup> W*∗ *ln* − (1 − *α*)*γi<sup>e</sup> γw<sup>e</sup> I* ∗ *en W*∗ *<sup>e</sup><sup>n</sup>* = *λw<sup>l</sup>* + *γw<sup>l</sup> γw<sup>e</sup> λw<sup>p</sup>* + *γw<sup>p</sup> γw<sup>l</sup>* " *TB*1*B*2*B*3*λ<sup>i</sup> f λw<sup>a</sup>* Λ*i<sup>e</sup>* (1 − *ρ*)*ρiγw<sup>p</sup>* # − *I* ∗ *an γw<sup>e</sup>* " *B*4(*λw<sup>l</sup>* + *γw<sup>l</sup>* ) *λw<sup>p</sup>* + *γw<sup>p</sup> γw<sup>l</sup>* + *λw<sup>l</sup>* + *γw<sup>l</sup>* (1 − *β*)*γ<sup>i</sup> l γw<sup>l</sup> λi<sup>a</sup> B*1 *βγ<sup>i</sup> l* (1 − *ρi*) + (1 − *α*)*γi<sup>e</sup> λia B*1*B*<sup>2</sup> *αγi<sup>e</sup>* (1 − *ρi*) #

(x) Λ*w<sup>e</sup> W*∗ *<sup>f</sup> W*<sup>∗</sup> *a T* − *λweW*<sup>∗</sup> *<sup>e</sup>* − *γweW*<sup>∗</sup> *<sup>e</sup>* = 0 *W*∗ *<sup>f</sup> W*<sup>∗</sup> *<sup>a</sup>* Λ*w<sup>e</sup> T* = Λ*w<sup>e</sup> ρB*4*γw<sup>p</sup> Tλw<sup>f</sup>* ! (*I* ∗ *a* ) 2 − Λ*w<sup>e</sup> ρB*1*B*2*B*3*λ<sup>i</sup> f ρi*Λ*Ieλw<sup>f</sup>* ! × *λw<sup>a</sup>* (1 − *ρ*) + *B*4*γw<sup>p</sup> I* ∗ *<sup>a</sup>* + Λ*w<sup>e</sup> ρB* 2 1 *B* 2 2 *B* 2 3 *λ* 2 *i f λw<sup>e</sup>* Λ<sup>2</sup> *Ie* (1 − *ρ*)*ρ* 2 *i λw<sup>f</sup>* (*λw<sup>e</sup>* + *γw<sup>e</sup>* )*W*∗ *<sup>e</sup>* = (*λw<sup>e</sup>* + *γw<sup>e</sup>* ) *λw<sup>p</sup>* + *γw<sup>p</sup> λw<sup>l</sup>* + *γw<sup>l</sup> γweγw<sup>l</sup>* − (*λw<sup>e</sup>* + *γw<sup>e</sup>* ) *γw<sup>e</sup>* × " *λw<sup>p</sup>* + *γw<sup>p</sup> λw<sup>l</sup>* + *γw<sup>l</sup> B*4 *γw<sup>l</sup>* + *λw<sup>l</sup>* + *γw<sup>l</sup>* (1 − *β*)*λi<sup>a</sup> B*1 *γw<sup>l</sup> β*(1 − *ρi*) + (1 − *α*)*λi<sup>a</sup> B*1*B*<sup>2</sup> *α*(1 − *ρi*) # *I* ∗ *a*

Λ*w<sup>e</sup> W*∗ *<sup>f</sup> W*<sup>∗</sup> *a T* − *λweW*<sup>∗</sup> *<sup>e</sup>* − *γweW*<sup>∗</sup> *<sup>e</sup>* = 0 Λ*w<sup>e</sup> ρB*4*γw<sup>p</sup> Tλw<sup>f</sup>* ! (*I* ∗ *a* ) <sup>2</sup> <sup>−</sup> (*λw<sup>e</sup>* + *γw<sup>e</sup>* ) *γw<sup>e</sup>* Λ*w<sup>e</sup> ρB*1*B*2*B*3*λ* ∗ *f ρi*Λ*Ieλw<sup>f</sup>* ! *λw<sup>a</sup>* (1 − *ρ*) + *B*4*γw<sup>p</sup>* × " *λw<sup>p</sup>* + *γw<sup>p</sup> λw<sup>l</sup>* + *γw<sup>l</sup> B*4 *γw<sup>l</sup>* + *λw<sup>l</sup>* + *γw<sup>l</sup>* (1 − *β*)*λi<sup>a</sup> B*1 *γw<sup>l</sup> β*(1 − *ρi*) + (1 − *α*)*λi<sup>a</sup> B*1*B*<sup>2</sup> *α*(1 − *ρi*) # *I* ∗ *<sup>a</sup>* + Λ*w<sup>e</sup> ρB* 2 1 *B* 2 2 *B* 2 3 *λ* 2 *i f λw<sup>e</sup>* Λ<sup>2</sup> *Ie* (1 − *ρ*)*ρ* 2 *i λw<sup>f</sup>* = 0

The above equation is a quadratic equation on *I* ∗ *an* . That is,

$$a\_1 I\_{a\_n}^{\*^2} + a\_2 I\_{a\_n}^\* + a\_3 \quad = \quad 0,$$

where,

$$\begin{split} a\_{1} &= \quad \frac{\Lambda \boldsymbol{w}\_{t} \rho \boldsymbol{B}\_{4} \boldsymbol{\gamma} \boldsymbol{w}\_{p}}{T \lambda \boldsymbol{w}\_{f}}; \\ a\_{2} &= \quad \left( \frac{\lambda\_{\boldsymbol{w}\_{t}} + \gamma\_{\boldsymbol{w}\_{t}}}{T \lambda \boldsymbol{w}\_{f}} \right) \left( \frac{\lambda\_{\boldsymbol{w}\_{t}} \lambda\_{i\_{f}} \rho \boldsymbol{B}\_{1} \boldsymbol{B}\_{2} \boldsymbol{B}\_{3}}{\rho\_{i} \lambda\_{i\_{d}} \lambda\_{\boldsymbol{w}\_{f}}} \right) \left( \frac{\lambda\_{\boldsymbol{w}\_{x}}}{(1-\rho)} + \boldsymbol{B}\_{4} \gamma\_{\boldsymbol{w}\_{p}} \right) \\ & \quad \left( \frac{(\lambda\_{\boldsymbol{w}\_{l}} + \gamma\_{\boldsymbol{w}\_{l}})(\lambda\_{\boldsymbol{w}\_{p}} + \gamma\_{\boldsymbol{w}\_{p}}) \boldsymbol{B}\_{4}}{\gamma\_{\boldsymbol{w}\_{l}}} + \frac{(\lambda\_{\boldsymbol{w}\_{l}} + \gamma\_{\boldsymbol{w}\_{l}})(1-\beta) \lambda\_{i\_{d}} \boldsymbol{B}\_{1}}{\gamma\_{\boldsymbol{w}\_{l}} \beta (1-\rho\_{i})} + \frac{(1-\alpha) \lambda\_{i\_{d}} \boldsymbol{B}\_{1} \boldsymbol{B}\_{2}}{\alpha (1-\rho\_{i})} \right); \\ a\_{3} &= \quad \frac{\Lambda\_{\boldsymbol{w}\_{t}} \rho \boldsymbol{T} \boldsymbol{B}\_{2}^{2} \boldsymbol{B}\_{3}^{2} \lambda\_{\boldsymbol{w}\_{t}}}{\Lambda\_{i\_{l}}^{2} (1-\rho) \rho\_{i}^{2} \lambda\_{\boldsymbol{w}\_{f}}}. \end{split}$$

#### These are the equilibrium points presented in Section 4.4.

#### **References**

