Dynamical analysis of a fractional-order avian-human influenza epidemic model with logistic growth for avian population

Abstract

In this paper, a fractional order avian-human influenza epidemic model with logistic growth for avian population is investigated. The dynamical behavior of this model is discussed. We first establish the existence, uniqueness, non-negativity and positive invariance of the solution. Then we analyze the existence of various equilibrium points, and some sufficient conditions are derived to ensure the global asymptotic stability of the disease free equilibrium point and endemic equilibrium point. Finally, we take some numerical simulations to validate the analytical results.

Introduction

Avian influenza is an infectious disease of birds caused by type A strains of the influenza virus. Most avian influenza viruses do not cause disease in humans. However, some can infect humans and cause disease. Influenza A virus subtype H5N1 is a highly pathogenic avian influenza virus, which was first detected in 1996 in geese in China and first detected in humans in 1997 during a poultry outbreak in Hong Kong. ¹ Differing from H5N1, the H7N9 virus is classified as a low pathogenicity avian influenza virus. ² However, the virus can infect humans and most of the reported cases of human H7N9 infection have resulted in severe respiratory illness. ³ So it is important for us to take effective measure to prevent the avian influenza and limit its damage at the minimum.

Mathematical modeling represents an important tool for analyzing the epidemiological characteristics of infectious diseases and can provide useful control measures. Various of mathematical compartmental models have been introduced to study the mechanism of the spread of avian influenza.^1,4–8 However, most of these studies focused on an integer-order compartmental epidemic model as a key tool for their investigation. In recent years, it has frequently been observed in many areas of engineering, physics, life sciences, and elsewhere that models based on fractional order derivatives can provide better agreement between measured and simulated data than classical models based on integer order derivatives^9,10 and the references therein.

In fact, fractional-order differential equations are naturally related to systems with memory which exists in most biological systems, and they have close relations to fractals which are abundant in biological systems. In recent past, fractional order differential equations have been used in several biological systems to explore the underlying dynamics. In Diethelm, ¹¹ Diethelm proposed a nonlinear fractional order differential equation model for the simulation of the dynamics of a dengue fever outbreak. Sardar et al. ¹² proposed two mathematical models to study the spread of a common single strain mosquito-transmitted diseases. Some fractional order models for HIV infection were studied in literature.^13,14 In Gonzalez-Parra et al., ¹⁵ the authors proposed a SEIR epidemic fractional differential model to explain and understand the outbreaks of influenza A(H1N1). Li et al. ¹⁶ investigated a fractional-order single-species model, which is composed of $n(n>2)$ patches connected by diffusion. A homogeneous-mixing population fractional model for human immunodeficiency virus (HIV) transmission was proposed in Huo et al. ¹⁷ A nonlinear fractional order epidemic model was proposed in Ye et al. ¹⁸ to investigate the spreading dynamical behavior of the avian influenza. Also we refer the interested reader in fractional biological systems^19–22 and the references therein for some other recent works in the subject. However, to the best of our knowledge, to this day, less work has investigated the dynamic behavior of fractional-order avian-human influenza model.

In this paper, we consider the following fractional-order avian-human influenza model:

$$\begin{array}{l}\frac{{\text{d}}^{\alpha }{S}_a}{\text{d}{t}^{\alpha }}=r_a{S}_a(1-\frac{S_a}{K_a})-{\beta }_a{S}_a{I}_a,\\ \frac{{\text{d}}^{\alpha }{I}_a}{\text{d}{t}^{\alpha }}={\beta }_a{S}_a{I}_a-({\mu }_a+{\delta }_a)I_a,\\ \frac{{\text{d}}^{\alpha }{S}_h}{\text{d}{t}^{\alpha }}={\Pi }_h-{\beta }_h{I}_a{S}_h-{\mu }_h{S}_h,\\ \frac{{\text{d}}^{\alpha }{I}_h}{\text{d}{t}^{\alpha }}={\beta }_h{I}_a{S}_h-({\mu }_h+{\delta }_h+\gamma )I_h,\\ \frac{{\text{d}}^{\alpha }{R}_h}{\text{d}{t}^{\alpha }}=\gamma I_h-{\mu }_h{R}_h\end{array}$$

(1)

with $\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}$ denoting the left-sided Caputo fractional derivative of order $\alpha \in (0,1],$ following Podlubny, ²³ defined as

$$\frac{{\text{d}}^{\alpha }y(t)}{\text{d}{t}^{\alpha }}=\frac{1}{\Gamma (1-\alpha )}{\displaystyle {\int }_0^t\frac{y\prime (\tau )\text{d}\tau }{{(t-\tau )}^{\alpha }}}$$

where $\Gamma (·)$ is the Gamma function. In model (1), $S_a(t)$ and $I_a(t)$ are the population densities or numbers of susceptible and infective avian population, respectively; $S_h(t),I_h(t)$ and $R_h(t)$ are the population densities or numbers of susceptible, infective and recovered/removed human population, respectively; r_a and K_a are the intrinsic growth rate and maximal carrying capacity of the avian population, respectively; β_a is the transmission rate from infective avian to susceptible avian; μ_a and δ_a are the natural death rate and disease death rate of the avian population, respectively. Π _h represents all new recruitments and newborns of the human population; β_h is the transmission rate from the infective avian to the susceptible human, μ_h is the natural death rate of the human population; δ_h is the disease related rate of the infected human population; γ is the recovery rate of the infective human. All parameters are assumed to be positive from biological point of view. The case of integer order model, i.e., the counterpart with α = 1 of (1) has been investigated by Liu et al. ²⁴

A simple dimensional analysis of model (1) shows that the left-hand side of the system (1) has dimension (time) ${\hspace{0.17em}}^{-\alpha }$ whereas the right-hand side dimension is (time)^–1. This dimension mismatch of the system (1) can be mathematically corrected using the procedure described in Diethelm. ¹¹ A corrected form corresponding to model (1) is written as

$$\begin{array}{l}\frac{{\text{d}}^{\alpha }{S}_a}{\text{d}{t}^{\alpha }}=r_a^{\alpha }{S}_a(1-\frac{S_a}{K_a})-{\beta }_a^{\alpha }{S}_a{I}_a,\\ \frac{{\text{d}}^{\alpha }{I}_a}{\text{d}{t}^{\alpha }}={\beta }_a^{\alpha }{S}_a{I}_a-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a,\\ \frac{{\text{d}}^{\alpha }{S}_h}{\text{d}{t}^{\alpha }}={\Pi }_h^{\alpha }-{\beta }_h^{\alpha }{I}_a{S}_h-{\mu }_h^{\alpha }{S}_h,\\ \frac{{\text{d}}^{\alpha }{I}_h}{\text{d}{t}^{\alpha }}={\beta }_h^{\alpha }{I}_a{S}_h-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h,\\ \frac{{\text{d}}^{\alpha }{R}_h}{\text{d}{t}^{\alpha }}={\gamma }^{\alpha }{I}_h-{\mu }_h^{\alpha }{R}_h\end{array}$$

(2)

