1. Introduction

Dengue is a vector-borne disease spread by the female mos-quitoes Aedes aegypti and Aedes albopictus. The mosquitoes Aedes aegypti were originated from Africa but have now been spread in tropical, subtropical and temperate regions of the world. The fast growth in human population, uncontrolled urbanization and inadequate waste management systems have led to abundance of mosquito breeding sites [6]. These are the main causes which bring this global distribution of mosquitoes and consequently spread of the mosquito-borne diseases. Before 1970, dengue was reported from only nine countries, now it has been widespread in the areas of North America, South America, Africa and Southeast Asia [7]. It is important to control dengue disease as more than 50 million people are affected every year [27]. Interestingly, only the female mosquitoes can transmit the virus. The female mosquitoes bite the human as they require blood for reproduction process. On biting the dengue infected person, these female mosquitoes become infected and can transmit the dengue virus to another person [7]. One of the ways to control dengue infection is to control mosquito population. Mosquito population can be chemically controlled by the use of insecticides and / or biological controlled by wolbochia [9,18,26]. The human awareness campaigns can also reduce the spread of infection.

One of the efficient ways to combat the infection is the use of sterile insect technique (SIT). In SIT, sterile male mosquitoes are released near the mosquito breeding sites. Female mosquitoes will not be able to fertilize when mating with these male mosquitoes. In this way, the mosquito population as well as the spread of infection is controlled. The idea of SIT was first conceived by American entomologist, Dr Edward F. Knipling and was successfully implemented to control the spread of screwworm fly in Florida [14]. After that, the SIT technique has been used for many flying insects by various countries [15,16]. To eradicate dengue infection, the government of China and Brazil have also released the sterile male mosquitoes to combat the dengue infection [10,11].

In literature, some mathematical models on control of vector-borne diseases are available [1,2,3,4,9,13,18,19,20,21,23,25]. Few mathematical models have also been formulated to control dengue infection by SIT [3,4,25]. In particular, Esteva and Yang proposed a mathematical model to control the mosquito density by releasing sterile male mosquitoes [3]. Optimal control analysis to control the Aedes aegypti female mosquitoes by SIT strategy has been performed by Thom$\acute{e}$, Yang and Esteva [25]. A pulsed spatial discrete time model using SIT method has also been proposed by Evans and Bishop [4]. They concluded that increased number of sterile will have no further benefit above a threshold. The effectiveness is highly reduced by density dependent mortality of sterile insects. Hendron and Bonsall have developed the $n$ patch dengue model with two control strategies namely vaccination and vector control by SIT method. The dynamics of two serotypes of dengue viruses with effect of cross immunity and antibody-dependent enhancement have been considered without explicitly incorporating the male/female vector (mosquito) dynamics [8].

In this paper, a host-vector model has been proposed to control the primary dengue transmission using SIT strategy. In SIT, the male mosquitoes (Aedes aegypti) are sterilized in the laboratory and are released into the environment to compete with wild male mosquitoes for mating with female mosquitoes. In section 2, a network model for $n$ patches has been formulated. The model analysis of an isolated patch is carried out in section 3. In section 4, the model analysis for $n$ patch network model has been carried out to explore the behavior of system. Conclusions are given in the last section.

---

2. Formulation of host-vector model

Consider A be the class of mosquitoes in aquatic stage (eggs, larvae or pupae). Let $U^{1}$ and $U^{2}$ be the normal male mosquitoes and sterile male mosquitoes population respectively. Let the female mosquitoes are divided in four compartments: $F^{1}$, $F^2, F^3$ and $F^4$. The unmated and uninfected mosquitoes are in $F^1$ compartment. Let $F^2$ be the class of fertilized female mosquitoes and are uninfected. The female mosquitoes $(F^3$) mated with sterile male mosquitoes are unfertilized and will not reproduce. The fertilized female mosquitoes need multiple blood meals to complete the reproduction cycle. The fertilized female mosquitoes who bite infected individuals become infected and transfer to $F^{4}$ class and play crucial role in the spread of dengue infection. The remaining female fertilized mosquitoes bite susceptible/exposed/recovered human population before laying eggs and remain uninfected. No separate class is considered for them. As such both female mosquitoes $F^{2}$ and $F^{4}$ contribute to aquatic stage. Let the susceptible $(S)$, exposed $(E)$, infected $(I)$ and recovered $(R)$ be the host population in the four compartments. Consider the SEIR dynamics in the host population interacting with susceptible female mosquitoes $F^{2}$ and infected female mosquitoes $F^{4}$. The model also incorporates the interstate transition between different classes of male and female mosquitoes. No vertical transmission has been assumed for vector dynamics. Let there be a network of $n$ patches. The host-vector dynamics in each patch is shown in the schematic diagram in Figure 1. Let the subscript $i$ represents the respective variable/parameter in the $i^{th}$ patch; $i$=1, ...$n$. The human migration is considered between patches while no interpatch migration between mosquitoes is considered. The superscripts $S$, $E$, $I$ and $R$ are used to incorporate the differential ability of migration in various compartments (e.g. $m_{ij}^{S}$ represents migration in susceptible compartment). The following $n$ patch metapopulation model is formulated according to the schematic diagram, Figure 1.

$\begin{array}{l} \frac{{d S_i }}{{dt}} &=& \omega_i - \beta^1_iS_iF_i^{4}- \mu_i S_i + \sum\limits^{n}_{j=1,j\neq i} m^{S}_{ij}S_j - \sum\limits^{n}_{j=1,j\neq i} m^{S}_{ji} S_i \label{3.01} \\\end{array}$

(1)

$\begin{array}{l}\frac{{d E_i }}{{dt}} &=& \beta^1_iS_iF_i^{4} - (k_i + \mu_i)E_i + \sum\limits^{n}_{j=1,j\neq i} m^{E}_{ij}E_j - \sum\limits^{n}_{j=1,j\neq i} m^{E}_{ji}E_i \label{3.02}\end{array}$

(2)

$\begin{array}{l}\frac{{d I_i}}{{dt}} &=& k_iE_i - (\alpha_i+\mu_i)I_i + \sum\limits^{n}_{j=1,j\neq i} m^{I}_{ij}I_j - \sum\limits^{n}_{j=1,j\neq i} m^{I}_{ji}I_i \label{3.03} \\\end{array}$

