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ˊ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 U1 and U2 be the normal male mosquitoes and sterile male mosquitoes population respectively. Let the female mosquitoes are divided in four compartments: F1, F2,F3 and F4. The unmated and uninfected mosquitoes are in F1 compartment. Let F2 be the class of fertilized female mosquitoes and are uninfected. The female mosquitoes (F3) 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 F4 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 F2 and F4 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 F2 and infected female mosquitoes F4. 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 ith 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. mSij represents migration in susceptible compartment). The following n patch metapopulation model is formulated according to the schematic diagram, Figure 1.
dSidt=ωi−β1iSiF4i−μiSi+n∑j=1,j≠imSijSj−n∑j=1,j≠imSjiSi
|
(1)
|
dEidt=β1iSiF4i−(ki+μi)Ei+n∑j=1,j≠imEijEj−n∑j=1,j≠imEjiEi
|
(2)
|
dIidt=kiEi−(αi+μi)Ii+n∑j=1,j≠imIijIj−n∑j=1,j≠imIjiIi
|
(3)
|
dRidt=αiIi−μiRi+n∑j=1,j≠imRijRj−n∑j=1,j≠imRjiRi
|
(4)
|
dAidt=ϕi(1−AiCi)(F2i+F4i)−(γi+di)Ai
|
(5)
|
dF1idt=piγiAi−β2iF1iU1iU1i+U2i−β2iF1iU2iU1i+U2i−d1iF1i
|
(6)
|
dF2idt=β2iF1iU1iU1i+U2i−β3iIiF2i−d1iF2i
|
(7)
|
dF3idt=β2iF1iU2iU1i+U2i−d1iF3i
|
(8)
|
dF4idt=β3iIiF2i−d1iF4i
|
(9)
|
dU1idt=(1−pi)γiAi−d1iU1i
|
(10)
|
All model parameters are defined in Table 1 and are assumed to be non-negative. The initial conditions associated with the above system are:
Si(0)≥0,Ei(0)≥0,Ii(0)≥0,Ri(0)≥0,0≤Ai(0)≤Ci,F1i(0)≥0,F2i(0)≥0,F3i(0)≥0,F4i(0)≥0,U1i(0)≥0,U2i(0)≥0
|
Table 1. Parameters of the Model.
Parameters |
Description of parameters |
α |
Human recovery rate |
β1 |
Transmission rate of infection from female mosquitoes to human |
β2 |
Mosquitoes mating rate |
β3 |
Transmission rate of infection from human to female mosquito |
γ |
Transition rate from aquatic stage to adult mosquito |
μ |
Natural death rate of human |
ω |
Birth rate of human |
ω1 |
Constant recruitment rate of sterile male mosquito |
ϕ |
Recruitment rate for aquatic mosquito |
C |
Carrying capacity for aquatic/adult mosquito |
d |
Natural death rate of mosquito at aquatic state |
d1 |
Natural death rate of mosquito |
k |
Rate at which exposed human become infectious |
mij |
Migration rate from patch j to patch i |
p |
Proportion of female mosquito |
In the next section, the dynamics of the model (1)-(11) is discussed for an isolated patch (i.e mSij=mEij=mIij=mRij=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:
dAdt=ϕ(1−AC)(F2+F4)−(γ+d)A
|
(16)
|
dF1dt=pγA−β2F1U1U1+U2−β2F1U2U1+U2−d1F1
|
(17)
|
dF2dt=β2F1U1U1+U2−β3IF2−d1F2
|
(18)
|
dF3dt=β2F1U2U1+U2−d1F3
|
(19)
|
Let X(t)=(S(t),E(t),I(t),R(t),A(t),F1(t),F2(t),F3(t),F4(t),U1(t),U2(t))⊺∈R11+. The system (12)-(22) with non-zero initial conditions can be written in the following form:
dXdt=Q(X(t),t),X(0)=X0≥0Q(X(t),t)=(Q1(X,t),Q2(X,t),Q3(X,t)....,Q11(X,t))⊺
|
(23)
|
3.1. Positivity and boundedness of the solution
Proposition 1. The positive cone Int(R11+) is invariant for the system (23).
Proof. Observe that the boundaries of non-negative cone R11+ are invariant for the system (23) and Q(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 Ω= {(S,E,I,R,A,F1,F2,F3,F4,U1,U2)∈R11+|S+E+I+R≤ωμ, A(t)≤C,F1+F2+F3+F4+U1+U2≤ω1+γCd1}.
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
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
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_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.
Parameters |
Parameters values |
\alpha |
0.3 |
\beta^1 |
0.02 |
\beta^2 |
0.7 |
\beta^3 |
0.03 |
\gamma |
0.075 |
\omega |
0.002 |
\phi |
5 |
\mu |
0.0000456 |
C |
450 |
d |
0.05 |
d^{1} |
0.0714 |
k |
0.1667 |
p |
0.5 |
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.
Parameters |
Parameters values |
\alpha_i |
0.5 |
\beta_i^1 |
0.001 |
\beta_i^2 |
0.7 |
\beta_i^3 |
0.001 |
\gamma_i |
0.075 |
\omega_i |
0.029 |
\phi_i |
5 |
\mu_i |
0.0000456 |
C_i |
450 |
d_i |
0.05 |
d_i^{1} |
0.0714 |
k_i |
0.1667 |
p_i |
0.5 |
\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)).
Table 4. When SIT is applied only in the patch-1 i.e.(\omega^{1}_1=70, \omega^{1}_2=0).
Cases |
Migration in patch-1 |
Migration in patch-2 |
Conclusions |
(a) |
m_{12}=0 |
m_{21}=0 |
Isolated patches with R_0 of patch-1 is 0.2515<1 (37). Disease is controlled in patch-1 only. |
(b) |
m_{12}=0.005 |
m_{21}=0.0006 |
R_0^2 for two patch network model is 0.7678<1 (43). Two patch network model may be disease-free. The time-series confirms that the infection level tends to zero in both the patches Figure 8a. |
(c) |
m_{12}=0.009 |
m_{21}=0.0025 |
R_0^2 for two-patch network model is 1.4763>1 (43). SIT method will not be able to control disease with this migration combination Figure 8b. |
Table 5. When SIT is applied only in the patch-1 i.e.(\omega^{1}_1=60, \omega^{1}_2=0).
Cases |
Migration in patch-1 |
Migration in patch-2 |
Conclusions |
(d) |
m_{12}=0 |
m_{21}=0 |
Isolated patches with R_0 of patch-1 is 0.2515<1. Disease is controlled in patch-1 only. |
(e) |
m_{12}=0.005 |
m_{21}=0.001 |
R_0^2 for two-patch network model is 1.1778>1. SIT method will not be able to control disease with this migration combination. Disease may persist in the network Figure 8c. |
(f) |
m_{12}=0.009 |
m_{21}=0.001 |
R_0^2 for two-patch network model is 0.7273<1. The network may now be disease free Figure 8d. |
\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).
Table 6. When SIT is applied only in the patch-1 i.e.(\omega^{1}_1=70, \omega^{1}_2=0 and \omega^{1}_3=0).
Cases |
Migration in patch-1 |
Migration in patch-2 |
Migration in patch-3 |
R^n_0 for network model |
(g) |
m_{12}=0.0001, m_{13}=0.002 |
m_{21}=0.0001 |
m_{31}=0.0001 |
0.6762 (<1) |
(h) |
m_{12}=0.0002, m_{13}=0.0002 |
m_{21}=0.001 |
m_{31}=0.001 |
2.6813 (>1) |
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.