In this paper we will investigate the above system, subject to an initial condition given at t = t₀. The remaining of this paper is organized as follows. In the next section, we derive the well-posedness of the problem. The dynamical behavior of the model are presented in ‘Dynamical behavior’ section. Finally, we present some numerical simulations in the penultimate section to verify our theoretical results. Some concluding remarks are given at the end of the paper.

Well-posedness

In this section, we will prove the existence, uniqueness, non-negativity and boundedness of the solutions. We first introduce some basic results.

Preliminaries

Lemma 2.1 (Generalized mean value theorem ²⁵ ) Suppose that $f(t)\in C[a,b]$ and $\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}f(t)\in C[a,b]$ with $0<\alpha \le 1,$ then we have

$$f(t)=f(a)+\frac{1}{\Gamma (\alpha )}\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}f(\xi ){(t-a)}^{\alpha }$$

where $a\le \xi \le t,\forall t\in [a,b].$

Lemma 2.2 ²⁰ Let u(t) be a continuous function on $[t_0,\infty ]$ and satisfy $\frac{{\text{d}}^{\alpha }u(t)}{\text{d}{t}^{\alpha }}\le -\lambda u(t)+\mu u(t_0)=u_{t_0}$ where $0<\alpha <1,(\lambda ,\mu )\in {ℝ}^2,\lambda \ne 0$ and $t_0\ge 0$ is the initial time. Then, its solution has the form

$$u(t)\le \left(u_{t_0}-\frac{\mu }{\lambda }\right)E_{\alpha }[-\lambda {(t-t_0)}^{\alpha }]+\frac{\mu }{\lambda }$$

$E_{\alpha }(z)$ is the Mittag–Leffler function with parameter α.

Lemma 2.3 ¹⁹ Let $x(t)\in {ℝ}^{+}$ be a continuous and derivable function. Then, for any time instant $t\ge t_0$ , $\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}\left[x(t)-x^{\ast }-x^{\ast }\text{ln}\frac{x(t)}{x^{\ast }}\right]\le \left(1-\frac{x^{\ast }}{x(t)}\right)\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}x(t)$

where $x^{\ast }\in {ℝ}^{+},\forall \alpha \in (0,1).$

Existence and uniqueness

Theorem 2.1 The system (2) possesses a unique continuous solution on $[0,T]\times \Omega$ , where $T<+\infty ,\Omega =\{(S_a,I_a,S_h,I_h,R_h)\in {ℝ}^5:\max \{\left|S_a\right|,\left|I_a\right|,\left|S_h\right|,\left|I_h\right|,\left|R_h\right|\}\le M\}$ and M is sufficiently large.

Proof. We denote $X=(S_a,I_a,S_h,I_h,R_h)$ and $\overline{X}=(\overline{S_a},\overline{I_a},\overline{S_h},\overline{I_h},\overline{R_h})$ . Consider a mapping $f(X)=(f_1(X),f_2(X),f_3(X),f_4(X),f_5(X))$ and

$$\begin{array}{l}{f}_1(X)=r_a^{\alpha }{S}_a(1-\frac{S_a}{K_a})-{\beta }_a^{\alpha }{S}_a{I}_a,\\ f_2(X)={\beta }_a^{\alpha }{S}_a{I}_a-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a,\\ f_3(X)={\Pi }_h^{\alpha }-{\beta }_h^{\alpha }{I}_a{S}_h-{\mu }_h^{\alpha }{S}_h,\\ f_4(X)={\beta }_h^{\alpha }{I}_a{S}_h-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h,\\ f_5(X)={\gamma }^{\alpha }{I}_h-{\mu }_h^{\alpha }{R}_h\end{array}$$

(3)

For any $X,\overline{X}\in \Omega$ , it follows from (3) that