(3)

$\begin{array}{l}\frac{{d R_i}}{{dt}} &=& \alpha_iI_i- \mu_iR_i+\sum\limits^{n}_{j=1,j\neq i} m^{R}_{ij}R_j -\sum\limits^{n}_{j=1,j\neq i} m^{R}_{ji}R_i\label{3.04}\\\end{array}$

(4)

$\begin{array}{l}\frac{{d A_i }}{{dt}} &=& \phi_i(1-\frac{A_i}{C_i})(F^{2}_i+F^{4}_i) - (\gamma_i+d_i)A_i\label{3.05}\\\end{array}$

(5)

$\begin{array}{l}\frac{{d F_i^{1} }}{{dt}} &=& p_i\gamma_iA_i - \frac{\beta^2_iF^{1}_iU^{1}_i}{U^{1}_i+U^{2}_i}-\frac{\beta^2_iF^{1}_iU^{2}_i}{U^{1}_i+U^{2}_i}-d^{1}_iF^{1}_i \label{3.06}\\\end{array}$

(6)

$\begin{array}{l}\frac{{d F_i^{2}}}{{dt}} &=& \frac{\beta^2_iF^{1}_iU^{1}_i}{U^{1}_i+U^{2}_i} -\beta^3_iI_iF^{2}_i-d^{1}_iF^{2}_i \label{3.07}\\\end{array}$

(7)

$\begin{array}{l}\frac{{d F_i^{3}}}{{dt}} &=& \frac{\beta^2_iF^{1}_iU^{2}_i}{U^{1}_i+U^{2}_i} - d^{1}_iF^{3}_i\label{3.08}\\\end{array}$

(8)

$\begin{array}{l}\frac{{d F_i^{4}}}{{dt}} &=&\beta^3_iI_iF^{2}_i- d^{1}_iF^{4}_i \label{3.09} \\\end{array}$

(9)

$\begin{array}{l}\frac{{d U_i^{1}}}{{dt}} &=& (1-p_i)\gamma_iA_i - d^{1}_iU^{1}_i \label{3.100}\\\end{array}$

(10)

$\begin{array}{l}\frac{{d U_i^{2}}}{{dt}} &=& \omega^{1}_i - d^{1}_iU^{2}_i\label{3.101}\end{array}$

(11)

All model parameters are defined in Table 1 and are assumed to be non-negative. The initial conditions associated with the above system are:

$S_i(0)\geq0, E_i(0)\geq 0, I_i(0)\geq 0, R_i(0)\geq 0, 0\leq A_i(0)\leq C_i, F_i^{1}(0)\geq 0, F_i^{2}(0)\geq\\ 0, F_i^{3}(0)\geq 0, F_i^{4}(0)\geq 0, U_i^{1}(0)\geq 0, U_i^{2}(0)\geq 0$

In the next section, the dynamics of the model (1)-(11) is discussed for an isolated patch (i.e $m^{S}_{ij}=m^{E}_{ij}=m^{I}_{ij}=m^{R}_{ij}=0$).

---

3. Model analysis for an isolated patch

In absence of interpatch migration, the patches are decoupled and isolated. The dynamics of each isolated patch is given below where the subscript $i$ is suppressed:

$\begin{array}{l} \frac{{d S }}{{dt}} &=& \omega - \beta^1SF^{4}- \mu S \label{3.1} \\\end{array}$

(12)

$\begin{array}{l}\frac{{d E }}{{dt}} &=& \beta^1SF^{4} - k E - \mu E\label{3.2} \\\end{array}$

(13)

$\begin{array}{l}\frac{{d I}}{{dt}} &=& k E - \alpha I - \mu I \label{3.3} \\\end{array}$

(14)

$\begin{array}{l}\frac{{d R}}{{dt}} &=& \alpha I - \mu R \label{3.4}\\\end{array}$

(15)

$\begin{array}{l}\frac{{d A}}{{dt}} &=& \phi(1-\frac{A}{C})(F^{2}+F^{4}) - (\gamma +d )A\label{3.5}\\\end{array}$

(16)

$\begin{array}{l}\frac{{d F^{1} }}{{dt}} &=& p\gamma A - \frac{\beta^2F^{1}U^{1}}{U^{1}+U^{2}}-\frac{\beta^2F^{1}U^{2}}{U^{1}+U^{2}}-d^{1}F^{1}\label{3.6}\\\end{array}$

(17)

$\begin{array}{l}\frac{{d F^{2}}}{{dt}} &=& \frac{\beta^2F^{1}U^{1}}{U^{1}+U^{2}} - \beta^3IF^{2}- d^{1}F^{2}\label{3.7}\\\end{array}$

(18)

$\begin{array}{l}\frac{{d F^{3}}}{{dt}} &=& \frac{\beta^2F^{1}U^{2}}{U^{1}+U^{2}} - d^{1}F^{3}\label{3.8}\\\end{array}$

(19)

$\begin{array}{l}\frac{{d F^{4}}}{{dt}} &=& \beta^3IF^{2} - d^{1}F^{4}\label{3.9}\\\end{array}$

(20)

$\begin{array}{l}\frac{{d U^{1}}}{{dt}} &=& (1-p)\gamma A - d^{1}U^{1}\label{3.10}\\\end{array}$

(21)

$\begin{array}{l}\frac{{d U^{2}}}{{dt}} &=& \omega^1 - d^{1}U^{2}\label{3.11}\end{array}$

(22)

Let $\mathbf{X}(t)=(S(t), E(t), I(t), R(t), A(t), F^{1}(t), F^{2}(t), F^{3}(t), F^{4}(t), U^{1}(t), U^{2}(t))^\intercal\in\mathbb R_{+}^{11}$. The system (12)-(22) with non-zero initial conditions can be written in the following form:

$\dfrac{{d\mathbf{X}}}{{dt}}=\mathbf{Q}(\mathbf{X}(t),t),\,\,\mathbf{X}(0)=\mathbf{X}_{0}\geq 0\,\, \\ \mathbf{Q}(\mathbf{X}(t),t)=(Q_1(\mathbf{X},t), Q_2(\mathbf{X},t), Q_3(\mathbf{X}, t)....,Q_{11}(\mathbf{X}, t))^\intercal \nonumber$

