A mathematical model of avian influenza with half-saturated incidence

Abstract

The widespread impact of avian influenza viruses not only poses risks to birds, but also to humans. The viruses spread from birds to humans and from human to human In addition, mutation in the primary strain will increase the infectiousness of avian influenza. We developed a mathematical model of avian influenza for both bird and human populations. The effect of half-saturated incidence on transmission dynamics of the disease is investigated. The half-saturation constants determine the levels at which birds and humans contract avian influenza. To prevent the spread of avian influenza, the associated half-saturation constants must be increased, especially the half-saturation constant H _m for humans with mutant strain. The quantity H _m plays an essential role in determining the basic reproduction number of this model. Furthermore, by decreasing the rate β _m at which human-to-human mutant influenza is contracted, an outbreak can be controlled more effectively. To combat the outbreak, we propose both pharmaceutical (vaccination) and non-pharmaceutical (personal protection and isolation) control methods to reduce the transmission of avian influenza. Vaccination and personal protection will decrease β _m, while isolation will increase H _m. Numerical simulations demonstrate that all proposed control strategies will lead to disease eradication; however, if we only employ vaccination, it will require slightly longer to eradicate the disease than only applying non-pharmaceutical or a combination of pharmaceutical and non-pharmaceutical control methods. In conclusion, it is important to adopt a combination of control methods to fight an avian influenza outbreak.

Appendices

Appendix 1: Feasibility of the solution

Since the model parameters are non-negative, it is important to show that all state variables remain non-negative for all non-negative initial conditions for t ≥ 0. From Eq. (2.1), we have

$$ \begin{aligned} \frac{{\rm d}N_{\rm b}}{{\rm d}t}&=\Uplambda_{\rm b}-\mu_{\rm b} N_{\rm b}-\delta_{\rm b} I_{\rm b} \leq \Uplambda_{\rm b}-\mu_{\rm b} N_{\rm b}\\ \frac{{\rm d}N_{\rm h}}{{\rm d}t}&=\Uplambda_{\rm h}-\mu_{\rm h}N_{\rm h}-{\rm d}I_{\rm a}-\alpha I_{\rm m} \leq \Uplambda_{\rm h}-\mu_{\rm h}N_{\rm h}. \\ \end{aligned} $$

The closed set

$$ {\mathcal{D}}=\left\{(S_{\rm b},I_{\rm b},S_{\rm h},I_{\rm a},I_{\rm m})\in{\mathbb{R}}_+ ^5:N_{\rm b}\leq \frac{\Uplambda_{\rm b}}{\mu_{\rm b}},N_{\rm h}\leq\frac{\Uplambda_{\rm h}}{\mu_{\rm h}}\right\} $$

is a feasible region of the model.

Proposition 1

The closed set ${\mathcal{D}}$ is bounded and positively invariant.

Proof

Because $\frac{{\rm d}N_{\rm b}}{{\rm d}t}\leq \Uplambda_{\rm b}-\mu_{\rm b} N_{\rm b}, \, N_{\rm b}$ is bounded above by $\frac{\Uplambda_{\rm b}}{\mu_{\rm b}}.$ Hence $\frac{{\rm d}N_{\rm b}}{{\rm d}t}<0$ whenever $N_{\rm b}(t)>\frac{\Uplambda_{\rm b}}{\mu_{\rm b}}.$ Using an integrating factor, we have

$$ N_{\rm b}(t)\leq N_{\rm b}(0){\rm e}^{-\mu_{\rm b}t}+\frac{\Uplambda_{\rm b}}{\mu_{\rm b}}(1-{\rm e}^{-\mu_{\rm b}t}). $$

As $t \rightarrow \infty,\; {\rm e}^{-\mu_{\rm b}t}\rightarrow 0$ and hence $\lim_{t \to \infty}N_{\rm b}(t)\leq \frac{\Uplambda_{\rm b}}{\mu_{\rm b}}.$

The other case is similar. Thus ${\mathcal{D}}$ is bounded and positively invariant in $\mathbb{R}_+ ^5.$

Appendix 2: Stability analysis of the avian-only model

We consider the avian-only model, given by the first two equations of the Eq. (2.1)

A feasible region is defined as

$${\mathcal{D}}_{\rm b}=\left\{(S_{\rm b},I_{\rm b})\in {\mathbb{R}}_+^2:S_{\rm b}+I_{\rm b}\leq\frac{\Uplambda_{\rm b}}{\mu_{\rm b}}\right\}. $$