$$\begin{array}{l}\parallel f(X)-f(\overline{X})||=\left|f_1(X)-f_1(\overline{X})\right|\\ +\left|f_2(X)-f_2(\overline{X})\right|+\left|f_3(X)-f_3(\overline{X})\right|\\ +\left|f_4(X)-f_4(\overline{X})\right|+\left|f_5(X)-f_5(\overline{X})\right|\\ =\left|r_a^{\alpha }{S}_a\right(1-\frac{S_a}{K_a})-{\beta }_a^{\alpha }{S}_a{I}_a-r_a^{\alpha }\overline{S_a}(1-\frac{\overline{S_a}}{K_a})\\ +{\beta }_a^{\alpha }\overline{S_a}\overline{I_a}|\\ +\left|{\beta }_a^{\alpha }{S}_a{I}_a-{\beta }_a^{\alpha }\overline{S_a}\overline{I_a}-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })(I_a-\overline{I_a})\right|\\ +\left|-{\beta }_h^{\alpha }{I}_a{S}_h+{\beta }_h^{\alpha }\overline{I_a}\overline{S_h}-{\mu }_h^{\alpha }(S_h-\overline{S_h})\right|\\ +\left|{\beta }_h^{\alpha }{I}_a{S}_h-{\beta }_h^{\alpha }\overline{I_a}\overline{S_h}-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })(I_h-\overline{I_h})\right|\\ +\left|{\gamma }^{\alpha }(I_h-\overline{I_h})-{\mu }_h^{\alpha }(R_h-\overline{R_h})\right|\\ =\left|r_a^{\alpha }(S_a-\overline{S_a})-\frac{r_a^{\alpha }}{K_a}(S_a^2-{\overline{S_a}}^2)-{\beta }_a^{\alpha }(S_a{I}_a-\overline{S_a}\overline{I_a})\right|\\ +\left|{\beta }_a^{\alpha }(S_a{I}_a-\overline{S_a}\overline{I_a})-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })(I_a-\overline{I_a})\right|\\ +\left|{\beta }_h^{\alpha }(I_a{S}_h-\overline{I_a}\overline{S_h})-{\mu }_h^{\alpha }(S_h-\overline{S_h})\right|\\ +\left|{\beta }_h^{\alpha }(I_a{S}_h-\overline{I_a}\overline{S_h})-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })(I_h-\overline{I_h})\right|\\ +\left|{\gamma }^{\alpha }(I_h-\overline{I_h})-{\mu }_h^{\alpha }(R_h-\overline{R_h})\right|\\ \le r_a^{\alpha }\left|S_a-\overline{S_a}\right|+\frac{r_a^{\alpha }}{K_a}\left|S_a^2-{\overline{S_a}}^2\right|\\ +2{\beta }_a^{\alpha }\left|S_a{I}_a-\overline{S_a}\overline{I_a}\right|+({\mu }_a^{\alpha }+{\delta }_a^{\alpha })\left|I_a-\overline{I_a}\right|\\ +2{\beta }_h^{\alpha }\left|I_a{S}_h-\overline{I_a}\overline{S_h}\right|+{\mu }_h^{\alpha }\left|S_h-\overline{S_h}\right|\\ +({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+2{\gamma }^{\alpha })\left|I_h-\overline{I_h}\right|\\ +{\mu }_h^{\alpha }\left|R_h-\overline{R_h}\right|\end{array}$$

(4)

Note that

$$\left|S_a^2-{\overline{S_a}}^2\right|=\left|S_a+\overline{S_a}\right|\left|S_a-\overline{S_a}\right|\le 2M\left|S_a-\overline{S_a}\right|$$

(5)

Furthermore, we obtain that

$$\begin{array}{l}\left|S_a{I}_a-\overline{S_a}\overline{I_a}\right|=\left|S_a{I}_a-\overline{S_a}{I}_a+\overline{S_a}{I}_a-\overline{S_a}\overline{I_a}\right|\\ \le \left|I_a\right|\left|S_a-\overline{S_a}\right|+\left|\overline{S_a}\right|\left|I_a-\overline{I_a}\right|\\ \le M\left|S_a-\overline{S_a}\right|+M\left|I_a-\overline{I_a}\right|\end{array}$$

(6)

Similarly, we get

$$\begin{array}{l}\left|I_a{S}_h-\overline{I_a}\overline{S_h}\right|=\left|I_a{S}_h-\overline{I_a}{S}_h+\overline{I_a}{S}_h-\overline{I_a}\overline{S_h}\right|\\ \le \left|S_h\right|\left|I_a-\overline{I_a}\right|+\left|\overline{I_a}\right|\left|S_h-\overline{S_h}\right|\\ \le M\left|I_a-\overline{I_a}\right|+M\left|S_h-\overline{S_h}\right|\end{array}$$

(7)

Plugging (5), (6) and (7) into (4), we have

$$\begin{array}{l}\parallel f(X)-f(\overline{X})||\le (r_a^{\alpha }+2M\frac{r_a^{\alpha }}{K_a}+2M{\beta }_a^{\alpha })\left|S_a-\overline{S_a}\right|\\ +({\mu }_a^{\alpha }+{\delta }_a^{\alpha }+2M{\beta }_a^{\alpha }+2M{\beta }_h^{\alpha })\left|I_a-\overline{I_a}\right|\\ +(2M{\beta }_h^{\alpha }+{\mu }_h^{\alpha })\left|S_h-\overline{S_h}\right|\\ +({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+2{\gamma }^{\alpha })\left|I_h-\overline{I_h}\right|+{\mu }_h^{\alpha }\left|R_h-\overline{R_h}\right|\\ \le L(\left|S_a-\overline{S_a}\right|+\left|I_a-\overline{I_a}\right|+\left|S_h-\overline{S_h}\right|\\ +\left|I_h-\overline{I_h}\right|+\left|R_h-\overline{R_h}\right|)\le L\parallel X-\overline{X}||\end{array}$$

(8)

where

$$\begin{array}{l}L=\max \{r_a^{\alpha }+2M(\frac{r_a^{\alpha }}{K_a}+{\beta }_a^{\alpha }),{\mu }_a^{\alpha }+{\delta }_a^{\alpha }\\ +2M({\beta }_a^{\alpha }+{\beta }_h^{\alpha }),{\mu }_h^{\alpha }+{\delta }_h^{\alpha }+2{\gamma }^{\alpha }\}\end{array}$$

Thus f(X) satisfies Lipschitz condition with respect to X, and it follows from the existence and uniqueness theorem (Theorem 3.4 in Podlubny ²³ ) that there exists a unique solution X(t) of the system (2) with the initial condition $X(0)=(S_a(0),I_a(0),S_h(0),I_h(0),R_h(0))$ . □

Nonnegativity and boundedness

Considering biological significance of the problem, we are only interested in solutions that are nonnegative and bounded. Denote ${ℝ}_{+}^5=\{x\in {ℝ}_{+}^5\}|x\ge 0\}$ and $x(t)={(S_a(t),I_a(t),S_h(t),I_h(t),R_h(t))}^T.$

Theorem 2.2 All solutions of the system (2) that start in ${ℝ}_{+}^5$ are nonnegative and uniformly bounded.

Proof. Firstly, we prove that the solutions $S_a(t)$ which start in ${ℝ}_{+}^5$ are non-negative. Assuming $S_a(0)>0$ for t = 0, we first prove that $S_a(t)\ge 0$ for all $t\ge 0$ . Let us suppose that $S_a(t)\ge 0,\forall t\ge 0$ is not true. Then, there exists a $t_1>0$ such that $S_a(t)>0$ for $0\le t<t_1,S_a(t)=0$ at t = t₁ and $S_a(t)<0$ for $t>t_1.$ From the first equation of (2), we have

$$\frac{{\text{d}}^{\alpha }{S}_a}{\text{d}{t}^{\alpha }}{|}_{t=t_1}=0$$

According to Lemma 2.1, we have $S_a(t_1^{+})=0$ which contradicts the fact $S_a(t_1^{+})<0$ , i.e., $S_a(t)<0$ for $t>t_1.$ Therefore, we have $S_a(t)\ge 0,t\ge 0$ . Using similar arguments, we can prove that $I_a(t),S_h(t),I_h(t)$ $R_h(t)$ are all non-negative for $t\ge 0$ . Next, we show that all solutions of system (2) that start in ${ℝ}_{+}^5$ are uniformly bounded.

Adding the first two equations of (2) and adding last three equations of (2), we get two fractional order differential equations representing the total avian and human population as follows:

$$\frac{{\text{d}}^{\alpha }{N}_a}{\text{d}{t}^{\alpha }}=r_a^{\alpha }{S}_a(1-\frac{S_a}{K_a})-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a$$

(9)

and

$$\frac{{\text{d}}^{\alpha }{N}_h}{\text{d}{t}^{\alpha }}={\prod }_h^{\alpha }-{\mu }_h^{\alpha }{N}_h-{\delta }_h^{\alpha }{I}_h$$

(10)

where $N_a=S_a+I_a$ and $N_h=S_h+I_h+R_h.$ Then from (9) we have

$$\begin{array}{l}\frac{{\text{d}}^{\alpha }{N}_a}{\text{d}{t}^{\alpha }}+({\mu }_a^{\alpha }+{\delta }_a^{\alpha })N_a\\ =r_a^{\alpha }{S}_a(1-\frac{S_a}{K_a})-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a+({\mu }_a^{\alpha }+{\delta }_a^{\alpha })(S_a+I_a)\\ =-\frac{r_a^{\alpha }}{K_a}{S}_a^2+(r_a^{\alpha }+{\mu }_a^{\alpha }+{\delta }_a^{\alpha })S_a\\ =-\frac{r_a^{\alpha }}{K_a}{\left(S_a-\frac{K_a(r_a^{\alpha }+{\mu }_a^{\alpha }+{\delta }_a^{\alpha })}{2r_a^{\alpha }}\right)}^2\\ +\frac{K_a{(r_a^{\alpha }+{\mu }_a^{\alpha }+{\delta }_a^{\alpha })}^2}{4r_a^{\alpha }}\le \frac{K_a{(r_a^{\alpha }+{\mu }_a^{\alpha }+{\delta }_a^{\alpha })}^2}{4r_a^{\alpha }}\end{array}$$

Applying Lemma 2.2, one gets

$$\begin{array}{l}{N}_a\left(t\right)\le \left(N_a\left(0\right)-\frac{K_a{(r_a^{\alpha }+{\mu }_a^{\alpha }+{\delta }_a^{\alpha })}^2}{4r_a^{\alpha }({\mu }_a^{\alpha }+{\delta }_a^{\alpha })}\right)\\ \times E_{\alpha }[-({\mu }_a^{\alpha }+{\delta }_a^{\alpha }\left)t^{\alpha }\right]+\frac{K_a{(r_a^{\alpha }+{\mu }_a^{\alpha }+{\delta }_a^{\alpha })}^2}{4r_a^{\alpha }({\mu }_a^{\alpha }+{\delta }_a^{\alpha })}\end{array}$$

Thus, $N_a\left(t\right)\to \frac{K_a{(r_a^{\alpha }+{\mu }_a^{\alpha }+{\delta }_a^{\alpha })}^2}{4r_a^{\alpha }({\mu }_a^{\alpha }+{\delta }_a^{\alpha })}$ as $t\to \infty$ .

Since $I_h\ge 0$ , we get from (10) that

$$\frac{{\text{d}}^{\alpha }{N}_h}{\text{d}{t}^{\alpha }}\le {\prod }_h^{\alpha }-{\mu }_h^{\alpha }{N}_h$$

Again by Lemma 2.2, we have

$$\begin{array}{l}{N}_H(t)\le \left(N_h(0)-\frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }}\right)E_{\alpha }(-{\mu }_h^{\alpha }{t}^{\alpha })\\ +\frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }}\to \frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }},t\to \infty \end{array}$$