(23)

---

3.1. Positivity and boundedness of the solution

Proposition 1. The positive cone Int($\mathbb R_{+}^{11}$) is invariant for the system (23).

Proof. Observe that the boundaries of non-negative cone $\mathbb R_{+}^{11}$ are invariant for the system (23) and $\mathbf{Q}(\mathbf{X}(t),t)$ is smooth enough. Applying the existence and uniqueness theorem [22] for differential equations, the system (23) will possess the positive solution.

Proposition 2. Solutions of system (12)-(22) are bounded in the domain $\Omega$= $\{(S, E, I, R, A, F^{1}, F^{2}, F^{3}, F^{4}, U^{1}, U^{2})\in \mathbb R_{+}^{11}\,\,|\,\,\,S+E+I+R\leq\dfrac{\omega}{\mu}$, $A(t)\leq C, F^{1}+ F^{2}+ F^{3}+F^{4}+U^{1}+U^{2}\leq \dfrac{\omega^{1}+\gamma C}{d^{1}}\}$.

Proof. Let $N(t)= S(t) + E(t)+ I(t)+ R(t)$ be the total host population at time $t$. Adding corresponding equations for host dynamics gives

$\dfrac{dN}{dt} = \omega- \mu N \Longrightarrow \limsup\limits_{t\rightarrow\infty}N(t)\leq\dfrac{\omega}{\mu}.$

For vector dynamics, it is clear that $A(t)\leq C \,\,\,\forall\,\,\ t$. Let us prove it by contradiction.

Let $t_0$ be the smallest value of $t$ such that $A(t_0)=C$. Assume $A(t)\leq C$ for $t\in [0, t_0)$ and $A(t)>C$, for $t\in (t_0, T)$ where, $T<\infty$. Accordingly, $\dot{A(t)}>0$ for $t\in (0, t_0)$ and $\dot{A(t)}<0$ for $t\in [t_0, T)$ from (16). By mean value theorem $\dot{A(t)}>0$ for $t\in [t_0, T)$ which is contradiction.

Let $M(t)$=$F^{1}+F^{2}+F^{3}+F^{4}+U^{1}+U^{2}$, now adding equations from (17)-(22) gives

$\begin{array}{c} \dfrac{dM}{dt} = \omega_1+\gamma A(t)-d^{1}(F^{1}(t)+F^{2}(t)+F^{3}(t)+F^{4}(t)+U^{1}(t)+U^{2}(t)) \\ \dfrac{dM}{dt}\leq \omega_1+\gamma C- d^{1}M(t)\\ \Longrightarrow\limsup_{t\rightarrow\infty}M(t)\leq\dfrac{\omega^{1}+\gamma C}{d^{1}}. \end{array}$

Therefore, the system (23) is bounded.

---

3.2. Equilibrium states

The equilibrium states for the system (12)-(22) are obtained by solving

$\begin{array}{c} \mathbf{Q}(\mathbf{X}t),t)=0. \end{array}$

The trivial disease-free state $\breve P_1$ always exists and it is given as

$\begin{array}{l} \breve{S}=\dfrac{\omega}{\mu},\breve{U^2}=\dfrac{\omega_1}{d^{1}}, \breve{E}=0, \breve{I}=0, \breve{R}=0 , \\ \breve{A}=0,\breve{F^1}=0,\breve{F^2}=0,\breve{F^3}=0,\breve{F^4}=0,\breve{U^1}=0. \end{array}$

Note that, this state is without native mosquitoes. Only the susceptible human and sterile male mosquitoes survive.

Let us define the non-dimensional number $T$, the basic offspring number, as

$T =\dfrac{\beta^2\phi p\gamma}{(\beta^2+d^{1})d^{1}(\gamma+d)}.$

(24)

It represents the average number of secondary female mosquitoes produced by single female mosquito, Esteva and Yang [3]. To maintain the mosquito population in environment, $T$ should be greater than 1. The existence of another non-trivial disease-free equilibrium state is established for $T>1$.

Proposition 3. The system admits two non-trivial disease-free states $P_2^+$ and $P_2^-$ when

$T>1\,\,\,\, and \,\,\,\ \omega^{1}<W(=\dfrac{(T-1)^2C\gamma(1-p)}{4T}).$

(25)

Proof. Let us denote $P_2^{\pm}$=$(\bar{S},\bar{E},\bar{I},\bar{R},\bar{A},\bar{F^1},\bar{F^2},\bar{F^3},\bar{F^4},\bar{U^1},\bar{U^2})$. For the disease-free state of the system (12)-(22), substituting $I=0$ gives

$\bar{S}=\dfrac{\omega}{\mu}, \bar{E}=0, \bar{R}=0, \bar{F^4}=0\,\,\, \hbox{and}\,\,\ \bar{U^2}=\dfrac{\omega_1}{d^{1}}.$

(26)

The other state variables $\bar{F^1}, \bar{F^2}, \bar{F^3}$ and $\bar{U^1}$ at equilibrium level are obtained in terms of $\bar A$:

$\bar{F^1}=\dfrac{p\bar A\gamma}{(\beta^2+d^{1})}, \bar{F^2}=\dfrac{\beta^2p\bar A\gamma \bar U^1}{(\bar U^1+\bar U^2)(\beta^2+d^{1})d^{1}},\\ \nonumber \bar{F^3}=\dfrac{\beta^2p\bar A\gamma\bar U^2}{d^{1}(\beta^2+ d^{1})(\bar U^1+\bar U^2)}, \bar{U^1}=\dfrac{\bar A\gamma (1-p)}{d^{1}}$

(27)

$\bar A$ is the root of following quadratic polynomial:

$\begin{array}{c} s(\bar A)=B_1\bar A^{2}+B_2\bar A+B_3=0\\ B_1=\dfrac{\beta^2p\phi\gamma}{(\beta^2+d^{1})(\gamma+d)d^{1}C}, B_2=-\dfrac{\beta^2\phi p\gamma}{(\beta^2+d^{1})d^{1}(\gamma+d)}+1, B_3=\dfrac{\omega_{1}}{\gamma(1-p)} \end{array}$

