Long-lasting insecticidal nets and the quest for malaria eradication: a mathematical modeling approach

Abstract

Recent dramatic declines in global malaria burden and mortality can be largely attributed to the large-scale deployment of insecticidal-based measures, namely long-lasting insecticidal nets (LLINs) and indoor residual spraying. However, the sustainability of these gains, and the feasibility of global malaria eradication by 2040, may be affected by increasing insecticide resistance among the Anopheles malaria vector. We employ a new differential-equations based mathematical model, which incorporates the full, weather-dependent mosquito lifecycle, to assess the population-level impact of the large-scale use of LLINs, under different levels of Anopheles pyrethroid insecticide resistance, on malaria transmission dynamics and control in a community. Moreover, we describe the bednet-mosquito interaction using parameters that can be estimated from the large experimental hut trial literature under varying levels of effective pyrethroid resistance. An expression for the basic reproduction number, ${\mathcal {R}}_0$, as a function of population-level bednet coverage, is derived. It is shown, owing to the phenomenon of backward bifurcation, that ${\mathcal {R}}_0$ must be pushed appreciably below 1 to eliminate malaria in endemic areas, potentially complicating eradication efforts. Numerical simulations of the model suggest that, when the baseline ${\mathcal {R}}_0$ is high (corresponding roughly to holoendemic malaria), very high bednet coverage with highly effective nets is necessary to approach conditions for malaria elimination. Further, while >50% bednet coverage is likely sufficient to strongly control or eliminate malaria from areas with a mesoendemic malaria baseline, pyrethroid resistance could undermine control and elimination efforts even in this setting. Our simulations show that pyrethroid resistance in mosquitoes appreciably reduces bednet effectiveness across parameter space. This modeling study also suggests that increasing pre-bloodmeal deterrence of mosquitoes (deterring them from entry into protected homes) actually hampers elimination efforts, as it may focus mosquito biting onto a smaller unprotected host subpopulation. Finally, we observe that temperature affects malaria potential independently of bednet coverage and pyrethroid-resistance levels, with both climate change and pyrethroid resistance posing future threats to malaria control.

Appendices

Appendix A: Proof of Lemma 2.1

Proof

(a) It should be noted, first of all, that the right-hand side of each of the equations of the model (2.1), (2.2),(2.4) is continuous and locally-Lipschitz at $t=0$. Hence, a solution of the model with non-negative initial conditions exists and is unique in $\Omega =\Omega _1\times \Omega _2\times \Omega _3$ for all time $t>0$ (see also Iboi et al. 2018; Okuneye et al. 2019). Furthermore, since ${\left( 1-\frac{E}{K_E} \right) }_+\ge 0,$ it follows from the first equation of the sub-system (2.1) that $E(t)\le K_E$ for all time $t>0$. Similarly, it follows from the second equation of the sub-system (2.1) that

$$\begin{aligned} {\dot{L}}_{1}= \sigma _E E-\left( \sigma _{L_1}+\mu _L \right) L_1\le \sigma _E K_E-\left( \sigma _{L_1}+\mu _L \right) L_1, \end{aligned}$$

so that $\displaystyle \limsup _{t\rightarrow \infty }L_1(t)\le \frac{\sigma _E K_E}{\sigma _{L_1}+\mu _L}=L^{\diamond }_1.$ Using a similar approach, it can be shown that $\displaystyle \limsup _{t\rightarrow \infty }L_2(t)\le \frac{\sigma _{L_1} L^{\diamond }_1}{\sigma _{L_2}+\mu _L}=L^{\diamond }_2$, $\displaystyle \limsup _{t\rightarrow \infty }L_3(t)\le \frac{\sigma _{L_2}L^{\diamond }_2}{\sigma _{L_3}+\mu _L}=L^{\diamond }_3$, $\displaystyle \limsup _{t\rightarrow \infty }L_4(t)\le \frac{\sigma _{L_3}L^{\diamond }_3}{\sigma _{L_4}+\mu _L}=L^{\diamond }_4$ and $\displaystyle \limsup _{t\rightarrow \infty }P(t)\le \frac{\sigma _{L_4}L^{\diamond }_4}{\sigma _P +\mu _P}=P^{\diamond }.$ That is, all solutions of the sub-system (2.1) are bounded for all time $t>0$.

For the boundedness of the solutions of the sub-system (2.2), we consider the following equation (for the rate of change of the total adult mosquito population):

$$\begin{aligned} {\dot{N}}_M= & {} f\sigma _P P-{\mu _{\mathbf{X}}} S_X-{\mu _{\mathbf{X}}} E_X-{\mu _{\mathbf{X}}} I_X-\mu _M N_M\nonumber \\&+b_H(Q_2+Q_3-Q_1+R_1+R_2) A_X, \end{aligned}$$

(A-1)

where, $N_M=A_X+A_Y+A_Z$ (with $A_X, A_Y, A_Z, Q_1, Q_2, Q_3, R_1$ and $R_2$ are as defined in Sect. 2). It can be shown that $Q_2+Q_3-Q_1+R_1+R_2<0$. Hence, Eq. (A-1) can be re-written as

$$\begin{aligned} {\dot{N}}_M= & {} f\sigma _P P-{\mu _{\mathbf{X}}} S_X-{\mu _{\mathbf{X}}} E_X-{\mu _{\mathbf{X}}} I_X-\mu _M N_M\nonumber \\&+b_H(Q_2+Q_3-Q_1+R_1+R_2) A_X\le f\sigma _P P-\mu _MN_M, \end{aligned}$$

(A-2)

so that,

$$\begin{aligned} \displaystyle \limsup _{t\rightarrow \infty } N_M(t)\le \frac{f\sigma _P P^{\diamond }}{\mu _M}. \end{aligned}$$

Hence, the solutions of the equations of the sub-system (2.2) are bounded for all time $t>0$. Similarly, consider the equation for the rate of change of the total human population, given by:

$$\begin{aligned} {\dot{N}}_H=\Pi -\mu _H N_H-\delta _H (I_{H_p}+I_{H_u})\le \Pi -\mu _H N_H, \end{aligned}$$

(A-3)

from which it follows that $\displaystyle \limsup _{t\rightarrow \infty } N(t)\le \frac{\Pi }{\mu _H}.$ Thus, the solutions of the sub-system (2.4) are bounded for all $t>0$. Since the solutions of the three sub-systems of the model (2.1), (2.2),(2.4) are bounded, it follows that the solutions of the model are bounded. This concludes the proof of Item (a).

(b) The proof for the invariance of the region $\Omega _1$ follows from the bounds established in Item (a) (i.e., $0<\displaystyle \limsup _{t\rightarrow \infty }L_1(t)\le L^{\diamond }_1$, $0<\displaystyle \limsup _{t\rightarrow \infty }L_j(t)\le L^{\diamond }_j$ and $0<\displaystyle \limsup _{t\rightarrow \infty }P(t)\le P^{\diamond }$) and the fact that ${\dot{E}}(t)<0$ whenever $E(t)>K_E$, ${\dot{L}}_1(t)<0$ whenever $L_1(t)>L^{\diamond }_1$ and ${\dot{L}}_j(t)<0$ whenever $L_j(t)>L^{\diamond }_1$ ($j=2,3,4$), respectively.

For the invariance of the region $\Omega _2$, it is convenient to consider the following equation for the rate of change of the total mosquito population given by:

$$\begin{aligned} {\dot{N}}_M= & {} f\sigma _P P-{\mu _{\mathbf{X}}} S_X-{\mu _{\mathbf{X}}} E_X-{\mu _{\mathbf{X}}} I_X-\mu _M N_M\\&+b_H(Q_2+Q_3-Q_1+R_1+R_2) A_X\le f\sigma _P P-\mu _M N_M. \end{aligned}$$

It follows that ${\dot{N}}_M<0$ whenever $N_M(t)> \frac{f\sigma _P P^{\diamond }}{\mu _M}$. Thus, the region $\Omega _2$ is invariant with respect to the sub-system (2.2) of the model (2.1), (2.2),(2.4).

Finally, consider the equation for the total human population given by

$$\begin{aligned} {\dot{N}}_H=\Pi -\mu _H N_H-\delta _H(I_{Hp}+I_{Hu})\le \Pi -\mu _H N_H. \end{aligned}$$

It follows that ${\dot{N}}_H<0$ whenever $N_H(t)> \frac{\Pi }{\mu _H}$. Thus, the region $\Omega _3$ is invariant with respect to the sub-system (2.4) of the model (2.1), (2.2),(2.4). Since the regions $\Omega _1$, $\Omega _2$ and $\Omega _3$ are positively-invariant and attracting, it follows that $\Omega =\Omega _1\times \Omega _2\times \Omega _3$ is positively-invariant and attracting for the model (2.1), (2.2),(2.4). This concludes the proof of Item (b). $\square$

Appendix B: Coefficients of Eq. (3.5)

$$\begin{aligned} C_1= & {} b_H\left\{ \pi _p(1 - \varepsilon _{deter})[(1 - \varepsilon _{die,p})\varepsilon _{bite|\sim die,p}+\varepsilon _{die,p}]\right. \\&\left. +\pi _u[(1 - \varepsilon _{die,u}) \varepsilon _{bite|\sim die,u}+\varepsilon _{die,u}]\right\} +{\mu _{\mathbf{X}}}+\mu _M>0, \\ C_2= & {} b_H[\pi _p(1 - \varepsilon _{deter})(1 - \varepsilon _{die,p}) \varepsilon _{bite|\sim die,p}+\pi _u(1 - \varepsilon _{die,u}) \varepsilon _{bite|\sim die,u}]>0, \\ C_3= & {} b_H\left\{ \pi _p(1 - \varepsilon _{deter})[(1 - \varepsilon _{die,p}) \varepsilon _{bite|\sim die,p}+\varepsilon _{die,p}]\right. \\&\left. +\pi _u[(1 - \varepsilon _{die,u}) \varepsilon _{bite|\sim die,u}+\varepsilon _{die,u}]\right\} +\kappa _V+{\mu _{\mathbf{X}}}+\mu _M>0, \\ C_4= & {} \left\{ \pi _p (1 - \varepsilon _{deter})[\varepsilon _{die,p} (1 - \varepsilon _{bite|\sim die,p})+\varepsilon _{bite|\sim die,p}]\right. \\&\left. +\pi _u [\varepsilon _{die,u} (1 - \varepsilon _{bite|\sim die,u})+\varepsilon _{bite|\sim die,u}]\right\} \kappa ^2_V>0, \\ C_5= & {} 2\kappa _V\left( \mu _M+\frac{\theta _Y}{2}+\frac{\varphi _Z}{2}\right) C_4>0, \\ C_6= & {} \pi _p(1 - \varepsilon _{deter})\left\{ [\varepsilon _{die,p} (1 - \varepsilon _{bite|\sim die,p})+\varepsilon _{bite|\sim die,p}]\mu ^2_M\right. \\&\left. +(\theta _Y+\varphi _Z)[\varepsilon _{die,u} (1 - \varepsilon _{bite|\sim die,u})+\varepsilon _{bite|\sim die,u}]\mu _M+\varepsilon _{bite|\sim die,p}\theta _Y\varphi _Z\right\}>0, \\ C_7= & {} \pi _u\left\{ [\varepsilon _{die,u} (1 - \varepsilon _{bite|\sim die,u})+\varepsilon _{bite|\sim die,u}]\mu ^2_M\right. \\&\left. +(\theta _Y+\varphi _Z)[\varepsilon _{die,u} (1 - \varepsilon _{bite|\sim die,u})+\varepsilon _{bite|\sim die,u}]\mu _M+\varepsilon _{bite|\sim die,u}\theta _Y\varphi _Z\right\}>0, \\ C_8= & {} (\mu _M+\kappa _V)K_{10}K_{12}, \\ C_{9}= & {} \varphi _Z+\theta _Y+b_H\left\{ \pi _p(1 - \varepsilon _{deter})[1-(1 - \varepsilon _{die,p}) (1 - \varepsilon _{bite|\sim die,p})]\right. \\&\left. +\pi _u[1-(1 - \varepsilon _{die,u}) (1 - \varepsilon _{bite|\sim die,u})]\right\}>0,\\ C_{10}= & {} \theta _Y\varphi _Z+b_H(\theta _Y+\varphi _Z)\left\{ \pi _p(1 - \varepsilon _{deter})[1-(1 - \varepsilon _{die,p}) (1 - \varepsilon _{bite|\sim die,p})]\right. \\&\left. +\pi _u[1-(1 - \varepsilon _{die,u}) (1 - \varepsilon _{bite|\sim die,u})]\right\}>0 , \\ C_{11}= & {} b_H\theta _Y\varphi _Z[\pi _p\varepsilon _{die,p}(1 - \varepsilon _{deter})+\varepsilon _{die,u}\pi _u]>0. \end{aligned}$$