(11)

Therefore, all the solutions of system (2) which start in ${ℝ}_{+}^5$ are confined to the region

$$\begin{array}{l}\Omega =\{(S_a,I_a,S_h,I_h,R_h)\in {ℝ}_{+}^5|S_a\\ +I_a\le \frac{K_a{(r_a^{\alpha }+{\mu }_a^{\alpha }+{\delta }_a^{\alpha })}^2}{4r_a^{\alpha }}+{\mathrm{ϵ}}_1,S_h\\ +I_h+R_h\le \frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }}+{\mathrm{ϵ}}_2\}\end{array}$$

with any ${\mathrm{ϵ}}_1>0,{\mathrm{ϵ}}_2>0$ . □

We will consider the stability of the equilibria of system (2) and its qualitative behaviors in the following sections.

Dynamical behavior

From (2), we can deduce two disease-free equilibria given by

$$E_1\left(0,0,\frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }},0,0\right)\text{and}{E}_2\left(K_a,0,\frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }},0,0\right)$$

Following the definition and computation procedure in Diekmann et al. ²⁶ and van den Driessche and Watmough, ²⁷ we can rewrite system (2) as follows:

$$\frac{{\text{d}}^{\alpha }X}{\text{d}{t}^{\alpha }}=\mathcal{F}-\mathcal{V}$$

where $X(t)={(I_a(t),I_h(t),S_a(t),S_h(t),R_h(t))}^T,$ and