By writing the quadratic $s(\bar A)$ in terms of $T$ gives

$s(\bar A)=\dfrac{T}{C}\bar A^{2}-(T-1)\bar A+\dfrac{\omega^1}{\gamma (1-p)}=0.$

(28)

Since $p<1$, no positive root of equation (28) is admissible for $T\leq1$. Clearly, the quadratic has real positive roots when the condition (25) is satisfied. The roots are given as

$\bar A^{\pm}=\dfrac{(T-1)C}{2T}(1\pm\sqrt{1-\dfrac{\omega^1}{{W}}});\,\,\ W=\dfrac{(T-1)^2C\gamma(1-p)}{4T}.$

(29)

Accordingly, there exists two non-trivial disease-free equilibrium states $P_2^-$ and $P_2^+$ corresponding to $\bar A^-$ and $\bar A^+$ respectively under condition (25).

These two roots $\bar A^-$ and $\bar A^+$ collapse at $\omega_1=W$. This gives a critical value $T=T^*$ where the quadratic has the unique solution:

$(\bar A=)\bar A^{*}=\dfrac{(T^{*}-1)C}{2T^{*}}$

(30)

$T^*=(1+\dfrac{2\omega^1}{C\gamma(1-p)})[1+\sqrt{(1-\dfrac{1}{(1+\dfrac{2\omega^1}{C\gamma(1-p)})^2})}]>1$

(31)

Therefore, the unique non-trivial disease-free equilibrium point $P^{*}_2$ exists at $T=T^{*}$.

Let us denote the endemic states by $P_3^{\pm}$ $(\hat{S},\hat{E},\hat{I},\hat{R},\hat{A},\hat{F^1},\hat{F^2},\hat{F^3},\hat{F^4},\hat{U^1},\hat{U^2})$. Now, the conditions for the existence of the endemic states $P_3^{\pm}$ are explored. The equilibrium level of state variables are obtained as $\hat{S}=\dfrac{\omega}{(\beta^1\hat{F^4}+\mu)}, \hat{E}=\dfrac{\beta^1\hat S\hat F^4}{(k+\mu)}, \hat{I}=\dfrac{\beta^1k\hat S\hat F^4}{(k+\mu)(\alpha+\mu)}, \hat{R}=\dfrac{\alpha\hat I}{\mu}$, $\hat {F^1}=\dfrac{p\gamma\hat A}{(\beta^2+d^{1})},$ $\hat {F^2}= \dfrac{((\gamma+d)\hat AC\beta^1+\mu\phi(C-\hat A))d^{1}(k+\mu)(\alpha+\mu)}{\phi(C-\hat A)(\beta^1\beta^3k\omega+\beta^1d^{1}(k+\mu)(\alpha+\mu))}$, $\hat{F^3}=\dfrac{\beta^2p\gamma\hat A\omega^1}{d^{1}(\beta^2+d^{1})(\hat{U^1}+\hat{U^2})}$, $\hat{F^4}=\dfrac{\beta^1\beta^3k\omega\hat{F^2}-\mu d^{1}(\alpha+\mu)(k+\mu)}{\beta^1d^{1}(\alpha+\mu)(k+\mu)}$, $\hat{U^1}=\dfrac{(1-p)\gamma\hat{A}}{d^{1}}, \hat{U^2}=\dfrac{\omega^1}{d^{1}}.$

Here, $\hat{A}$ is the root of same quadratic polynomial $s(\hat A)=0$ as given in (28). Accordingly, the condition (25) is necessary for the existence of both the non-trivial disease-free states $P_2^{\pm}$ as well as the endemic states $P_3^{\pm}$. Another condition $\hat{F^4}>0$ is required for existence of $P_3^{\pm}$ states.

Remark 1. It is to be noted that for the non-trivial disease-free states $P_2^{\pm}$ and endemic states $P_3^{\pm}$, the equilibrium densities of following variables are found to be the same:

$\begin{array}{c} \bar A=\hat A, \bar F{^{1}}=\hat F{^{1}}, \bar U{^{1}}=\hat U{^{1}}, \bar U{^{2}}=\hat U{^{2}} \end{array}$

Remark 2. Simplifying the expression for $\hat F^{2}$ gives

$\begin{array}{c} \hat {F_2}=\dfrac{\beta^2\gamma p\hat A\hat {U_1}}{(\hat {U_1}+\hat {U_2})(\beta^2+d^1)(\beta^3\hat I+d^{1})}. \end{array}$

Comparing it with $\bar{F_2}$ yields,

$\bar F^{2}>\hat F^{2}$

(32)

Let us define

$R^{+}_0=\dfrac{\beta^1\beta^3k\omega\bar F^{2+}}{\mu d^{1}(\alpha+\mu)(k+\mu)} \,\,\ \hbox{and}\,\,\ R^{-}_0=\dfrac{\beta^1\beta^3k\omega\bar F^{2-}}{\mu d^{1}(\alpha+\mu)(k+\mu)}.$

(33)

Further, for positive $\hat{F^4}^{+}$ and $\hat{F^4}^{-}$, the following conditions should be satisfied:

$\dfrac{\beta^1\beta^3k\omega\hat F^{2+}}{\mu d^{1}(\alpha+\mu)(k+\mu)}>1\,\,\ \hbox{and}\,\,\ \dfrac{\beta^1\beta^3k\omega\hat F^{2-}}{\mu d^{1}(\alpha+\mu)(k+\mu)}>1$

(34)

Using (32) gives the conditions for the existence of the endemic states $P_3^{\pm}$.

$R^{+}_0>1\,\,\ \hbox{and}\,\,\ R^{-}_0>1$

(35)

Remark 3. It can be seen from (29) that $\bar A^{+}>\bar A^{-}$. Similarly, it can be easily proved that $\bar F^{2+}>\bar F^{2-}$. Accordingly, it can be concluded that:

$\begin{array}{c} R^{+}_0>R^{-}_0 \end{array}$

Let us define $R_0$ as

$R_0=\max (R^{+}_0, R^{-}_0)(>1).$

(36)

Accordingly, the following proposition is established for the existence of endemic states $P_3^{\pm}$:

Proposition 4. When condition (25) is satisfied, the two positive endemic equilibrium states $P_3^{-}$ and $P_3^{+}$ will exist provided the equation (36) is satisfied.

Consider the following choice of parameters [3]:

$\begin{array}{c} d=0.05, d^{1}=0.0714, p=0.5, \gamma=0.075 \end{array}$

The existence of $P_2^{\pm}$ and $P_3^{\pm}$ states depend on the $T$ and $W$. The parameters $\beta^2$ and $\phi$ are considered to be critical for their existence. The curve $T=1$ in Figure 2 divides the $\beta^2-\phi$ plane into two regions. On the right of this curve ($T>1$) both the equilibrium states $P_2^{\pm}$ and $P_3^{\pm}$ may exist depending on the other conditions. In particular, the states $P_2^{\pm}$ and $P_3^{\pm}$ will not exist on the left of the curve $T=1$ $(T<1)$.

---

3.3. The basic reproduction number

The basic reproduction number has been computed by next generation approach [12]. It is defined as the average number of secondary infections produced by single infected individual in susceptible population. Using the Remark 3, it is given as

$R_0=\dfrac{\beta^1\beta^3k\omega \bar F^{2+} }{(\alpha + \mu)d^{1}\mu(\mu + k)}.$

(37)

---

3.4. Stability of trivial disease-free state ($P_1$)

For local stability, the Jacobian matrix $J[P]$ of the system (12)-(22) about a state $P$ is given as:

$\tiny \begin{bmatrix} -\beta^1F^4-\mu&0&0&0&0&0&0&0&-\beta^1S&0&0 \\ \beta^1F^4& -k-\mu&0&0&0&0&0&0&\beta^1S&0&0 \\ 0&k&-\alpha-\mu&0&0&0&0&0&0&0&0\\ 0&0&\alpha&-\mu&0&0&0&0&0&0&0 \\ 0&0&0&0&j_{5,5}&0&j_{5,7}&0&j_{5,9} &0&0\\ 0&0&0&0&p\gamma &-\beta^2-d^{1}&0&0&0&0&0\\ 0&0&-\beta^3F^{2}&0&0 &j_{7,6} &-d^{1}-\beta^3I& 0&0&j_{7,10}&j_{7,11}\\ 0&0&0&0&0&j_{8,6}& 0&-d^{1}&0&j_{8,10}& j_{8,11}\\ 0&0&\beta^3F^{2}&0&0&0&\beta^3I&0 &-d^{1}&0&0\\ 0&0&0&0&(1-p)\gamma&0&0&0&0 &-d^{1}&0\\ 0&0&0&0&0&0&0&0&0&0 &-d^{1}\\ \end{bmatrix}$

$\begin{array}{c} j_{5,5}=-\gamma-d+\frac{-\phi}{C}(F^2+F^4), j_{7,6}=\frac{\beta^2U^1}{U^1+U^2}, j_{7,10}=\frac{\beta^2F^{1}U^2}{(U^1+U^2)^2}, \\ j_{7,11}= \frac{-\beta^2 F^{1}U^1}{(U^{1}+U^2)^2}, j_{8,6}=\frac{\beta^2U^2}{(U^{1}+U^2)}, j_{8,10}=-j_{7,10}, j_{8,11}=-j_{7,11} \\ j_{5,7}=\phi(1-\frac{A}{C})=j_{5,9} \end{array}$

The eleven eigenvalues of the Jacobian matrix about the trivial disease-free state $P_1$ are obtained as

$\begin{array}{c} -\alpha-\mu, -\mu\text{(multiplicity 2)}, -\beta^2-d^{1}, -d^{1}\text{(multiplicity 5)}, -\mu-k, -\gamma-d^{1}. \end{array}$

Since all the eigenvalues have negative real part, the state $P_1$ will always be locally asymptotically stable.

Proposition 5. The locally asymptotically stable trivial disease-free state $P_1$ is also globally stable for

$T<1.$

(38)

Proof. Consider the positive definite function $L(A, F^{1}, F^{2}, F^{4})$ for arbitrarily chosen positive constants $C_1, C_2$, $C_3$ and $C_4$:

$\begin{array}{l} L(A, F^{1}, F^{2}, F^{3}, F^{4}) = C_1A+ C_2F^{1} + C_3F^{2} + C_4F^{4} \end{array}$

Computing its time derivative $\dot L(A, F^{1}, F^{2}, F^{4})$ along the trajectories of system (12)-(22),

$\begin{aligned} \dot L(A, F^{1}, F^{2}, F^{4}) =& C_1\dot A+ C_2\dot F^{1}+ C_3\dot F^{2}+ C_4\dot F^{4} \\ \dot L(A, F^{1}, F^{2}, F^{4}) =& C_1(\phi(1-\dfrac{A}{C})(F^{2}+ F^{4})-(\gamma +d)A)+ C_2(p\gamma A-\beta^2F^{1}-d^{1}F^{1})\\ +& C_3(\dfrac{\beta^2F^{1}U_1}{(U_1+U_2)}-\beta^3IF^{2}-d^{1}F^{2})+ C_4(\beta^3IF^{2}-d^{1}F^{4})\\ \dot L(A, F^{1}, F^{2}, F^{4}) \leq& -F^{2}(C_3d^{1}-C_1\phi)-F^{4}(C_4d^{1}-C_1\phi)- A(C_1(\gamma+d)-C_2p\gamma)\\ &- F^{1}(C_2(\beta^2+d^{1})-C_3\beta^2) \end{aligned}$

For the derivative of the function $L(A, F^{1}, F^{2}, F^{4})$ to be negative definite

$\begin{array}{c} {d^{1}}C_{3}>{\phi}C_{1}; d^{1}C_{4}>{\phi}C_{1}; (\gamma+d)C_{1}>{p\gamma}C_{2}; (\beta^2+d^{1})C_{2}>\beta^2C_{3} \end{array}$

or

$\begin{array}{c} \dfrac{C_{3}}{C_{1}}>\dfrac{\phi}{d^{1}}, \dfrac{C_{4}}{C_{1}}>\dfrac{\phi}{d^{1}}, \dfrac{C_{1}}{C_{2}}>\dfrac{p\gamma}{(\gamma+d)}, \dfrac{C_{2}}{C_{3}}>\dfrac{\beta_2}{(\beta^2+d^{1})} \end{array}$