It can be shown from Proposition 1 that ${\mathcal{D}}_{\rm b}$ is bounded and positively invariant over $\mathbb{R} _+^2.$ The DFE (disease-free equilibrium) is

$$ E_{\rm b}^0=\left(S_{\rm b}^0, I_{\rm b}^0 \right)=\left(\frac{\Uplambda_{\rm b}}{\mu_{\rm b}},0 \right) $$

and the EE (endemic equilibrium) is

$$ E_{\rm b}^*=(S_{\rm b}^*,I_{\rm b}^*)=\left(\frac{\Uplambda_{\rm b}+H_{\rm b}(\mu_{\rm b}+\delta_{\rm b})}{\mu_{\rm b}+\beta_{\rm b}}, \frac{\beta_{\rm b}\Uplambda_{\rm b}-\mu_{\rm b}H_{\rm b}(\mu_{\rm b}+\delta_{\rm b})}{(\mu_{\rm b}+\beta_{\rm b}) (\mu_{\rm b}+\delta_{\rm b})} \right). $$

The basic reproduction number (see Li et al. 2011; van den Driessche and Watmough 2002 for further details and some complications) for the avian-only model is thus

$$ R_{\rm b} = \frac{\beta_{\rm b}\Uplambda_{\rm b}}{\mu_{\rm b}H_{\rm b}(\mu_{\rm b}+\delta_{\rm b})}. $$

Next, we would like to determine whether or not the model achieves global stability of the DFE or EE with respect to positive initial conditions.

Theorem 2

(Global stability of the DFE for the avian-only model) Let $E_{\rm b}^{0}=(S_{\rm b}^{0},I_{\rm b}^{0})=\left(\frac{\Uplambda_{\rm b}}{\mu_{\rm b}}, 0\right).$ If R _b < 1, then the DFE, E ⁰_b , is globally stable in the interior $\Upgamma_{\rm b}=\{(S_{\rm b},I_{\rm b})\in \mathbb{R}^2_{+}: S_{\rm b}, I_{\rm b}\le N_{\rm b}\}.$

Proof

Consider a Lyapunov function, $L(S_{\rm b},I_{\rm b})=\frac{\Uplambda_{\rm b}}{\mu_{\rm b}H_{\rm b}(\mu_{\rm b}+\delta_{\rm b})}I_{\rm b}.$ At the DFE, we have $L(S_{\rm b}^{0},I_{\rm b}^{0})=L\left(\frac{\Uplambda_{\rm b}}{\mu_{\rm b}}, 0\right)=0.$ Its derivative is

$$ \begin{aligned} \frac{{\rm d}L}{{\rm d}t}&= \frac{\Uplambda_{\rm b}}{\mu_{\rm b}H_{\rm b}(\mu_{\rm b}+\delta_{\rm b})} \left[\frac{\beta_{\rm b}S_{\rm b}I_{\rm b}}{H_{\rm b}+I_{\rm b}}-(\mu_{\rm b}+\delta_{\rm b})I_{\rm b} \right] \\ &= R_{\rm b}S_{\rm b}\frac{I_{\rm b}}{H_{\rm b}+I_{\rm b}}-\frac{\Uplambda_{\rm b}I_{\rm b}}{\mu_{\rm b}H_{\rm b}} \\ & \le R_{\rm b}S_{\rm b} \frac{I_{\rm b}}{H_{\rm b}}-\frac{\Uplambda_{\rm b}I_{\rm b}}{\mu_{\rm b}H_{\rm b}} \\ & \le R_{\rm b}\frac{\Uplambda_{\rm b}}{\mu_{\rm b}} \left(\frac{I_{\rm b}}{H_{\rm b}}\right)-\frac{\Uplambda_{\rm b}I_{\rm b}}{\mu_{\rm b}H_{\rm b}} \quad {\hbox {where at the DFE, we have}} \; N_{\rm b} \le \frac{\Uplambda_{\rm b}}{\mu_b} \Rightarrow S_{\rm b} \le \frac{\Uplambda_{\rm b}}{\mu_{\rm b}} \\ &= \frac{\Uplambda_{\rm b}I_{\rm b}}{\mu_{\rm b}H_{\rm b}}(R_{\rm b}-1) \\ &< 0 \quad {\hbox {if }} R_{\rm b}<1 \end{aligned} $$

Thus a periodic solution for this avian-only model does not exist for $(S_{\rm b}, I_{\rm b}) \in \Upgamma_{\rm b}.$ Therefore, the global stability of the DFE is satisfied.□