$$\mathcal{F}=\left(\begin{array}{c}{\beta }_a^{\alpha }{S}_a{I}_a\\ {\beta }_h^{\alpha }{I}_a{S}_h\\ 0\\ 0\\ 0\end{array}\right),\mathcal{V}=\left(\begin{array}{c}({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a\\ ({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h\\ {\beta }_a^{\alpha }{S}_a{I}_a-r_a^{\alpha }{S}_a(1-\frac{S_a}{K_a})\\ {\beta }_h^{\alpha }{I}_a{S}_h-{\mu }_h^{\alpha }{S}_h-{\Pi }_h^{\alpha }\\ {\mu }_h^{\alpha }{R}_h-{\gamma }^{\alpha }{I}_h\end{array}\right)$$

Let

$$F_1=\left(\begin{array}{l}{\beta }_a^{\alpha }{S}_a{I}_a\\ {\beta }_h^{\alpha }{I}_a{S}_h\end{array}\right),\hspace{1em}{V}_1=\left(\begin{array}{c}({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a\\ ({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h\end{array}\right)$$

then

$$\begin{array}{l}F={\left(\frac{\partial F_1}{\partial I_a},\frac{\partial F_1}{\partial I_h},\right)|}_{E_2}=\left(\begin{array}{cc}{\beta }_a^{\alpha }{K}_a& 0\\ \frac{{\beta }_h^{\alpha }{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }}& 0\end{array}\right),\\ V={\left(\frac{\partial V_1}{\partial I_a},\frac{\partial V_1}{\partial I_h},\right)|}_{E_2}=\left(\begin{array}{cc}{\mu }_a^{\alpha }+{\delta }_a^{\alpha }& 0\\ 0& {\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha }\end{array}\right)\end{array}$$

and

$$F{V}^{-1}=\left(\begin{array}{cc}\frac{{\beta }_a^{\alpha }{K}_a}{{\mu }_a^{\alpha }+{\delta }_a^{\alpha }}& 0\\ \frac{{\beta }_h^{\alpha }{S}_h^{*}}{{\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha }}& 0\end{array}\right)$$

Hence, we derive the basic reproduction number as follows

$$R_0=\rho (F{V}^{-1})=\frac{{\beta }_a^{\alpha }{K}_a}{{\mu }_a^{\alpha }+{\delta }_a^{\alpha }}$$

(12)

where ρ denotes the spectral radius or largest eigenvalue of $F{V}^{-1}$ .

If $R_0>1$ , we can also derive a unique endemic equilibrium given by $E_3(S_a^{*},I_a^{*},S_h^{*},I_h^{*},R_h^{*})$ where

$$\begin{array}{l}{S}_a^{*}=\frac{{\mu }_a^{\alpha }+{\delta }_a^{\alpha }}{{\beta }_a^{\alpha }},I_a^{*}=\frac{r_a^{\alpha }}{{\beta }_a^{\alpha }}(1-\frac{1}{R_0}),\\ S_h^{*}=\frac{{\Pi }_h^{\alpha }}{{\beta }_h^{\alpha }{I}_a^{*}+{\mu }_h^{\alpha }},I_h^{*}=\frac{{\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}}{{\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha }},R_h^{*}=\frac{{\gamma }^{\alpha }{I}_h^{*}}{{\mu }_h^{\alpha }}\end{array}$$

(13)

Before analyzing the dynamical behavior of the full model (2), we study the dynamical behavior of the avian-only subsystem.

Analysis of the avian-only subsystem

In this subsection, we will consider the dynamics of fractional order avian-only subsystem, which is independent of the human system, as follows:

$$\begin{array}{l}\frac{{\text{d}}^{\alpha }{S}_a}{\text{d}{t}^{\alpha }}=r_a^{\alpha }{S}_a(1-\frac{S_a}{K_a})-{\beta }_a^{\alpha }{S}_a{I}_a,\\ \frac{{\text{d}}^{\alpha }{I}_a}{\text{d}{t}^{\alpha }}={\beta }_a^{\alpha }{S}_a{I}_a-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a\end{array}$$

(14)

The avian-only subsystem (14) always has two disease-free equilibria given by $E_{1_a}(0,0)$ and $E_{2_a}(K_a,0)$ . If $R_0>1$ , the system also has a unique endemic equilibrium given by $E_{3_a}(S_a^{*},I_a^{*})$ .

Lemma 3.1 (i) The disease-free equilibrium $E_{1_a}(0,0)$ of (14) is always unstable; (ii) If $R_0<1$ , then the disease-free equilibrium $E_{2_a}(K_a,0)$ of (14) is locally asymptotically stable, otherwise it is unstable. (iii) If $R_0>1$ , the endemic equilibrium $E_{3_a}(S_a^{*},I_a^{*})$ of (14) is locally asymptotically stable.

Proof. (i) The Jacobian matrix of system (14) evaluated at $E_{1_a}$ is given by

$$J(E_{1_a})=\left(\begin{array}{cc}{r}_a^{\alpha }& 0\\ 0& -({\mu }_a^{\alpha }+{\delta }_a^{\alpha })\end{array}\right)$$

By solving the characteristic equation $\mathit{det}(J(E_{1_a})-\lambda I_2)=0$ , with I₂ being the unit matrix, we get the following eigenvalues

$${\lambda }_1=r_a^{\alpha }>0,{\lambda }_2=-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })<0$$

Note that $\left|\mathit{arg}({\lambda }_1)\right|=0$ and $\left|\mathit{arg}({\lambda }_2)\right|=\pi$ . Since the first eigenvalue λ₁ does not satisfy $\left|\mathit{arg}({\lambda }_1)\right|>\frac{\alpha \pi }{2}$ for all $\alpha \in (0,1],$ therefore $E_{1_a}(0,0)$ is always unstable.

(ii) The Jacobian matrix $J(E_{2_a})$ is computed as

$$J(E_{2_a})=\left(\begin{array}{cc}-r_a^{\alpha }& -{\beta }_a^{\alpha }{K}_a\\ 0& {\beta }_a^{\alpha }{K}_a-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })\end{array}\right)$$

The corresponding eigenvalues are

$$\begin{array}{l}{\lambda }_1=-r_a^{\alpha }<0,{\lambda }_2={\beta }_a^{\alpha }{K}_a-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })\\ =({\mu }_a^{\alpha }+{\delta }_a^{\alpha })(R_0-1)\end{array}$$

Here, two cases arise depending on whether $R_0>1$ or $R_0<1$ .

Case 1. If $R_0<1$ , then we can see that $\left|\mathit{arg}({\lambda }_i)\right|=\pi >\frac{\alpha \pi }{2},\forall \alpha \in (0,1],i=1,2$ . Therefore, the equilibrium $E_{2_a}$ is locally asymptotically stable.

Case 2. If $R_0>1$ , then it is easy to see that $\left|\mathit{arg}({\lambda }_2)\right|=0$ . In this case, $E_{2_a}$ is unstable.

(iii) The Jacobian matrix $J(E_{3_a})$ is computed as

$$J(E_{3_a})=\left(\begin{array}{ll}\frac{-r_a^{\alpha }{S}_a^{*}}{K_a}& -{\beta }_a^{\alpha }{S}_a^{*}\\ {\beta }_a^{\alpha }{I}_a^{*}& 0\end{array}\right)$$

The corresponding characteristic equation becomes

$${\lambda }^2+\frac{r_a^{\alpha }{S}_a^{*}}{K_a}\lambda +{({\beta }_a^{\alpha })}^2{S}_a^{*}{I}_a^{*}=0$$

Since $S_a^{*}>0$ and $I_a^{*}>0$ if $R_0>1$ , all eigenvalues have negative real parts. Therefore, $\left|\mathit{arg}({\lambda }_i)\right|\in (\frac{\pi }{2},\pi ]>\frac{\alpha \pi }{2},\forall \alpha \in (0,1],i=1,2$ . Hence, the equilibrium $E_{3_a}$ is locally asymptotically stable. □

Lemma 3.2 If $R_0<1$ , then the disease-free equilibrium $E_{2_a}(K_a,0)$ of (14) is globally asymptotically stable.

Proof. We consider the following positive definite function $V_1(S_a,I_a)$ given by

$$V_1(S_a,I_a)=S_a-K_a-K_a\text{ln}\frac{S_a}{K_a}+I_a$$

Here, V₁ is a C¹ function such that $V_1>0$ for all values of $(S_a,I_a)\ne (K_a,0)$ and $V_1=0$ only at $(K_a,0)$ . Calculating the α order derivative of $V_1(S_a,I_a)$ along the solution of system (14), it follows from Lemma 2.3 that

$$\begin{array}{l}{\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_1(S_a,I_a)|}_{(3.3)}\le \frac{S_a-K_a}{S_a}\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{S}_a+\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{I}_a\\ =(S_a-K_a)\left(r_a^{\alpha }(1-\frac{S_a}{K_a})-{\beta }_a^{\alpha }{I}_a\right)\\ +{\beta }_a^{\alpha }{S}_a{I}_a-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a\\ =-\frac{r_a^{\alpha }}{K_a}(S_a-K_a{)}^2+\left({\beta }_a^{\alpha }{K}_a-{\mu }_a^{\alpha }-{\delta }_a^{\alpha }\right)I_a\\ =-\frac{r_a^{\alpha }}{K_a}(S_a-K_a{)}^2+({\mu }_a^{\alpha }+{\delta }_a^{\alpha })(R_0-1)I_a\end{array}$$

Note that if $R_0\le 1$ then $\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_1(S_a,I_a)\le 0$ for all $(S_a,I_a)\in {ℝ}_{+}^2$ , and $\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_1(S_a,I_a)=0$ at $E_{2_a}$ . Therefore, the only invariant set on which $\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_1(S_a,I_a)=0$ is the singleton $\left\{E_{2_a}\right\}$ . Then using Lemma 4.6 in Huo et al. ¹⁷ which generalizes the integer-order LaSalle Invariance Principle to fractional-order system, it follows that every nonnegative solution tends to $E_{2_a}$ when $R_0<1$ . Thus, $E_{2_a}$ is globally asymptotically stable if $R_0<1$ . □

Lemma 3.3 If $R_0>1$ , then the endemic equilibrium $E_{3_a}(S_a^{*},I_a^{*})$ of (14) is globally asymptotically stable.

Proof. Since $(S_a^{*},I_a^{*})$ is the equilibrium point of system (14), we have

$$\begin{array}{l}{r}_a^{\alpha }{S}_a^{*}(1-\frac{S_a^{*}}{K_a})-{\beta }_a^{\alpha }{S}_a^{*}{I}_a^{*}=0,\\ {\beta }_a^{\alpha }{S}_a^{*}{I}_a^{*}-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a^{*}=0\end{array}$$

(15)

Next, the global asymptotic stability of the equilibrium point $E_{3_a}(S_a^{*},I_a^{*})$ of system (14) will be demonstrated.

Let us define the Lyapunov function as

$$V_2(S_a,I_a)=S_a-S_a^{*}-S_a^{*}\text{ln}\frac{S_a}{S_a^{*}}+I_a-I_a^{*}-I_a^{*}\text{ln}\frac{I_a}{I_a^{*}}$$

Here, V₂ is a C¹ function such that $V_2>0$ for all values of $(S_a,I_a)\ne (S_a^{*},I_a^{*})$ and $V_2=0$ only at $(S_a^{*},I_a^{*})$ . Calculating the α order derivative of $V_2(S_a,I_a)$ along the solution of system (14), it follows from Lemma 2.3 that

$$\begin{array}{l}{\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_2(S_a,I_a)|}_{(3.3)}\le \frac{S_a-S_a^{*}}{S_a}\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{S}_a+\frac{I_a-I_a^{*}}{I_a}\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{I}_a\\ =(S_a-S_a^{*})\left(r_a^{\alpha }\right(1-\frac{S_a}{K_a})-{\beta }_a^{\alpha }{I}_a)\\ +(I_a-I_a^{*})({\beta }_a^{\alpha }{S}_a-({\mu }_a^{\alpha }+{\delta }_a^{\alpha }))\\ =(S_a-S_a^{*})\left[-\frac{r_a^{\alpha }}{K_a}(S_a-S_a^{*})-{\beta }_a^{\alpha }(I_a-I_a^{*})\right]\\ +{\beta }_a^{\alpha }(I_a-I_a^{*})(S_a-S_a^{*})\\ =-\frac{r_a^{\alpha }}{K_a}(S_a-S_a^{*}{)}^2\end{array}$$

Clearly, ${\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_2(S_a,I_a)|}_{\left(14\right)}\le 0$ for all $(S_a,I_a)\in {ℝ}_{+}^2$ , and ${\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_2(S_a,I_a)|}_{\left(14\right)}=0$ implies that $(S_a,I_a)=(S_a^{*},I_a^{*})$ . In fact, ${\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_2(S_a,I_a)|}_{\left(14\right)}=0$ implies that $S_a=S_a^{*}$ . When $S_a=S_a^{*}$ , from the first equation of system (14), we obtain

$$0=\frac{{\text{d}}^{\alpha }{S}_a^{*}}{\text{d}{t}^{\alpha }}=r_a^{\alpha }{S}_a^{*}(1-\frac{S_a^{*}}{K_a})-{\beta }_a^{\alpha }{S}_a^{*}{I}_a$$

(16)

it follows from the first equation of (15) and (16) that $I_a=I_a^{*}$ . That is, $\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_2(S_a,I_a)=0$ implies that $(S_a,I_a)=(S_a^{*},I_a^{*})$ . Therefore, the only invariant set on which $\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_2(S_a,I_a)=0$ is the singleton $\left\{E_{3_a}\right\}$ . By Lemma 4.6 in Huo et al., ¹⁷ it follows that the equilibrium point $E_{3_a}(S_a^{*},I_a^{*})$ of (14) is globally asymptotically stable if $R_0>1$ . This completes the proof of Lemma 3.3.

Analysis of the full system

Since the first four equations of system (2) are independent of the variable R_h, we only need to analyze the dynamical behavior of the following equivalent reduced system

$$\begin{array}{l}\frac{{\text{d}}^{\alpha }{S}_a}{\text{d}{t}^{\alpha }}=r_a^{\alpha }{S}_a(1-\frac{S_a}{K_a})-{\beta }_a^{\alpha }{S}_a{I}_a,\\ \frac{{\text{d}}^{\alpha }{I}_a}{\text{d}{t}^{\alpha }}={\beta }_a^{\alpha }{S}_a{I}_a-({\mu }_a^{\alpha }+{\delta }_a^{\alpha })I_a,\\ \frac{{\text{d}}^{\alpha }{S}_h}{\text{d}{t}^{\alpha }}={\Pi }_h^{\alpha }-{\beta }_h^{\alpha }{I}_a{S}_h-{\mu }_h^{\alpha }{S}_h,\\ \frac{{\text{d}}^{\alpha }{I}_h}{\text{d}{t}^{\alpha }}={\beta }_h^{\alpha }{I}_a{S}_h-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h\end{array}$$

(17)

System (17) always has two disease-free equilibria given by $E_{1_{ah}}\left(0,0,\frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }},0\right)$ and $E_{2_{ah}}\left(K_a,0,\frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }},0\right).$ If $R_0>1$ , then system (17) also has a unique endemic equilibrium given by $E_{3_{ah}}(S_a^{*},I_a^{*},S_h^{*},I_h^{*})$ .