Now, choosing $C_3$=$C_4$=1 and using above inequalities give

$\begin{array}{c} \dfrac{\beta^2\phi p\gamma}{(\beta^2+d^{1})d^{1}(\gamma+d)}(=T)<1 . \end{array}$

Accordingly, the function $L(A, F^{1}, F^{2}, F^{4})$ is a Lyapunov function for the condition (38). Since $\{P_1\}$ is the largest invariant set that contains the subset in which $\dot L=0$ for $A=0, F^{1}=0, F^{2}=0, F^{4}=0$. By applying LaSalle's invariance principle [17], locally asymptotically stable disease-free state $P_1$ is also globally asymptotically stable under (38).

---

3.5. Local stability of non-trivial disease-free states ($P_2^{\pm}$)

Proposition 6. For the local stability of non-trivial disease-free states $P_2^{\pm}$, the following conditions should be satisfied:

$R{_0}(=\dfrac{\beta^1\beta^3\bar {F_2}^{+} k\omega}{(\alpha + \mu)\mu d^{1}(\mu + k)})<1$

(39)

and

$k(\bar A)(=D_4)=\dfrac{T\bar A^2}{C}-\dfrac{\omega^1}{\gamma(1-p)}>0$

(40)

Proof. For the local stability of $P_2^{\pm}$ states, the 11 eigenvalues of Jacobian matrix of the system (12)-(22) about the states $P_2^{\pm}$ are given as: $-\mu$ (multiplicity 2), -$d^{1}$ (multiplicity 2) and the remaining eigenvalues are the roots of cubic polynomial ($q_3(\lambda)$) and fourth degree polynomial ($q_4(\lambda)$), given as

$\begin{array}{c} q_3(\lambda)=\lambda^3+ B_1\lambda^2+ B_2\lambda + B_3=0 \end{array}$

$\begin{aligned} B_1=&\alpha + 2\mu + d^{1} + k\\ B_2=&\mu(\mu + 2d^{1}) + (\mu + d^{1})k + \alpha(\mu + d^{1} + k)\\ B_3=&(\alpha + \mu)d^{1}(\mu + k) - \beta^1\beta^3\bar F_2 k \bar S \end{aligned}$

The fourth degree polynomial is given as

$\begin{array}{c} q_4(\lambda)=\lambda^4+D_1\lambda^3+ D_2\lambda^2 +D_3\lambda+ D_4 =0. \end{array}$

$\begin{aligned} D_1=&\beta^2 + d + 3d^{1} + \dfrac{\bar F^2\phi}{C} + \gamma >0\\ D_2=&\dfrac{(\bar F_2\phi+ C(d + d^{1} + \gamma))(\beta^2+3d^{1})}{C}+\beta^2d^{1}>0\\ D_3=&d^{1}\Big((\gamma+d)(\beta^2+2d^{1})+d^{1}(\beta^2+d^{1})\\+&\dfrac{(\gamma+d)(2\beta^2+3d^{1})}{d^{1}(\bar U^1+\bar U^2)}\Big(\dfrac{T\bar A^2}{C}-\dfrac{\omega^1d^{1}}{\gamma(1-p)(2\beta^2+3d^{1})}\Big)\Big)\\ D_4=&\dfrac{T\bar A^2}{C}-\dfrac{\omega^1}{\gamma(1-p)} \end{aligned}$

By Routh-Hurwitz criteria, all the roots of $q_3(\lambda)$ will be negative provided

$\begin{array}{c} R{_0}(=\dfrac{\beta^1\beta^3\bar {F^2}^{+} k\omega}{(\alpha + \mu)\mu d^{1}(\mu + k)})<1. \end{array}$

Further, all the roots of fourth degree polynomial $q_4(\lambda)$ will have negative real part iff

$D_i>0, i=1...4$

(41)

$D_1D_2D_3>D_3^2+D_1^2D_4.$

(42)

Note that $D_1$ and $D_2$ are positive and also $D_4$ is positive provided the condition (40) is satisfied. This condition also ensures the positivity of $D_3$. The condition (42) is simplified and found to be satisfied for $D_4>0$.

All eigenvalues have negative real part provided the conditions (39) and (40) are satisfied simultaneously and this completes the proof.

Corollary 1. The non-trivial disease-free state $P_2^{-}$ remains unstable always while the state $P_2^{+}$ is found to be locally asymptotically stable under the condition (39).

Proof. For the stability of the two non-trivial disease-free states $P_2^{\pm}$ explicitly namely $P_2^-$ and $P_2^+$, the sign of $k(A^-)$ and $k(A^+)$ are critical. The $k(\bar A)$ is defined in (40). For the point $A^*$ where $A^+$ and $A^-$ collapse, we have $A^-\leq A^*\leq A^+$. Using the expression for $A^*$ as given in (30), $k(A^*)$ is found to be zero. Accordingly, $k(A^-)$ will be of negative sign giving instability of $P_2^-$ and $k(A^+)$ will be of positive sign. Hence, $P_2^{-}$ is always unstable while $P_2^{+}$ will be locally asymptotically stable for (39) to be satisfied.

Remark 4. The condition (40) is always satisfied for $P_2^+$. Therefore its stability depends on the condition (39).

Remark 5. When condition (39) is not satisfied, the non-trivial disease-free state $P_2^+$ will also become unstable. Consequently, the endemic states $P_3^+$ and $P_3^-$ will start to exist by (36).

The Figure 3 is two parameter bifurcation diagram with respect to $\beta^1$ and $\beta^3$ showing the region of stability of the states $P_2^{\pm}$ as well as the existence of the $P_3^{\pm}$ states. Consider the following choice of parameters with varying transmission rates $\beta^1$ and $\beta^3$:

$\begin{array}{c} \omega=0.039, \beta^2=0.5, k=0.1667, \phi=0.5, \mu=0.00004, d=0.1096, d^{1}=0.0714, \alpha=\\0.7, p=0.5, \omega^1=0.97, \gamma=0.7, C=500 \end{array}$