Theorem 3

(Global stability of the EE for the avian-only model) Let $E_{\rm b}^*=(S_{\rm b}^*,I_{\rm b}^*)$. If R _b > 1, then the EE, E ^*_b , is globally stable in ${\Upgamma_{\rm b} = \{(S_{\rm b},I_{\rm b})\in \mathbb{R}_{+}^{2}: S_{\rm b}\le N_{\rm b}, I_{\rm b}\le N_{\rm b}\}. }$

Proof

Let $f_{1}=\Uplambda_{\rm b}-\mu_{\rm b}S_{\rm b}-\frac{\beta_{\rm b} S_{\rm b} I_{\rm b}}{H_{\rm b}+I_{\rm b}},\,f_{2}=\frac{\beta_{\rm b} S_{\rm b}I_{\rm b}}{H_{\rm b}+I_{\rm b}}-(\mu_{\rm b}+\delta_{\rm b})I_{\rm b}$ and $B(S_{\rm b},I_{\rm b})=\frac{1}{I_{\rm b}}.$

$$\nabla(Bf)=\frac{\partial}{\partial S_{\rm b}} \left(Bf_{1}\right)+\frac{\partial}{\partial I_{\rm b}}(Bf_{2})= - \left[\frac{\mu_{\rm b}}{I_{\rm b}}+\frac{\beta_{\rm b}}{H_{\rm b}+I_{\rm b}}+\frac{\beta_{\rm b}S_{\rm b}}{(H_{\rm b}+I_{\rm b})^{2}} \right]\,< 0 \quad \forall\ (S_{\rm b}, I_{\rm b}) \in \Upgamma_{\rm b}.$$

As $\nabla(Bf)<0 \forall\ (S_{\rm b}, I_{\rm b}) \in \Upgamma_{\rm b},$ then by Bendixson’s Negative Criterion (Perko 2001), no periodic orbit can lie entirely in $\Upgamma_{\rm b}.$ Since R _b > 1, the DFE is unstable and hence, in a two-dimensional system, the EE is globally stable in $\Upgamma_{\rm b}.$□

Appendix 3: Stability analysis of the avian–human model

Since R _h decouples from the remaining equations in model (2.1), we consider the first five equations of model (2.1). We denote this as the avian–human model.

The transmission matrix, F, and transition matrix, V, of this model are

$$ \begin{aligned} F&=\left( \begin{array}{llll} \frac{\beta_{\rm b}H_{\rm b}S_{\rm b}}{(H_{\rm b}+I_{\rm b})^2} & 0 & 0 \\ \frac{\beta_{\rm bh}H_{\rm hb}S_{\rm h}}{(H_{\rm bh}+I_{\rm b})^2} & \frac{\beta_{\rm a}H_{\rm a}S_{\rm h}}{(H_{\rm a}+I_{\rm a})^2} & 0\\ 0 & 0 & \frac{\beta_{\rm m}H_{\rm m}S_{\rm h}}{(H_{\rm m}+I_{\rm m})^2}\\ \end{array} \right)\\ V&=\left( \begin{array}{lll} \mu_{\rm b}+\delta_{\rm b} & 0 & 0 \\ 0 & \mu_{\rm h}+d+\epsilon+\gamma_{\rm a} & 0\\ 0 & -\epsilon &\mu_{\rm h}+\alpha+\gamma_{\rm m}\\ \end{array} \right)\end{aligned}. $$

Thus we have

$$ FV^{-1}=\left( \begin{array}{llllll} \frac{\beta_{\rm b}H_{\rm b}S_{\rm b}}{(\mu_{\rm b}+\delta_{\rm b}) (H_{\rm b}+I_{\rm b})^2} & 0 & 0 \\ \frac{\beta_{\rm bh}H_{\rm hb}S_{\rm h}}{(\mu_{\rm b}+\delta_{\rm b}) (H_{\rm bh}+I_{\rm b})^2} & \frac{\beta_{\rm a}H_{\rm a}S_{\rm h}}{(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a}) (H_{\rm a}+I_{\rm a})^2} & 0\\ 0 & \frac{\epsilon\beta_{\rm m}H_{\rm m}S_{\rm h}}{(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a}) (\mu_{\rm h}+\alpha+\gamma_{\rm m}) (H_{\rm m}+I_{\rm m})^2} & \frac{\beta_{\rm m}H_{\rm m}S_{\rm h}}{(\mu_{\rm h}+\alpha+\gamma_{\rm m}) (H_{\rm m}+I_{\rm m})^2}\\ \end{array} \right). $$