Lemma 3.4 (i) The disease-free equilibrium $E_{1_{ah}}$ is always unstable. (ii) If $R_0<1$ , then the disease-free equilibrium $E_{2_{ah}}$ is locally asymptotically stable, otherwise it is unstable. (iii) If $R_0>1$ , then the disease-free equilibrium $E_{2_{ah}}$ is unstable and the endemic equilibrium $E_{3_{ah}}$ is locally asymptotically stable.

Proof. The characteristic equation of the Jacobian matrix at an arbitrary equilibrium $(S_a,I_a,S_h,I_h)$ takes the form

$$\begin{array}{l}(\lambda +{\beta }_h^{\alpha }{I}_a+{\mu }_h^{\alpha })(\lambda +{\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })\\ \times [\left(\lambda -(r_a^{\alpha }-2\frac{r_a^{\alpha }{S}_a}{K_a}-{\beta }_a^{\alpha }{I}_a)\right)\\ \times \left(\lambda -({\beta }_a^{\alpha }{S}_a-{\mu }_a^{\alpha }-{\delta }_a^{\alpha })\right)+({\beta }_a^{\alpha }{)}^2{S}_a{I}_a]=0\end{array}$$

(18)

• If $(S_a,I_a,S_h,I_h)=\left(0,0,\frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }},0\right)$ , the eigenvalues are ${\lambda }_1=-{\mu }_h^{\alpha },{\lambda }_2=-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha }),{\lambda }_3=r_a^{\alpha }>0$ and ${\lambda }_4=-({\mu }_a^{\alpha }+{\delta }_a^{\alpha }).$ Note that $\left|\mathit{arg}({\lambda }_1)\right|=\left|\mathit{arg}({\lambda }_2)\right|=\left|\mathit{arg}({\lambda }_4)\right|=\pi$ , and $\left|\mathit{arg}({\lambda }_3)\right|=0$ . Since the third eigenvalue λ₃ does not satisfy $\left|\mathit{arg}({\lambda }_3)\right|>\frac{\alpha \pi }{2}$ for all $\alpha \in (0,1)$ , therefore $E_{1_{ah}}$ is always unstable.

• If $(S_a,I_a,S_h,I_h)=\left(K_a,0,\frac{{\Pi }_h^{\alpha }}{{\mu }_h^{\alpha }},0\right)$ , we get from (18) that the eigenvalues are ${\lambda }_1=-{\mu }_h^{\alpha }<0,{\lambda }_2=-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })<0,{\lambda }_3=-r_a^{\alpha }<0$ , and ${\lambda }_4={\beta }_a^{\alpha }{K}_a-{\mu }_a^{\alpha }-{\delta }_a^{\alpha }=({\mu }_a^{\alpha }+{\delta }_a^{\alpha })(R_0-1)$ . If $R_0<1$ , then we get that $\left|\mathit{arg}({\lambda }_i)\right|=\pi >\frac{\alpha \pi }{2},\forall \alpha \in (0,1],i=1,2,3,4.$ Therefore, the equilibrium $E_{2_{ah}}$ is locally asymptotically stable. If $R_0>1$ , then it is easy to see that $\left|\mathit{arg}({\lambda }_4)\right|=0$ , hence $E_{2_a}$ is unstable.