In Figure 3, the curve $R_0=1$ bifurcates the parameter plane $\beta^1-\beta^3$ into two regions. The non-trivial disease-free state $P_2^+$ is locally stable in the region $R_0<1$. The existence of endemic states $P_3^+$ and $P_3^-$ is possible in the region $R_0>1$.

---

3.6. Numerical simulations

The local stability of $P_3^{\pm}$ states could not be achieved analytically due to large and complex expressions. To discuss the local stability of $P_3^{\pm}$ states, the numerical simulations have been performed for the data given in Table 2.

Further, the ratio between the carrying capacity to the net sterile male mosquitoes i.e. ($\frac{C}{\omega^1/d^{1}}$) is assumed to be 0.8 [3]. It is found that for the data set given in Table 2, the basic offspring number $T(=19.06388689)>1$ and the threshold $R_0$ =8.179178686($>1$). Accordingly, the states $P_3{^{+}}$(0.76383733, 0.01181548, 0.00656461, 50.31770846,394.752323, 19.19005978, 46.75477253,141.25412, 0.128961,207.32790,624.64985) and $P_3{^{-}}$(15.250814, 0.008414931, 0.00467529, 35.836020, 79.020, 3.8413960, 2.3417136, 35.314431, 0.00460008, 41.502141, 624.64985) exist by condition (36). For the stability of $P^-_3$ state, numerical simulations have been performed for the four initial conditions $Y_1$, $Y_2$, $Y_3$ and $Y_4$ (defined in caption) in the neighborhood of the state $P^-_3$.

It is observed that the trajectories starting from the initial conditions $Y_1$ and $Y_2$ converge to the point $Y$ which is the projection of the state $P^+_3$ on the hyperplane. Moreover, it is checked that all other state variables also tend to corresponding values of $P^+_3$. Thus, the solution trajectories converge to $P^+_3$. Similarly, the trajectories starting from the initial conditions $Y_3$ and $Y_4$ converge to the state $P_1$. It shows the unstable behavior of $P^-_3$ state. The phase plot in 3D hyperplane $S-I-F^{4}$ has been drawn in Figure 4a. Further, for the stability of $P^+_3$ state, numerical simulations have been carried out for the four initial conditions (defined in caption) in the neighborhood of the state $P^+_3$.

It is found that starting with the initial conditions $Z_1$, $Z_2$, $Z_3$ and $Z_4$, the solution converges to the $P^+_3$ state. The phase plot in 3D hyperplane $S-I-F^{4}$ has been drawn in Figure 4b.

A bifurcation diagram has been drawn for the system (12)-(22) with respect to sterile male mosquitoes rate $\omega^{1}$ in Figure 5. The transmission rates from vector to host($\beta_1$) and host to vector ($\beta_3$) are considered to be 0.002 nd 0.003 and rest of the parameter values are taken from Table 2. It is observed that the non-trivial disease-free states $P^+_2$ and $P^-_2$ exist from $\omega^{1}$=0 to $W$=80.0699 and collapse in $P_2^*$ at $\omega^{1}$=$W$(=80.0699). The expression for $W$ is given in (29). For $W>80.0699$, only the state $P_1$ exists and stable. Further, it is found that at $\omega^{1}=\omega^{**1}$ =30, $R_0$= 1.002868566. The value of $\omega^{**1}$ can be obtained using the condition (37). Therefore, for $\omega^{1}>30$, the non-trivial disease-free state $P^+_2$ is found to be stable. However, for $\omega^{1}<30$, the endemic states $P^+_3$ and $P^-_3$ start to exist while the disease-free state $P^+_2$ becomes unstable. Thus, for $\omega^{1}<W$, the control on natural mosquitoes will depend on the initial conditions and beyond this $W$, the control will work successfully irrespective of any initial condition.

---

3.7. Discussion

In this paper, a host-vector model has been formulated to analyze the effect of SIT control. It is observed from the analysis that when basic offspring number $(T)$ is less than one then all the natural mosquitoes will be eliminated and only the sterile mosquitoes will survive in the long run. On the other hand if $T$ is greater than $T^*$ (given in (31)), the two non-trivial disease-free states start to exist in which bigger one is found to be stable for $R_0<1$ while the lower one is unstable always. At $T=T^*$, the non-trivial disease-free states collapse. However, for $R_0>1$, two endemic states exist. Numerical simulations have been carried out to analyze the stability of endemic states. It is found that one state gets stable while other remains unstable for choice of relevant data from literature. To control the disease, a threshold has been computed for the rate $\omega^{1}$ at which the sterile mosquitoes be introduced. If the initial population size is in neighborhood of disease-free state then SIT control will work effectively and disease will die out for ${\omega^{**}}^{1}<\omega^{1}<W$. However, when $\omega^{1}<{\omega^{**}}^{1}$, the disease will persist, as after this value endemic states start to exist (shown in Figure 5).

To study the effect of human migration in disease transmission, the model (1)-(11) has been analyzed in next section.

---

4. $n$ patch network model analysis

---

4.1. Basic reproduction number for $n$ patch model

When no movement of human is considered, $R_0$ is given as

$\begin{array}{c} R_0= \max{(R_{0i})}; i=1, 2, 3...n\\ \text{where}, R_{0i}=\dfrac{k_i \beta_i^{1}\beta_i^{3}\omega_iF^{20}_i}{(\mu_i+\alpha_i)(k_i+\mu_i)\mu_id^{1}_i} \end{array}$

When the patches are coupled due to movement of human population, the basic reproduction number for the network model (1)-(11) is computed below:

Let us denote the susceptible human population and fertilized female population at disease free state for the $i^{th}$ node be $S^{0}_i$ and $F^{20}_i$ and define

$\begin{array}{c} q_i=\beta_i^{1}S^{0}_i\text{ and }l_i=\beta_i^{3}F^{20}_i \end{array}$

Further, define $n\times n$ diagonal matrices $D_q$, $D_l$, $D_k$ and $D_d^{1}$ as

$D_q$= diag$(q_1, q_2, ....q_n)$, $D_l$= diag$(l_1, l_2, ....l_n)$, $D_k$= diag$(k_1, k_2, ....k_n)$ and $D_d^{1}$= diag$(d^1_1, d^1_2, ....d^1_n)$