The DFE of this model is

$$ E_{\rm ah}^0=(S_{\rm b}^0, I_{\rm b}^0, S_{\rm h}^0, I_{\rm a}^0, I_{\rm m}^0)=\left(\frac{\Uplambda_{\rm b}}{\mu_{\rm b}}, 0, \frac{\Uplambda_{\rm h}}{\mu_{\rm h}}, 0, 0\right). $$

At the DFE, we have

$$ FV^{-1}=\left( \begin{array}{llllll} \frac{\beta_{\rm b}\Uplambda_{\rm b}}{H_{\rm b}\mu_{\rm b}(\mu_{\rm b}+\delta_{\rm b})} & 0 & 0 \\ \frac{\beta_{\rm bh}\Uplambda_{\rm h}}{H_{\rm bh}\mu_{\rm h}(\mu_{\rm b}+\delta_{\rm b})} & \frac{\beta_{\rm a}\Uplambda_{\rm h}}{H_{\rm a}\mu_{\rm h}(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a})} & 0\\ 0 & \frac{\epsilon\beta_{\rm m}\Uplambda_{\rm h}}{H_{\rm m}\mu_{\rm h}(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a}) (\mu_{\rm h}+\alpha+\gamma_{\rm m})} & \frac{\beta_{\rm m}\Uplambda_{\rm h}}{H_{\rm m}\mu_{\rm h}(\mu_{\rm h}+\alpha+\gamma_{\rm m})}\\ \end{array} \right) $$

and

$$ \begin{aligned} R_{\rm ah} &\equiv \max\left \{\frac{\beta_{\rm b}\Uplambda_{\rm b}}{H_{\rm b}\mu_{\rm b}(\mu_{\rm b}+\delta_{\rm b})} , \frac{\beta_{\rm a}\Uplambda_{\rm h}}{H_{\rm a}\mu_{\rm h}(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a})}, \frac{\beta_{\rm m}\Uplambda_{\rm h}}{H_{\rm m}\mu_{\rm h}(\mu_{\rm h}+\alpha+\gamma_{\rm m})}\right \} \\ \Rightarrow R_{\rm ah}&=\max\{R_{\rm b}, R_{{\rm h}1}, R_{{\rm h}2}\} \end{aligned}$$

(5.3)

where $R_{{\rm h}1}=\frac{\beta_{\rm a}\Uplambda_{\rm h}}{H_{\rm a}\mu_{\rm h}(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a})} \quad \hbox {and} \; \quad R_{{\rm h}2}=\frac{\beta_{\rm m}\Uplambda_{\rm h}}{H_{\rm m}\mu_{\rm h}(\mu_{\rm h}+\alpha+\gamma_{\rm m})}.$

Theorem 4

(Global stability of the DFE for the avian–human model) Let $E_{\rm ah}^0=(S_{\rm b}^0, I_{\rm b}^0, S_{\rm h}^0, I_{\rm a}^0, I_{\rm m}^0)=\left (\frac{\Uplambda_{\rm b}}{\mu_{\rm b}}, 0, \frac{\Uplambda_{\rm h}}{\mu_{\rm h}}, 0, 0\right).$ If R _ah < 1, then the DFE, E ⁰_ah , is globally stable in the interior ${\Upgamma_{\rm ah}=\{(S_{\rm b},I_{\rm b},S_{\rm h},I_{\rm a},I_{\rm m})\in \mathbb{R}^5_{+}: S_{\rm b}, I_{\rm b}\le N_{\rm b}, S_{\rm h}, I_{\rm a}, I_{\rm m}\le N_{\rm h}\}. }$

Proof

Consider a Lyapunov function,

$$ L(S_{\rm b},I_{\rm b},S_{\rm h},I_{\rm a},I_{\rm m})=\frac{\Uplambda_{\rm b}}{H_{\rm b}\mu_{\rm b}(\mu_{\rm b}+\delta_{\rm b})}I_{\rm b}+\frac{\Uplambda_{\rm h}}{H_{\rm a}\mu_{\rm h}(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a})} I_{\rm a} + \frac{\Uplambda_{\rm h}}{H_{\rm m}\mu_{\rm h}(\mu_{\rm h}+\alpha+\gamma_{\rm m})} I_{\rm m}. $$