• If $R_0>1$ and $(S_a,I_a,S_h,I_h)=E_{3_{ah}}(S_a^{*},I_a^{*},S_h^{*},I_h^{*})$ , by plugging (13) into (18), one gets

  $$\begin{array}{l}(\lambda +{\beta }_h^{\alpha }{I}_a^{\ast }+{\mu }_h^{\alpha })(\lambda +{\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })\\ \times \left[\left(\lambda +\frac{r_a^{\alpha }{S}_a^{\ast }}{K_a}\right)\lambda +{({\beta }_a^{\alpha })}^2{S}_a^{\ast }{I}_a^{\ast }\right]=0\end{array}$$

  (19)

It is easy to see that ${\lambda }_1=-({\beta }_h^{\alpha }{I}_a^{\ast }+{\mu }_h^{\alpha })<0,{\lambda }_2=-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })<0$ are two of the eigenvalues. The remaining two eigenvalues, λ₃ and λ₄, are the roots of the equation

$${(\lambda )}^2+\frac{r_a^{\alpha }{S}_a^{\ast }}{K_a}\lambda +{({\beta }_a^{\alpha })}^2{S}_a^{\ast }{I}_a^{\ast }=0$$

(20)

Since $S_a^{*}>0$ and $I_a^{*}>0$ if $R_0>1$ , all eigenvalues have negative real parts. Therefore, $\left|\mathit{arg}({\lambda }_i)\right|\in (\frac{\pi }{2},\pi ]>\frac{\alpha \pi }{2},\forall \alpha \in (0,1],i=1,2,3,4$ . Hence, the equilibrium $E_{3_{ah}}$ is locally asymptotically stable. □

Theorem 3.1 If $R_0<1$ , then the disease-free equilibrium $E_{2_{ah}}$ of the full reduced system (17) is globally asymptotically stable.

Proof. According to Lemma 3.2, the disease-free equilibrium $E_{2_a}(K_a,0)$ of (14) is globally asymptotically stable if $R_0\le 1$ . To prove the global stability of $E_{2_{ah}}$ , we only need to consider system (17) with the avian components already at the disease-free steady state, given by

$$\begin{array}{l}\frac{{\text{d}}^{\alpha }{S}_h}{\text{d}{t}^{\alpha }}={\prod }_h^{\alpha }-{\mu }_h^{\alpha }{S}_h,\\ \frac{{\text{d}}^{\alpha }{I}_h}{\text{d}{t}^{\alpha }}=-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h\end{array}$$

(21)

Clearly, we can obtain that $S_h\to S_h^{\ast },I_h\to 0$ if $t\to \infty$ . Hence, the disease-free equilibrium $E_{2_{ah}}$ is globally asymptotically stable. □

Theorem 3.2 If $R_0>1$ , then the endemic equilibrium $E_{3_{ah}}$ of the full reduced system (17) is globally asymptotically stable.

Proof. By Lemma 3.3, the endemic equilibrium $E_{3_a}(S_a^{*},I_a^{*})$ of (14) is globally asymptotically stable if $R_0>1$ . To prove the global stability of the equilibrium $E_{3_{ah}}$ , we only need to consider system (17) with the avian components already at the endemic steady state, given by

$$\begin{array}{l}\frac{{\text{d}}^{\alpha }{S}_h}{\text{d}{t}^{\alpha }}={\prod }_h^{\alpha }-{\beta }_h^{\alpha }{I}_a^{*}{S}_h-{\mu }_h^{\alpha }{S}_h,\\ \frac{{\text{d}}^{\alpha }{I}_h}{\text{d}{t}^{\alpha }}={\beta }_h^{\alpha }{I}_a^{*}{S}_h-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h\end{array}$$

(22)

We can easily deduce that subsystem (22) has a unique positive equilibrium $(S_h^{*},I_h^{*})$ which is locally asymptotically stable.

To prove the global stability of the positive equilibrium $(S_h^{*},I_h^{*})$ of subsystem (22), we choose a Lyapunov function as follows

$$V_3(S_h,I_h)=S_h-S_h^{*}-S_h^{*}\text{ln}\frac{S_h}{S_h^{*}}+I_h-I_h^{*}-I_h^{*}\text{ln}\frac{I_h}{I_h^{*}}$$

Here, V₃ is a C¹ function such that $V_3>0$ for all values of $(S_h,I_h)\ne (S_h^{*},I_h^{*})$ and $V_3=0$ only at $(K_a,0)$ . Calculating the α order derivative of $V_3(S_h,I_h)$ along the solution of system (22), it follows from Lemma 2.3 that

$$\begin{array}{l}{\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_3(S_h,I_h)|}_{(3.11)}\\ \le \frac{S_h-S_h^{*}}{S_h}\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{S}_h+\frac{I_h-I_h^{*}}{I_h}\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{I}_h\\ =\left(1-\frac{S_h^{*}}{S_h}\right)({\Pi }_h^{\alpha }-{\beta }_h^{\alpha }{I}_a^{*}{S}_h-{\mu }_h^{\alpha }{S}_h)\\ +\left(1-\frac{I_h^{*}}{I_h}\right)({\beta }_h^{\alpha }{I}_a^{*}{S}_h-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h)\end{array}$$

(23)

Since $(S_h^{*},I_h^{*})$ is the equilibrium point of system (22), we have