The matrices $Q_E$ and $Q_I$ are also defined as

$Q_E = \begin{bmatrix} \mu_1+k_1+\sum_{j=1}^{1=n} m^{E}_{j1}&-m^{E}_{12}&...& -m^{E}_{1n} \\ -m^{E}_{21}&\mu_2+k_2+\sum_{j=1}^{1=n} m^{E}_{j2}&...& -m^{E}_{2n} \\ ...&......&....&....&\\ -m^{E}_{n1} &.....& ..&\mu_n+k_n+\sum_{j=1}^{1=n} m^{E}_{jn}\\ \end{bmatrix}$

and

$Q_I = \begin{bmatrix} \mu_1+\alpha_1+\sum_{j=1}^{1=n} m^{I}_{j1}&-m^{I}_{12}&...& -m^{I}_{1n} \\ -m^{I}_{21}&\mu_2+\alpha_2+\sum_{j=1}^{1=n} m^{I}_{j2}&...& -m^{I}_{2n} \\ ...&......&.....&.....&\\ -m^{I}_{n1} &...&...&\mu_n+\alpha_n+\sum_{j=1}^{1=n} m^{I}_{jn}\\ \end{bmatrix}$

The basic reproduction number is computed by next generation approach [12]. The Jacobian matrices of the system (1)-(11) for the new infections $(F)$ and transfer from one compartment to another $(Y)$ are given below:

$F = \begin{bmatrix} 0&0&D_q \\ 0&0&0 \\ 0&D_l&0 \\ \end{bmatrix}$

and

$Y = \begin{bmatrix} Q_E&0&0\\ -D_k&Q_I&0 \\ 0&0&D_{d^{1}}\\ \end{bmatrix}$

$FY^{-1} = \begin{bmatrix} 0&0&D_qD^{-1}_{d^{1}}\\ 0&0&0 \\ Q^{-1}_{E}Q^{-1}_{I}D_kD_l&D_lQ^{-1}_{E}&0\\ \end{bmatrix}$

The dominant eigenvalue of the next generation matrix $FY^{-1}$ is the basic reproduction number for $n$ patch model, $R_{0}^n$(say)

$R_0^{n}= \rho(Q^{-1}_E Q^{-1}_I D_kD_lD_qD^{-1}_{d^{1}}).$

(43)

Particularly, for single node $(n=1)$, $R_0$ is obtained as

$\begin{array}{c} R_0= \dfrac{k_1\beta_1^{1}\beta_1^{3}S^{0}_1F^{20}_1 }{(\mu+\alpha_1)(k_1+\mu)d^{1}}. \end{array}$

---

4.2. Numerical examples

Let us perform the numerical experiments to control the disease in a network using SIT. For $n=2$ and $n=3$, the network topologies are given in Figure 6 and Figure 7 respectively. Let us assume that $m^{S}_{ij}= m^{E}_{ij}= m^{I}_{ij}= m^{R}_{ij}= m_{ij}$ and $\omega_i^1$ and $m_{ij}$ are varying in different patches. For simplicity, rest of the data is considered to be the same in all patches given in Table 3.

$\textbf{Network with }n=\textbf{2}$

Observe that $R_0=3.6819 (>1)$ in all the isolated patches for the data set given in Table 3. Accordingly, the disease is endemic in patches. Let us apply SIT in patch-1 so that the disease is controlled in it. The role of network coupling in controlling the disease in the network is explored for different combinations of migration parameters $m_{ij}$ in Table 4. Further, if SIT control is reduced in patch-1 i.e. $\omega^{1}_1=60$ and no control is applied in patch-2 then for different combinations of migration parameters, the cases have been discussed in Table 5. It may be noted from Table 4, case (a) that the disease is controlled in the first patch where SIT is applied. The migration from second patch may bring infection in first patch. If the adequate number of sterile mosquitoes are present in the first patch, this additional infection will be eliminated. Consequently, the network is disease-free and SIT is successful. When $\omega^{1}_1$ is reduced (Table 5), the number of sterile mosquitoes present in the first patch are not sufficient enough to control the disease. Accordingly, the disease persists in the network (case (e)). However, changing the coupling (migration) suitably may again control the disease in the network (case (f)).

$\mathbf{Network with n=3}$

In this case, the following two cases have been considered in Table 6: From the Table 6, it is observed that the basic reproduction number for case (g) is less than one. The disease may be controlled in 3 node network model for this case. By numerically solving the network model (1)-(11) for $n$=3, it is found that the infection level converges to zero as shown in Figure 8e. Further, for case (h), $(R_0^n)>1$. The time series for infective population ($I_1$, $I_2$ and $I_3$) have been drawn in Figure 8f. It is evident that the infection persists. Thus, disease is controlled in case(g), but there is a failure of SIT in case (h).

Keeping the results of above numerical experiments, the following conclusions can be drawn:

It is not necessary to apply SIT in the whole network. The disease can be controlled in network by applying SIT in one patch only. This is possible with suitable coupling of patches. Further, this selective way of applying SIT is more economical (cost effective). The success of SIT depends on the coupling strength of the network (migration parameters $m_{ij}$) and the recruitment rate of sterile mosquitoes $\omega^{1}_1$. Accordingly, they are suggested to be critical parameters.

---

5. Conclusion

In this paper, $n$ patch network model has been formulated to control disease transmission by using SIT. First, the dynamics of an isolated patch has been analyzed. The basic reproduction number has been computed. For the existence and stability of disease-free and endemic states, the two critical parameters namely basic offspring number $(T)$ and basic reproduction number ($R_0$) have been identified. The bifurcation diagram has been plotted to show the existence and stability regions of disease-free and endemic states for an isolated patch. The critical level of sterile male mosquitoes has been obtained to control the disease. For $n$ patch model, the basic reproduction number has been computed. The numerical simulations have been performed by considering two and three nodes with different combinations of migration and SIT recruitment to study the effects in disease elimination and persistence. It is concluded that the cost effective way to control the disease in suitably coupled network is to apply SIT in one patch only.

---

Acknowledgments

We would like to thank IFCAM for visiting support for this research collaboration. Authors also would like to thank anonymous reviewers for their valuable comments. Their comments have enhanced the clarity and the quality of results and text.

---