At the DFE, we have $L(S_{\rm b}^0,I_{\rm b}^0,S_{\rm h}^{0},I_{\rm a}^0,I_{\rm m}^0)=L\left(\frac{\Uplambda_{\rm b}}{\mu_{\rm b}},0,\frac{\Uplambda_{\rm h}}{\mu_{\rm h}}, 0, 0\right)=0.$ Its derivative is

$$ \begin{aligned} \frac{{\rm d}L}{{\rm d}t}&= \frac{\Uplambda_{\rm b}}{H_{\rm b}\mu_{\rm b}(\mu_{\rm b}+\delta_{\rm b})}I'_{\rm b}+\frac{\Uplambda_{\rm h}}{H_{\rm a}\mu_{\rm h}(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a})} I'_{\rm a} + \frac{\Uplambda_{\rm h}}{H_{\rm m}\mu_{\rm h}(\mu_{\rm h}+\alpha+\gamma_{\rm m})} I'_{\rm m}\\ &=\frac{\Uplambda_{\rm b}}{H_{\rm b}\mu_{\rm b}(\mu_{\rm b}+\delta_{\rm b})} \left [\frac{\beta_{\rm b}S_{\rm b}I_{\rm b}}{H_{\rm b}+I_{\rm b}}-(\mu_{\rm b}+\delta_{\rm b})I_{\rm b} \right]+\frac{\Uplambda_{\rm h}}{H_{\rm a}\mu_{\rm h}(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a})} \Bigg[\frac{\beta_{\rm bh}S_{\rm h}I_{\rm b}}{H_{\rm bh}+I_{\rm b}}+\frac{\beta_{\rm a} S_{\rm h} I_{\rm a}}{H_{\rm a}+I_{\rm a}}\\ &\quad -(\mu_{\rm h}+d+\epsilon+\gamma_{\rm a})I_{\rm a} \Bigg]+\frac{\Uplambda_{\rm h}}{H_{\rm m}\mu_{\rm h}(\mu_{\rm h}+\alpha+\gamma_{\rm m})} \left[\frac{\beta_{\rm m} S_{\rm h} I_{\rm m}}{H_{\rm m}+I_{\rm m}}+\epsilon I_{\rm a}-(\mu_{\rm h}+\alpha+\gamma_{\rm m})I_{\rm m} \right]\\ & \le R_{\rm b}S_{\rm b}\frac{I_{\rm b}}{H_{\rm b}+I_{\rm b}}-\frac{\Uplambda_{\rm b}I_{\rm b}}{H_{\rm b}\mu_{\rm b}}+R_{{\rm h}1}S_{\rm h}\frac{I_{\rm a}}{H_{\rm a}+I_{\rm a}}-\frac{\Uplambda_{\rm h}I_{\rm a}}{H_{\rm a}\mu_{\rm h}}+R_{{\rm h}2}S_{\rm h}\frac{I_{\rm m}}{H_{\rm m}+I_{\rm m}}-\frac{\Uplambda_{\rm h}I_{\rm m}}{H_{\rm m}\mu_{\rm h}}\\ \le R_{\rm b}\frac{\Uplambda_{\rm b}I_{\rm b}}{\mu_{\rm b}H_{\rm b}}-\frac{\Uplambda_{\rm b}I_{\rm b}}{H_{\rm b}\mu_{\rm b}}+R_{{\rm h}1} \frac{\Uplambda_{\rm h}I_{\rm a}}{\mu_{\rm h}H_{\rm a}}-\frac{\Uplambda_{\rm h}I_{\rm a}}{H_{\rm a}\mu_{\rm h}}+R_{{\rm h}2} \frac{\Uplambda_{\rm h}I_m}{\mu_{\rm h}H_{\rm m}}-\frac{\Uplambda_{\rm h}I_{\rm m}}{H_{\rm m}\mu_{\rm h}}\\&= \frac{\Uplambda_{\rm b}I_{\rm b}}{\mu_{\rm b}H_b}(R_{\rm b}-1)+\frac{\Uplambda_{\rm h}I_{\rm a}}{H_{\rm a}\mu_{\rm h}} (R_{{\rm h}1}-1)+\frac{\Uplambda_{\rm h}I_m}{H_{\rm m}\mu_{\rm h}} (R_{{\rm h}2}-1)\\ &< 0 \quad \hbox {if } \;R_{\rm b}, R_{{\rm h}1}, R_{{\rm h}2}<1 \end{aligned} $$

Thus, a periodic solution for this avian–human model does not exist for $(S_{\rm b},I_{\rm b},S_{\rm h}, I_{\rm a}, I_{\rm m}) \in \Upgamma_{\rm ah}.$□