$$\begin{array}{l}{\prod }_h^{\alpha }-{\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}-{\mu }_h^{\alpha }{S}_h^{*}=0,\\ {\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h^{*}=0\end{array}$$

which leads to

$${\prod }_h^{\alpha }={\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}+{\mu }_h^{\alpha }{S}_h^{*},{\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha }=\frac{{\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}}{I_h^{*}}$$

It follows that

$$\begin{array}{l}\left(1-\frac{S_h^{*}}{S_h}\right)({\Pi }_h^{\alpha }-{\beta }_h^{\alpha }{I}_a^{*}{S}_h-{\mu }_h^{\alpha }{S}_h)\\ =\left(1-\frac{S_h^{*}}{S_h}\right)({\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}+{\mu }_h^{\alpha }{S}_h^{*}-{\beta }_h^{\alpha }{I}_a^{*}{S}_h-{\mu }_h^{\alpha }{S}_h)\\ ={\mu }_h^{\alpha }{S}_h^{*}\left(2-\frac{S_h}{S_h^{*}}-\frac{S_h^{*}}{S_h}\right)+2{\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}\\ -{\beta }_h^{\alpha }{I}_a^{*}{S}_h-{\beta }_h^{\alpha }{I}_a^{*}\frac{{(S_h^{*})}^2}{S_h}\end{array}$$

(24)

and

$$\begin{array}{l}\left(1-\frac{I_h^{*}}{I_h}\right)({\beta }_h^{\alpha }{I}_a^{*}{S}_h-({\mu }_h^{\alpha }+{\delta }_h^{\alpha }+{\gamma }^{\alpha })I_h)\\ =\left(1-\frac{I_h^{*}}{I_h}\right)\left({\beta }_h^{\alpha }{I}_a^{*}{S}_h-\frac{{\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}}{I_h^{*}}{I}_h\right)\\ ={\beta }_h^{\alpha }{I}_a^{*}{S}_h+{\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}-{\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}\frac{I_h}{I_h^{*}}-{\beta }_h^{\alpha }{I}_a^{*}{S}_h\frac{I_h^{*}}{I_h}\end{array}$$

(25)

Plugging (24) and (25) into (23), we obtain

$$\begin{array}{l}{\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_3(S_h,I_h)|}_{(3.11)}\le {\mu }_h^{\alpha }{S}_h^{*}\left(2-\frac{S_h}{S_h^{*}}-\frac{S_h^{*}}{S_h}\right)\\ +{\beta }_h^{\alpha }{I}_a^{*}{S}_h^{*}\left(3-\frac{S_h^{*}}{S_h}-\frac{I_h}{I_h^{*}}-\frac{S_h}{S_h^{*}}\frac{I_h^{*}}{I_h}\right)\end{array}$$

(26)

Since the arithmetic mean exceeds the geometric mean, we have

$$\begin{array}{l}2-\frac{S_h}{S_h^{*}}-\frac{S_h^{*}}{S_h}\le 0,\\ 3-\frac{S_h^{*}}{S_h}-\frac{I_h}{I_h^{*}}-\frac{S_h}{S_h^{*}}\frac{I_h^{*}}{I_h}\le 0\end{array}$$

Hence, ${\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_3(S_h,I_h)|}_{\left(22\right)}\le 0$ for all $(S_h,I_h)\in {ℝ}_{+}^2$ , and ${\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }}{V}_3(S_h,I_h)|}_{\left(22\right)}=0$ implies that $(S_h,I_h)=(S_h^{*},I_h^{*})$ . Again, by Lemma 4.6 in Huo et al., ¹⁷ it follows that $S_h\to S_h^{*}$ and $$ as $t\to \infty$ . Therefore, the endemic equilibrium $E_{3_{ah}}$ of the full system (17) is globally asymptotically stable. □

Now we can state our results for the full avian-human model (2).

Theorem 3.3 (i) The disease-free equilibrium E₁ of model (2) is always unstable; (ii) If $R_0<1$ , then the disease-free equilibrium E₂ of model (2) is globally asymptotically stable. (iii) If $R_0>1$ , then then the disease-free equilibrium E₂ of model (2) is unstable but the endemic equilibrium E₃ of model (2) is globally asymptotically stable.

Remark 3.1 The above stability results we obtained for the fractional order model is also true for α = 1, which has been obtained by Liu et al. ²⁴

Numerical simulations

In this section, numerical simulations are provided to substantiate the theoretical results established in the previous section of this paper. In addition, we also show the effects of fractional order α on the stability of the equilibrium points, these numerical simulations are very important from biological point of view.

The existing numerical methods for fractional differential equations include finite difference, Galerkin finite element, spectral method and so on. The time stepping scheme employed in our simulations is the so-called $2-\alpha$ order finite difference method, proposed by Lin and Xu ²⁸ for the temporal discretization of the time fractional diffusion equation. It was proved by Lin and Xu that this scheme is unconditionally stable, and the convergence rate is of $2-\alpha$ order, where $\alpha$ is the order of the fractional derivative. In what follows, all the parameter values are chosen hypothetically due to the unavailability of real world data.

We consider the parameter values as $r_a=0.02,K_a=50,{\beta }_a=0.002,{\mu }_a=0.1,{\delta }_a=0.2,{\Pi }_h=60,{\beta }_h=0.006,{\mu }_h=0.003,{\delta }_h=0.08,\gamma =0.1$ and initial point $S_a(0)=80,I_a(0)=10,S_h(0)=60,I_h(0)=20,R_h(0)=10$ .

In Figure 1, we plot the variation of the value of the reproduction number for different values of α. It is observed that as α is decreased from 1, R₀ goes to values more than 1. This indicates that at an order $\alpha \in (0.7,0.8)$ there is a bifurcation from a disease-free equilibrium ( $R_0<1$ ) to a stable endemic equilibrium ( $R_0>1$ ) in the model. Thus, analytically, α is a bifurcation parameter of the model.

Figure 1.

Variation of the reproduction number, R₀, versus the order of the fractional derivative α.

Let $\alpha =0.8,0.9,1.$ The corresponding reproduction numbers are $0.7978,0.5160,0.3333$ , respectively. Following Theorem 3.1, the disease-free equilibrium $E_{2_{ah}}$ of (2) is stable. In Figure 2, we plot the solutions of the reduced system (17) with different values of α. It shows that all populations remain stable for all values of α though solutions reach to equilibrium value more slowly for larger value of α.

Figure 2.

Dynamics of model (17) for $\alpha =0.8,0.9,1$ .

Now we choose $\alpha =0.4,0.5,0.6$ . The corresponding reproduction numbers are 4.5080, 2.9289, 1.9008, respectively. Following Theorem 3.2, the endemic equilibrium $E_{3_{ah}}$ of (2) is stable. In Figures 3 to 5, we report the dynamics of model (17) for $\alpha =0.4,0.5,0.6$ . As expected, the disease is endemic, since all the solutions of model (17) approach to the endemic equilibrium $E_{3_{ah}}$ .

Figure 3.

Global asymptotic stability of the endemic equilibrium point for $\alpha =0.4$ . Upper: Time series of $S_a,I_a,S_h,I_h$ ; Lower: phase portrait of the avian system (left) and the reduced human system (right) with different initial value.

Figure 4.

Global asymptotic stability of the endemic equilibrium point for $\alpha =0.5$ . Upper: Time series of $S_a,I_a,S_h,I_h$ ; Lower: phase portrait of the avian system (left) and the reduced human system (right) with different initial value.

Figure 5.

Global asymptotic stability of the endemic equilibrium point for $\alpha =0.6$ . Upper: Time series of $S_a,I_a,S_h,I_h$ ; Lower: phase portrait of the avian system (left) and the reduced human system (right) with different initial value.

In Figures 6 and 7, we plot the variation of the avian infective individuals and the human infective individuals for different values ν. One can see from Figures 6 and 7 that both I_a and I_h tends to a stationary state over a longer period of time with decreasing the value of α.

Figure 6.

Numbers of avian infective individuals for different values α.

Figure 7.

Numbers of human infective individuals for different values α.

Concluding remarks

In this paper, we have proposed a fractional order avian-human influenza epidemic model with logistic growth for avian population. We prove that the solution of the model exists uniquely and all solutions remain positive and bounded whenever they start with positive initial value, thus justifying the well-posedness of a biological model. We show that our system contains three equilibrium points. The trivial equilibrium point is always unstable, implying that all populations can not go to extinction simultaneously. The disease free equilibrium is locally and globally asymptotically stable if $R_0<1$ , and the endemic equilibrium is locally asymptotically and globally stable if $R_0>1.$ Our numerical results show that solutions of fractional order system converge to the respective equilibrium more slowly as the order of the differential equation becomes smaller.