Appendix C: Equations of the model (2.1), (2.2), (2.4) without bednets intervention

In the absence of the bednet-based intervention, the sub-systems of the model involving the adults and human dynamics, given by Eqs. (2.2) and (2.4), reduce, respectively, to the following sub-systems:

$$\begin{aligned} \begin{aligned} {\mathrm{Stage}} \,\,\,\,\text {I} {\left\{ \begin{array}{ll} {\dot{S}}_{X} &{}=f\sigma _P P+\varphi _Z S_Z+ b_H Q_3 S_X-\left( b_H+{\mu _{\mathbf{X}}}+\mu _M\right) S_{X}, \\ {\dot{E}}_{X} &{}= \varphi _Z E_Z+ b_H Q_3 E_X -\left( b_H+\kappa _V+{\mu _{\mathbf{X}}}+\mu _M\right) E_{X}, \\ {\dot{I}}_{X} &{}=\varphi _Z I_Z+\kappa _V E_{X}+ b_H Q_3 I_X-\left( b_H+{\mu _{\mathbf{X}}}+\mu _M\right) I_{X}.\\ \end{array}\right. }\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \end{aligned} \end{aligned}$$

$$\begin{aligned} \begin{aligned} {\mathrm{Stage}}\,\,\,\, \text {II} {\left\{ \begin{array}{ll} {\dot{S}}_{Y} &{}= (1 - \beta _V \omega ) R_2 S_{X} -\left( \theta _Y+\mu _M\right) S_{Y}, \\ {\dot{E}}_{Y} &{}= b_H \beta _V \omega R_2 S_{X} + b_H R_2 E_X-\left( \theta _Y+\kappa _V+\mu _M\right) E_{Y}, \\ {\dot{I}}_{Y} &{}=\kappa _V E_Y+b_H R_2 I_X-\left( \theta _Y+\mu _M\right) I_{Y}.\\ \end{array}\right. }\quad \end{aligned} \end{aligned}$$

(C-1)

$$\begin{aligned} \begin{aligned} {\mathrm{Stage}}\,\,\,\, \text {III} {\left\{ \begin{array}{ll} {\dot{S}}_{Z} &{}=\theta _YS_{Y} -\left( \varphi _Z+\mu _M\right) S_{Z}, \\ {\dot{E}}_{Z} &{}= \theta _Y E_{Y}-\left( \varphi _Z+\kappa _V+\mu _M\right) E_{Z}, \\ {\dot{I}}_{Z} &{}=\theta _Y I_{Y}+\kappa _V E_Z-\left( \varphi _Z+\mu _M\right) I_{Z},\\ \end{array}\right. } \end{aligned} \end{aligned}$$

$$\begin{aligned} \begin{aligned} {\mathrm{Human}}\,\,\,\, {\left\{ \begin{array}{ll} {\dot{S}}_{H} &{}=\Pi - (\lambda _{VH}+\mu _{H})S_{H}+\eta _H R_{H}, \\ {\dot{E}}_{H} &{}=\lambda _{VH}S_{H}-(\gamma _H+\mu _H)E_{H}, \\ {\dot{I}}_{H} &{}=\gamma _HE_{H}-(\alpha _H+\delta _H+\mu _H)I_{H}, \\ {\dot{R}}_{H} &{}=\alpha _HI_{H}-(\eta _H+\mu _{H})R_{H}. \\ \end{array}\right. } \end{aligned} \end{aligned}$$

(C-2)

The equations for the aquatic dynamics, given by (2.1), remain unchanged. Hence, the reduced (no-bednets) model consist of the Eqs. (2.1), (C-1), (C-2).

It can be shown, using the next generation operator method (as in Sect. 3), that the basic reproduction number of the reduced model (2.1), (C-1), (C-2) is given by

$$\begin{aligned} \tilde{{\mathcal {R}}}_{0*}=\sqrt{\frac{b_H\beta _VJ_1J_3S^0_X\beta _M\gamma _H\theta _Y\varphi _Z\kappa _V[(K_9+K_{12})J_5+K_9K_{11}]}{K_{13}K_{14}N^*_H(J_5K_{10}K_{12}-J_6\theta _Y\varphi _Z)(J_4K_{9}K_{11}-J_6\theta _Y\varphi _Z)}}, \end{aligned}$$

(C-3)

where,

$$\begin{aligned} \begin{aligned} J_1&= b_H\left[ \varepsilon _{bite|die,u} \,\varepsilon _{die,u} + \varepsilon _{bite|\sim die,u} (1 - \varepsilon _{die,u})\right] ,\,\,\\ J_2&= (1 - \varepsilon _{die,u}) (1 - \varepsilon _{bite|\sim die,u}),\,\,N^*_H=\frac{\Pi }{\mu _H},\\ J_3&= (1 - \varepsilon _{die,u}) \varepsilon _{bite|\sim die,u},\,\,J_4=(b_H+{\mu _{\mathbf{X}}}+\mu _M)-b_HJ_2,\,\,\\ J_5&=(b_H+\kappa _V+{\mu _{\mathbf{X}}}+\mu _M)-b_HJ_2,\\ J_6&=b_HJ_3,\,\, {\mathcal {N}}_{0*}=\frac{(\psi _E\varphi _Z\sigma _Ef\sigma _P\theta _YJ_6)\prod \limits ^4_{i=1}\sigma _{L_{i}}}{\left( J_4K_9K_{11}-J_6\theta _Y\varphi _Z\right) \prod \limits ^6_{i=1} K_i},\,\,\\ S^0_X&=\frac{\left[ f\sigma _E\sigma _PK_E\left( 1-\frac{1}{{\mathcal {N}}_{0*}}\right) K_9K_{11}\right] \prod \limits ^4_{i=1}\sigma _{L_{i}}}{(J_4K_9K_{11}-J_6\theta _Y\varphi _Z)\prod \limits ^6_{i=2} K_i}. \end{aligned} \end{aligned}$$

Substituting the baseline parameter values in Table 4 (for the holo-endemic setting) shows that the worst-case scenario basic reproduction number $(\tilde{{\mathcal {R}}}_0)$ of the model (2.1), (2.2), (2.4), or, equivalently, the reduced model (C-1), (C-2), given by (C-3), is $\tilde{{\mathcal {R}}}_{0*}$ $=$ 11.4.

Appendix D: Proof of Theorem 3.2

Proof

The proof of Theorem 3.2 is based on using center manifold theory (Carr 1981; Castillo-Chavez and Song 2004). It is convenient to define the following change of variables for the model (2.1), (2.2), (2.4): $E=x_1$, $L_{1}=x_2$, $L_{2}=x_3$, $L_{3}=x_4$, $L_{4}=x_5$, $P=x_6$, $S_{X}=x_7$, $E_{X}=x_8$, $I_{X}=x_9$, $S_{Y}=x_{10}$, $E_{Y}=x_{11}$, $I_{Y}=x_{12}$, $S_{Z}=x_{13}$, $E_{Z}=x_{14}$, $I_{Z}=x_{15}$, $S_{Hp}=x_{16}$, $E_{Hp}=x_{17}$, $I_{Hp}=x_{18}$, $R_{Hp}=x_{19}$, $S_{Hu}=x_{20}$, $E_{Hu}=x_{21}$, $I_{Hu}=x_{22}$, $R_{Hu}=x_{23}$. Using the vector notation $X=(x_1,\ldots x_{23})^T$ and $F=(f_1,\ldots f_{23})^T$, the model can then be written in the form $\frac{dX}{dt} = (f_1,\ldots f_{23})^T$, as follows:

$$\begin{aligned} \begin{aligned} {\dot{x}}_1 \equiv f_1&=\psi _E\varphi _Z{\left( 1-\frac{x_1}{K_E} \right) }_+(x_{13}+x_{14}+x_{15}) -\left( \sigma _E+\mu _E\right) x_1,\\ {\dot{x}}_2 \equiv f_2&= \sigma _E x_1-\left( \sigma _{L_1}+\mu _L \right) x_2,\\ {\dot{x}}_3 \equiv f_3&= \sigma _{L_1}x_2-\left( \sigma _{L_2}+\mu _L \right) x_3,\\ {\dot{x}}_4 \equiv f_4&= \sigma _{L_2}x_3-\left( \sigma _{L_3}+\mu _L \right) x_4, \\ {\dot{x}}_5 \equiv f_5&= \sigma _{L_3}x_4-\left( \sigma _{L_4}+\mu _L \right) x_5,\\ {\dot{x}}_6 \equiv f_6&= \sigma _{ L_{4}}x_5-\left( \sigma _P +\mu _P \right) x_6,\\ {\dot{x}}_7 \equiv f_7&= f\sigma _P x_6+\varphi _Z x_{13}+ b_H (Q_2+Q_3) x_7-\left( b_H Q_1+{\mu _X}+\mu _M\right) x_{7},\\ {\dot{x}}_8 \equiv f_8&= \varphi _Z x_{14}+ b_H (Q_2+Q_3) x_8 -\left( b_H Q_1+\kappa _V+{\mu _X}+\mu _M\right) x_{8}, \\ {\dot{x}}_9 \equiv f_9&= \varphi _Z x_{15}+\kappa _V x_{8}+ b_H (Q_2+Q_3) x_9-\left( b_H Q_1+{\mu _X}+\mu _M\right) x_{9},\\ {\dot{x}}_{10} \equiv f_{10}&= b_H [(1 - \beta _V \omega _p) R_1 + (1 - \beta _V \omega _u) R_2] x_{7} -\left( \theta _Y+\mu _M\right) x_{10}, \\ {\dot{x}}_{11 } \equiv f_{11}&= b_H (\beta _V \omega _p R_1 + \beta _V \omega _u R_2) x_{7}\\&\quad + b_H (R_1+R_2) x_8-\left( \theta _Y+\kappa _V+\mu _M\right) x_{11}, \\ {\dot{x}}_{12} \equiv f_{12}&= \kappa _V x_{11}+b_H(R_1+R_2) x_9-\left( \theta _Y+\mu _M\right) x_{12},\\ {\dot{x}}_{13} \equiv f_{13}&= \theta _Y x_{10} -\left( \varphi _Z+\mu _M\right) x_{13}, \\ {\dot{x}}_{14} \equiv f_{14}&= \theta _Y x_{11}-\left( \varphi _Z+\kappa _V+\mu _M\right) x_{14}, \\ {\dot{x}}_{15} \equiv f_{15}&= \theta _Y x_{12}+\kappa _V x_{14}-\left( \varphi _Z+\mu _M\right) x_{15},\\ {\dot{x}}_{16} \equiv f_{20}&= \Pi \pi _p - (\lambda _{VH_p}+\mu _{H})x_{16}+\eta _H x_{19}, \\ {\dot{x}}_{17} \equiv f_{21}&= \lambda _{VH_p}x_{16}-(\gamma _H+\mu _H)x_{17}, \\ {\dot{x}}_{18} \equiv f_{22}&= \gamma _H x_{17}-(\alpha _H+\delta _H+\mu _H)x_{18}, \\ {\dot{x}}_{19} \equiv f_{23}&=\alpha _H x_{18}-(\eta _H+\mu _{H})x_{19}, \\ {\dot{x}}_{20} \equiv f_{24}&= \Pi \pi _u - (\lambda _{VH_u}+\mu _{H})x_{20}+\eta _H x_{23}, \\ {\dot{x}}_{21} \equiv f_{25}&= \lambda _{VH_u}x_{20}-(\gamma _H+\mu _H)x_{21}, \\ {\dot{x}}_{22} \equiv f_{26}&= \gamma _Hx_{21}-(\alpha _H+\delta _H+\mu _H)x_{22}, \\ {\dot{x}}_{23} \equiv f_{27}&= \alpha _Hx_{22}-(\eta _H+\mu _{H})x_{23}, \end{aligned} \end{aligned}$$

(D-1)

where,

$$\begin{aligned} \text {EIR}_p= & {} b_H \frac{x_9}{N_{H_p}} \pi _p(1 - \varepsilon _{deter}) \left[ \varepsilon _{bite|die,p} \varepsilon _{die,p} + \epsilon _{bite|\sim die,p} (1 - \varepsilon _{die,p})\right] , \\ \text {EIR}_u= & {} b_H \frac{x_9}{N_{H_u}}\pi _u \left[ \varepsilon _{bite|die,u} \varepsilon _{die,u} + \varepsilon _{bite|\sim die,u} (1 - \varepsilon _{die,u})\right] ,\\ \lambda _{VH_p}= & {} {\beta _M\,\text {EIR}_p},\\ \lambda _{VH_u}= & {} \beta _M\,\text {EIR}_u,\\ \omega _p= & {} \frac{x_{18}}{N_{H_p}},\quad \omega _u= \frac{x_{22}}{N_{H_u}}. \\ \end{aligned}$$

Let ${\mathcal {R}}_0=1$ and suppose, further, that $\beta _M=\beta _M^*$ is chosen as a bifurcation parameter. Solving for $\beta _M=\beta _M^*$ from ${\mathcal {R}}_0=1$ gives

$$\begin{aligned} \beta _M=\beta _M^*=\frac{\Pi \pi _p\pi _u\left( C_3K_{10}K_{12}-C_2\theta _Y\varphi _Z\right) \left( C_1K_9K_{11}-C_2\theta _Y\varphi _Z\right) }{\mu ^2_Hb_H x^*_7\kappa _V \varphi _Z \theta _Y[(K_9+K_{12})C_3+K_9K_{11}]({\mathcal {R}}_{H_pV}+{\mathcal {R}}_{H_uV})}. \end{aligned}$$

The Jacobian of the transformed system (D-1), evaluated at the DFE $({\mathcal {T}}_2)$ with $\beta _M=\beta _M^*$, is given by

$$\begin{aligned} {J(\beta _M^*)}= \left[ \begin{array}{*{20}c} J_1&{}\quad J_2\\ J_3&{}\quad J_4\ \end{array} \right] , \end{aligned}$$

where,

$$\begin{aligned} {J_1}= & {} \left[ \begin{array}{*{20}c} -\frac{\psi _E\varphi _Zx^*_{13}}{K_E}-K_1&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ \sigma _E&{}\quad -K_2&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad \sigma _{L_1}&{}\quad -K_3&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad \sigma _{L_2}&{}\quad -K_4&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad \sigma _{L_3}&{}\quad -\mu _H&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \sigma _{L_4}&{}\quad -K_5&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad f\sigma _{P}&{}\quad -C_1&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad -C_3&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \kappa _V&{}-C_1&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad C_2&{}\quad 0&{}\quad 0&{}\quad -K_9&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad C_2&{}\quad 0&{}\quad 0&{}\quad -K_{10}&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad C_2&{}\quad 0&{}\quad \kappa _V&{}\quad -K_9\ \end{array} \right] ,\\ {J_2}= & {} \left[ \begin{array}{*{20}c} \psi _E\varphi _Z(1-\frac{x^*_1}{K_E})&{}\quad \psi _E\varphi _Z(1-\frac{x^*_1}{K_E})&{}\quad \psi _E\varphi _Z(1-\frac{x^*_1}{K_E})&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ \varphi _Z&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad \varphi _Z&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad \varphi _Z&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad -\frac{b_HR_1\beta _V x^*_7}{x^*_{16}}&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad -\frac{b_HR_2\beta _V x^*_7}{x^*_{20}}&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \frac{b_HR_1\beta _V x^*_7}{x^*_{16}}&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \frac{b_HR_2\beta _V x^*_7}{x^*_{20}}&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\ \end{array} \right] ,\\ {J_3}= & {} \left[ \begin{array}{*{20}c} 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \theta _Y&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \theta _Y&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \theta _Y\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad -\beta _MQ_p&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \beta _MQ_p&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad -\beta _MQ_u&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \beta _MQ_u&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\ \end{array} \right] , \end{aligned}$$

and,

$$\begin{aligned} {J_4}= \left[ \begin{array}{*{20}c} -K_{11}&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad -K_{12}&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad \kappa _V&{}\quad -K_{11}&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad -\mu _H&{}\quad 0&{}\quad 0&{}\quad \eta _H&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad -K_{13}&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \gamma _H&{}\quad -K_{14}&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \alpha _H&{}\quad -K_{15}&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad -\mu _H&{}\quad 0&{}\quad 0&{}\quad \eta _H\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad -K_{13}&{}\quad 0&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \gamma _H&{}\quad -K_{14}&{}\quad 0\\ 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad 0&{}\quad \alpha _H&{}\quad -K_{15}\ \end{array} \right] . \end{aligned}$$

The Jacobian $J(\beta _M^*)$ has a simple zero eigenvalue (and all other eigenvalues having negative real parts). Hence, the center manifold theory Carr (1981); Castillo-Chavez and Song (2004) can be used to analyse the dynamics of (D-1) near $\beta _M=\beta _M^*$. This entails carrying out the following computations.

Eigenvectors of$J({\mathcal {T}}_2)\mid _{\beta _M=\beta _M^*}$ The Jacobian of the transformed system (D-1), evaluated at the DFE $({\mathcal {T}}_2)$ with $\beta _M=\beta _M^*$, has a right and left eigenvectors (associated with the zero eigenvalue) (the expression for the eigenvector $w_i$ and $v_i$, $i=1,2,..,23$, are given in the Supplementary Material).

Computations of bifurcation coefficients of a and b

By computing the associated non-zero partial derivatives of F(x) evaluated the the DFE, it follows from Theorem 4.1 in Castillo-Chavez and Song (2004) that the associated bifurcation coefficients, a and b, are given, respectively, by

$$\begin{aligned} \begin{aligned} a=&-2 b_H \left[ \frac{\left( -w_{{22}}R_{{2}}w_{{7}}v_{{11}}\beta _{{V}}+v_{{21}}Q_{{u}} \beta _{{M}}w_{{9}} \left( w_{{22}}+w_{{21}}+w_{{23}} \right) \right) x^*_{{20}}+w_{{22}}\beta _{{V}}R_{{2}}x^*_{{7}}v_{{11}} \left( w_{{20}}+w_{ {21}}+w_{{22}}+w_{{23}} \right) }{(x^*_{20})^2}\right] \\&-2 b_H \left[ \frac{\left( -w_{{18}}\beta _{{V}}R_{{1}}w_{{7}}v_{{11}}+Q_{{p}}\beta _{{M}}w _{{9}} \left( w_{{18}}+w_{{19}}+1 \right) \right) x^*_{{16}}+w_{{18}}{(x^* _{{20}}})^2\beta _{{V}}R_{{1}}x^*_{{7}}v_{{11}} \left( w_{{16}}+w_{{18}} +w_{{19}}+1 \right) }{(x^*_{16})^2}\right] , \end{aligned}\nonumber \\ \end{aligned}$$

(D-2)

and,

$$\begin{aligned} b=b_{{H}}w_{{9}} \left( Q_{{u}}v_{{21}}+Q_{{p}} \right) >0. \end{aligned}$$

(D-3)

Hence, it follows from Theorem 4.1 of Castillo-Chavez and Song (2004) that the transformed model (D-1) undegoes a backward bifurcation at ${\mathcal {R}}_0=1$ if the bifurcation coefficient a (given by (D-2)) is positive. This concludes the proof. $\square$

Appendix E: Proof of Theorem 3.3.

Proof

Consider the special case of the model (2.1), (2.2), (2.4) without disease-induced mortality in the host population (i.e., $\delta _H=0$). Further, let ${\mathcal {N}}_0>1$ (so that the NDFE, ${\mathcal {T}}_2$, exists) and $\tilde{{\mathcal {R}}}_0\le 1$. Setting $\delta _H=0$ in the model (2.1), (2.2), (2.4) gives $N_H(t)\rightarrow \frac{\Pi }{\mu _H}$ as $t \rightarrow \infty$, and ${\bar{K}}_{14}=\alpha _H+\mu _H$. Hence, from now on, $N_H(t)$ is replaced by its limiting value, $\frac{\Pi }{\mu _H}$. Furthermore, consider the following linear Lyapunov function:

$$\begin{aligned} \mathcal {{F}}= & {} g_1 E_{Hp}+g_2 {I}_{Hp}+g_3 {E}_{Hu}+g_4 {I}_{Hu}+g_5 {E}_{X}+g_6 {I}_{X}\\&+g_7 {E}_{Y}+g_8 {I}_{Y}+g_9 {E}_{Z}+g_{10} {I}_{Z}, \end{aligned}$$

where,

$$\begin{aligned} \begin{aligned} g_1&=\beta _{{V}}N_{{{ Hu}}}\gamma _{{H}}b_{{H}}R_{{1}}S_{{X}}\kappa _{{V}} \varphi _{{Z}}\theta _{{Y}} \left( C_{{3}} \left( K_{{9}}+K_{{12}} \right) +K_{{9}}K_{{11}} \right) ,\\ g_2&=\beta _{{V}}N_{{{ Hu}}}b_{{H}}R_{{1}}S_{{X}}\kappa _{{V}}\varphi _{{Z} }\theta _{{Y}} \left( C_{{3}} \left( K_{{9}}+K_{{12}} \right) +K_{{9}}K _{{11}} \right) K_{{13}} ,\\ g_3&=\beta _{{V}}N_{{{ Hp}}}\gamma _{{H}}b_{{H}}R_{{2}}S_{{X}}\kappa _{{V}} \varphi _{{Z}}\theta _{{Y}} \left( C_{{3}} \left( K_{{9}}+K_{{12}} \right) +K_{{9}}K_{{11}} \right) ,\\ g_4&=\beta _{{V}}N_{{{ Hp}}}b_{{H}}R_{{2}}S_{{X}}\kappa _{{V}}\varphi _{{Z} }\theta _{{Y}} \left( C_{{3}} \left( K_{{9}}+K_{{12}} \right) +K_{{9}}K _{{11}} \right) K_{{13}} ,\\ g_5&={\frac{C_{{2}}\beta _{{M}}\gamma _{{H}}{\theta ^2_{{Y}}}{\varphi ^2_{{Z}} }{\kappa ^2_{{V}}}S_{{X}}\beta _{{V}} \left( C_{{3}} \left( K_{{9} }+K_{{12}} \right) +K_{{9}}K_{{11}} \right) ^{2} \left( N_{{{ Hp}}} Q_{{u}}R_{{2}}S_{{{ Hu}}}+N_{{{ Hu}}}Q_{{p}}R_{{1}}S_{{{ Hp}} } \right) }{ \left( C_{{3}}K_{{10}}K _{{12}}-C_{{2}}\theta _{{Y}}\varphi _{{Z}} \right) \left( C_{{1}}K_{{9}}K_{{11}}-C_{{2}}\theta _{{Y}} \varphi _{{Z}} \right) C_{{3}}}}\\&\quad +{\frac{K_{{9}}K_{{11}}K_{{13}}{\bar{K}}_{{14}}N_{{{ Hp}}}N_{{{ Hu}}} \left( C_{{3}}K_{{10}}K_{{12}}-C_{{2}}\theta _{{Y}}\varphi _{{Z}} \right) \kappa _{{V}}}{C_{{3}}}} ,\\ g_6&=K_{{9}}K_{{11}}K_{{13}}{\bar{K}}_{{14}}N_{{{ Hp}}}N_{{{ Hu}}} \left( C_{{3}}K_{{10}}K_{{12}}-C_{{2}}\theta _{{Y}}\varphi _{{Z}} \right) ,\\ g_7&={\frac{\beta _{{M}}\gamma _{{H}}{\theta ^2_{{Y}}}{\varphi ^2_{{Z}}}{ \kappa ^2_{{V}}}S_{{X}}\beta _{{V}} \left( C_{{3}} \left( K_{{9}}+K_{{ 12}} \right) +K_{{9}}K_{{11}} \right) ^{2} \left( N_{{{ Hp}}}Q_{{u} }R_{{2}}S_{{{ Hu}}}+N_{{{ Hu}}}Q_{{p}}R_{{1}}S_{{{ Hp}}} \right) }{ \left( C_{{3}}K_{{10}}K_ {{12}}-C_{{2}}\theta _{{Y}}\varphi _{{Z}} \right) \left( C_{{1}}K_{{9}}K_{{11}}-C_{{2}}\theta _{{Y}} \varphi _{{Z}} \right) }} ,\\ g_8&=\theta _{{Y}}K_{{13}}{\bar{K}}_{{14}}\varphi _{{Z}}N_{{{ Hp}}}N_{{{ Hu}}} \left( C_{{3}}K_{{10}}K_{{12}} -C_{{2}}\theta _{{Y}}\varphi _{{Z}}\right) ,\\ g_9&={\frac{C_{{2}}\beta _{{M}}\gamma _{{H}}{\theta ^{2}_{{Y}}}{\varphi ^{3}_{{Z}} }{\kappa ^{2}_{{V}}}S_{{X}}\beta _{{V}} \left( C_{{3}} \left( K_{{9} }+K_{{12}} \right) +K_{{9}}K_{{11}} \right) ^{2} \left( N_{{{ Hp}}} Q_{{u}}R_{{2}}S_{{{ Hu}}}+N_{{{ Hu}}}Q_{{p}}R_{{1}}S_{{{ Hp}} } \right) }{ \left( C_{{3}}K_{{10}}K _{{12}}-C_{{2}}\theta _{{Y}}\varphi _{{Z}} \right) \left( C_{{1}}K_{{9}}K_{{11}}-C_{{2}}\theta _{{Y}} \varphi _{{Z}} \right) K_{{12}}C_{{3}}}}\\&\quad +{\frac{\varphi _{{Z}}K_{{9}}K_{{13}}{\bar{K}}_{{14}}N_{{{ Hp}}}N_{{{ Hu }}} \left( C_{{3}}K_{{10}}K_{{12}}-C_{{2}}\theta _{{Y}}\varphi _{{Z}} \right) \kappa _{{V}} \left( K_{{11}}+C_{{3}} \right) }{K_{{12}}C_{{3} }}} ,\,\,\\ \mathrm{and,}\,\, g_{10}&=\varphi _{{Z}}K_{{9}}K_{{13}}{\bar{K}}_{{14}}N_{{{ Hp}}}N_{{{ Hu}}} \left( C_{{3}}K_{{10}}K_{{12}}-C_{{2}}\theta _{{Y}}\varphi _{{Z}} \right) . \end{aligned} \end{aligned}$$

The Lyapunov derivative of ${\mathcal {F}}$ (where a dot represents differentiation with respect to t) is given by:

$$\begin{aligned} \mathcal {{\dot{F}}}&=g_1 {\dot{E}}_{Hp}+g_2 {\dot{I}}_{Hp}+g_3 {\dot{E}}_{Hu}+g_4 {\dot{I}}_{Hu}+g_5 {\dot{E}}_{X}+g_6 {\dot{I}}_{X}\\&\quad +g_7 {\dot{E}}_{Y}+g_8 {\dot{I}}_{Y}+g_9 {\dot{E}}_{Z}+g_{10} {\dot{I}}_{Z},\\&=g_1 (\lambda _{VH_p}S_{H_p}-K_{13}E_{H_p})+g_2 (\gamma _HE_{H_p}-{\bar{K}}_{14}I_{H_p})+g_3 (\lambda _{VH_u}S_{H_u}-K_{13}E_{H_u})\\&\quad +g_4 (\gamma _HE_{H_u}-{\bar{K}}_{14}I_{Hu})+g_5 (\varphi _Z E_Z-C_3E_{X})+g_6 (\varphi _Z I_Z+\kappa _V E_{X}-C_1I_{X})\\&\quad +g_7 [(b_H\beta _V \omega _p R_1 +b_H \beta _V \omega _u R_2) S_{X} + C_2 E_X-K_{10}E_{Y}]\\&\quad +g_8 (\kappa _V E_Y+C_2 I_X-K_9I_{Y})\\&\quad +g_9 (\theta _Y E_{Y}-K_{12}E_{Z})+g_{10} (\theta _Y I_{Y}+\kappa _V E_Z-K_{11}I_{Z}),\\&=\left( \frac{g_1\beta _MQ_pS_{Hp}}{N_{Hp}}+\frac{g_3\beta _MQ_uS_{Hu}}{N_{Hu}}-g_6C_1+g_8C_2\right) I_X\\&\quad +\left( \frac{g_7b_H\beta _VR_1S_X}{N_{Hp}}-g_2{\bar{K}}_{14}\right) I_{Hp}\\&\quad +\left( \frac{g_7b_H\beta _VR_2S_X}{N_{Hu}}-g_4{\bar{K}}_{14}\right) I_{Hu}+\left( -g_1K_{13}+g_2\gamma _H\right) E_{Hp}\\&\quad +\left( -g_3K_{13}+g_4\gamma _H\right) E_{Hu}\\&\quad +\left( -g_5C_3+g_6\kappa _V+g_7C_2\right) E_{X}+\left( -g_7K_{10}+g_8\kappa _V+g_9\theta _Y\right) E_{Y}\\&\quad +\left( -g_9K_{12}+g_{10}\kappa _V+g_5\varphi _Z\right) E_{Z}\\&\quad +(-g_8K_9+g_{10}\theta _Y)I_Y+(-g_{10}K_{11}+g_6\varphi _Z)I_Z.\\ \end{aligned}$$

Since $S_{Hp}(t)\le N_{Hp}$, $S_{Hu}(t)\le N_{Hu}$, $N_{Hp}(t)=\frac{\Pi \,\pi _p}{\mu _H}$ and $N_{Hu}(t)=\frac{\Pi \,\pi _u}{\mu _H}$ in $\Omega$ for all $t>0$, it follows that:

$$\begin{aligned} \begin{aligned} \mathcal {{\dot{F}}}&=\le K_{13}{\bar{K}}_{14}N_{Hp}N_{Hu}(C_1K_9K_{11}-C_2\theta _Y\varphi _Z)(C_3K_{10}K_{12}-C_2\theta _Y\varphi _Z)(\tilde{{\mathcal {R}}}^2_0-1)I_X\\&\quad +\beta _{{V}}N_{{{ Hu}}}b_{{H}}R_{{1}}S_{{X}}\kappa _{{V}}\varphi _{{Z} }\theta _{{Y}} \left( C_{{3}} \left( K_{{9}}+K_{{12}} \right) +K_{{9}}K _{{11}} \right) K_{{13}}K_{{14}} (\tilde{{\mathcal {R}}}^2_0-1)I_{Hp}\\&\quad +\beta _{{V}}N_{{{ Hp}}}b_{{H}}R_{{2}}S_{{X}}\kappa _{{V}}\varphi _{{Z} }\theta _{{Y}} \left( C_{{3}} \left( K_{{9}}+K_{{12}} \right) +K_{{9}}K _{{11}} \right) K_{{13}}K_{{14}} (\tilde{{\mathcal {R}}}^2_0-1)I_{Hu}\\&\quad +{\frac{{\varphi ^{2}_{{Z}}} \left( N_{{{ Hp}}}Q_{{u}}R_{{2}}S_{{{ Hu}}}+N_{{{ Hu}}}Q_{{p}}R_{{1}}S_{{{ Hp}}} \right) {\kappa ^{2}_{ {V}}}\beta _{{V}} \left( C_{{3}} \left( K_{{9}}+K_{{12}} \right) +K _{{9}}K_{{11}} \right) ^{2}{\theta ^{2}_{{Y}}}S_{{X}}\gamma _{{H}}\beta _ {{M}}}{b_HK_{{12}}C_{{3}} \left( C_{{1}}K_{{9}}K_{{11}}-C_{{2}}\theta _{{Y }}\varphi _{{Z}} \right) }}\\&\quad \left( 1-\frac{1}{\tilde{{\mathcal {R}}}^2_0}\right) E_{Y}. \end{aligned} \end{aligned}$$

Hence, ${\dot{\mathcal {F}}}\le 0$ if $\tilde{{\mathcal {R}}}_0<1$ with ${\dot{\mathcal {F}}}= 0$ if and only if $I_X=I_{Hp}=I_{Hu}=E_Y=0.$ Therefore, $\mathcal {{F}}$ is a Lyapunov function in $\Omega$ and it follows from LaSalle’s Invariance Principle (1976) that every solution to the equations in (2.1), (2.2), (2.4) (with $\delta _H=0$ and initial conditions in $\Omega$) converges to ${\mathcal {T}}_2$ as $t\rightarrow \infty$ That is,

$$\begin{aligned}&(E_X(t),I_X(t),E_Y(t),I_Y(t),E_Z(t),I_Z(t),E_{Hp}(t),I_{Hp}(t),R_{Hp}(t),\\&\quad E_{Hu}(t),I_{Hu}(t),R_{Hu}(t))\\&\quad \rightarrow (0,0,0,0,0,0,0,0,0,0,0,0)\,\,\mathrm{as}\,\, \mathrm{t}\rightarrow \infty .\\&\mathrm{Thus},\,\,{\mathcal {X}}(t) \rightarrow \left( E^*,L^*_1,L^*_2,L^*_3,L^*_4,P^*,S^*_X,0,0,S^*_Y,0,0,S^*_Z,0,0,\right. \\&\qquad \left. \frac{\Pi \,\pi _p}{\mu _H},0,0,0,\frac{\Pi \,\pi _u}{\mu _H},0,0,0\right) \end{aligned}$$

as $t\rightarrow \infty$ for $\tilde{{\mathcal {R}}}_0\le 1$. Hence, the NDFE, ${\mathcal {T}}_2$, is globally-asymptotically stable in $\Omega$ if $\tilde{{\mathcal {R}}}_0\le 1$ for the special case of the model (2.1), (2.2), (2.4) with $\delta _H=0$. This completes the proof. $\square$