MODELLING THE DYNAMICS OF AVIAN INFLUENZA WITH NONLINEAR RECOVERY RATE AND PSYCHOLOGICAL EFFECT∗

In this paper, a SI-SEIR type avian influenza epidemic model with psychological effect, nonlinear recovery rate and saturation inhibition effect is formulated to study the transmission and control of avian influenza virus. By setting the basic reproductive number as the threshold parameter and constructing Lyapunov function, Dulac function and using the Li-Muldowney’s geometry approach, we prove the local and global stability of disease-free equilibria and endemic equilibrium. Theoretical analysis are carried out to show the role of the saturation inhibition effect, psychological effect and effective medical resources in this model, and numerical simulations are also given to verify the results.


Introduction
Avian influenza refers to the disease which caused by infection with avian (bird) influenza (flu) Type A viruses. These viruses occur naturally among wild aquatic birds worldwide and can infect domestic poultry, other birds and animal species. Avian flu viruses do not normally infect humans. However, sporadic human infections with avian flu viruses have occurred [7].
Avian influenza A viruses are designated as highly pathogenic avian influenza (HPAI) or low pathogenic avian influenza (LPAI). The other three subtypes AH9, AH5 and AH7 can simultaneously infect humans and birds. Viruses (H1-H9) found in poultry and wild birds worldwide belong to the same category (LPAI). Rarely, sporadic cases of human infection with H9N2 leading to mild upper respiratory disease have been reported. Reports from more than 15 countries indicate that the H5N1 viruses can infect humans, most of which directly cause severe pneumonia, with a mortality rate close to 60%. The human infected H7 viruses are not common, but it has also been documented that direct contact with infected birds can infect H7 viruses, especially during outbreaks of H7 viruses in poultry. LPAI (H7N2, H7N3 and H7N7) viruses would result in mild-to-moderate diseases, while HPAI(H7N3,H7N7) virus infection can result in mild-to-moderate or even lethal diseases. Human cases of H7N9 virus infection were first detected in 2013, including death and severe respiratory illnesses [8].
In recent years, avian influenza has been rampant, and its mutation rate is very fast, which has a non-negligible impact on human life, health as well as social economy. In view of this, according to the epidemiological characteristics of avian influenza, researchers have established a large number of mathematical models to study, and combined with its dynamical behavior, remarkable achievements have been made in the prevention and control of avian influenza [1,19,23]. Iwami et al. [18] proposed an Ordinary Differential Equation (ODE) model in 2007 to describe the dynamics of avian influenza between human and avian populations. Then, many different mathematical models were designed for the H5N1 avian influenza virus [2,12,30,33]. A series of analysis has also been carried out on the epidemiology of influenza A H7N9 influenza in recent years [5,21,27,36]. However, most of them have not considered the latent state of avian influenza virus in human body, and it does exist according to reported cases, so it is very necessary to take it into account in models.
Human behavior and social reaction have great significance to the spread of infectious disease [4,[13][14][15], since it is a key factor in disease control efforts. Wang et al. [34] found that 70% urban respondents who participated in the survey said that since the first case of H7N9 were detected in China in March 2013, their visit to live poultry market have been relatively reduced. In addition, the research of Wu et al. [35] also suggested that people experience changes in protective behaviors such as reducing access to live poultry markets and buying live poultry in the context of continuous H7N9 outbreaks. Although human behavior and social responses during the transmission of infectious diseases are often reported, few systematic studies have been conducted on their effects. It must be admitted that it is difficult to combine social behavioral responses with human behavior and social responses in mathematical models.
As we know, incidence rate of a disease is one of crucial factors in the transmission of disease. It is worth noting that in almost all models related to avian influenza, the incidence between susceptible birds and infected birds (susceptible humans and infected birds) has taken in form of bilinear interaction, which is constantly increasing unbounded, whereas as survey [37] revealed that if the case was reported by the media, people would go to live poultry market less, incidence of human will decline. In order to better explain the practical significance of infectious disease model, the non-monotonic incidence rate is used here which was proposed by Liu et al. [22] perfectly describes this phenomenon caused by psychological effects(also see [20,21,38]); Similarly, with the increasing number of infected birds, poultry farmers will be highly vigilant and take the corresponding protective measures to make the poultry incidence rate reach saturation, which can be described by the saturation incidence function used in [6,31].
There is one thing in common in the classical avian influenza model, that is, the recovery rate is always assumed to be constant, which is equivalent to the default public medical resources are always sufficient, and it is obviously unreasonable by the following two points. First of all, hospital resources (such as doctors, medicines, beds) are limited for the public; Secondly, according to the cases reported by the CDC [7], there are some similarities between human infection with avian influenza virus and common influenza virus in terms of infection time and early clinical manifestations, then some of the available hospital resources have already occupied. Proportion between the number of beds and the population, ie., the number of hospital beds available per 10,000 people, is widely used by health planners to estimate public resource availability [17]. Abdelrazec et al. [3] established a dengue propagation dynamics model which considered a recovery rate function restricted by hospital bed-population proportion and the number of infected. They proved that the model has oscillations and backward bifurcation due to the limited resources. At this time, controlling the basic regeneration number R 0 < 1 is not enough to ensure the eradication of the disease. Therefore, the recovery rate function which is defined by [32] will be introduced in this paper, further learning the impact of the limited availability of existing medical resources on the spread of avian influenza virus.
The rest of this paper is outlined as follows: Section 2 constructs an avian influenza epidemic model based on the above discussion, and gives a general explanation of parameters; In Section 3, the avian-only sub-model and its detailed mathematical analysis are given; Section 4 presents the whole mathematical analysis of the model in details; Section 5 carries out the numerical simulation to the full avian influenza infectious disease model; Finally, conclusions and discussions are also given in Section 6.

Model discription
In this section, we assume that the avian influenza virus does not spread from person to person and mutate. The avian population is classified into two subclasses: susceptible S a (t) and infected I a (t), respectively. The human population is classified into four subclasses: susceptible S h (t), exposed E h (t), infected I h (t) and recovered human R h (t), respectively. Then we have the model (2.1) as follows: with the following assumptions: (1) The net growth rate function of susceptible poultry is subject to Logistic growth g(S a ) = m a S a (1 − S a K a ) : R + → R is continuous, where m a and K a are the intrinsic growth rate and the maximum environmental capacity of poultry, respectively. R = (−∞, ∞) and R + = (0, ∞); (2) All new population (poultry) are subordinate to susceptible populations (poultry), and new population recruitment rate is represented by Λ h ; (3) The avian influenza virus is transmitted only between infected poultry and susceptible populations; (4) Infected poultry can't recover and maintain the disease state until death. Infected population can recover and gain permanent immunity; (5) The incidence rate between susceptible poultry (or population) and infected poultry is where µ 1 is the highest per capita recovery rate, µ 0 is the minimum per capita recovery rate due to lack of clinical resources, b is the ratio between the number of hospital bed and the population.
Thus, we obtain the above model, where β a is the transmission rate from infective avian to susceptible avian, µ a is the natural death rate of the avian population, δ a is the disease-related death rate of the infected avian,β 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 death rate of the infected human. The parameters α and c are constants measuring the inhibition effect and psychological effect, respectively.
We consider system (2.1) with initial conditions The parameters of system (2.1) are shown in Table 1. Human recruitment and birth rates 30 [23] βh The rate of disease transmission between infected poultry and susceptible human 5 * 10 −9 ∽ µh The natural mortality of population Human disease-related mortality 0.077 The covert rate of disease from latent to infected state The lowest rate of recovery in humans (0.067 ∼ 0.100) The highest recovery rate in humans (µ0 ∼ 10) Hospital beds-population ratio (0, 20) Proof. Setting N a = S a + I a and adding the first two equations of system (2.1), we have Then, we obtain then it follows Therefore, D is positively invariant. Next, we will calculate the basic reproduction number of system (2.1) .
Using the method of the next generation matrix in [10,11], system (2.1) can be rewritten as Then, we have Then the basic reproduction number of system (2.1) is

Analysis of Avian-only sub-model
We first discuss the dynamic behavior of the following avian sub-model: (3.1)

Existence of equilibria in system (3.1)
Obviously, A a (0, 0) and B a (K a , 0) are two disease-free equilibria of system (3.1). Now we consider the existence of endemic equilibrium C a (S * a , I * a ). Then S * a and I * Substituting S * a into (3.2), we obtain 3) has two distinct real roots. According to the Vieta theorem, we have .
(3) If R 0 < 1, then A > 0, B > 0 and C > 0. It follows from Descartes' rule of signs that equation (3.3) does not have a positive root. From the above discussions, we know that system (3.1) exists a unique endemic equilibrium C a (S * a , I * a ) when R 0 > 1.

Local stability of equilibria in system (3.1)
Then the corresponding characteristic equation is 0) is a saddle node bifurcation point which is locally asymptotically stable.
Actually, let S a = S 1 + K a and I a = I 1 . System (3.1) is transformed into Using S 1 = x − kβy m a , I 1 = y, t = − τ m a , the above system can be rewritten as and expanding its Taylor series, we get According to [29, p147-p152], we can get that B a (K a , 0) is a saddle-node point when R 0 = 1. Further, it is locally asymptotically stable.
(1 + αI * a ) 3 . The first equation of (3.2) can be deduced to Then we conclude that all eigenvalues have strictly negative real parts. Therefore, the endemic equilibrium C a is locally asymptotically stable when R 0 > 1.

Global stability of equilibria in system (3.1)
We can also get {(S a , I a ) ∈ intD 1 )| dV dt | (3) = 0}={(S a , I a )|S a = K a , I a = 0}={B a }. According to LaSalle's invariant principle, the disease-free equilibrium B a = (K a , 0) is globally asymptotically stable if R 0 ≤ 1. S a I a . By simple calculation, we obtain By using Bendixon Dulac criterion, we know that system (3.1) can't have a closed orbit in D 1 , namely, the endemic equilibrium C a (S * a , I * a ) is globally asymptotically stable.

Analysis of full influenza epidemic model
Since the first five equations of system (2.1) are independent of R h , we only need to consider the following system:

Substituting (4.3) into the third equation of (4.2)
, we obtain It is clear that f (0) = m 0 > 0 and f (+∞) < 0. Applying the interval-value theorem of continuous functions, (4.4) has at least one positive root. Suppose that (4.4) has four real roots I h1 , I h2 , I h3 , I h4 . By the Vieta theorem, we have From (4.5) and (4.6), we can deduce that there are one negative root and three positive roots (or one positive root and three negative roots). Assume that (4.4) has a negative root I h1 and three positive roots I h2 , I h3 , I h4 . It follows from (4.6) that I h1 + I h2 < 0, I h1 + I h3 < 0, I h1 + I h4 < 0.
Proof. The Jacobian matrix of system (4.1) is (1) When (S a , I a , S h , E h , I h ) = (0, 0, Λ h µ h , 0, 0), we know that the characteristic equation corresponding to J(A ah ) always has a positive root λ = m a . Hence, A ah is always unstable.
By some calculations, we can get a 2 a 1 − a 3 > 0. Then by using the Routh-Hurwitz criterion, we obtain that the endemic equilibrium C ah is locally asymptotically stable.

Global stability of equilibria in system (4.1)
We first use the method in [9] to prove the global stability of B ah (K a , 0, Λ h µ h , 0, 0). Rewrite system (4.1) as Z).
The disease-free equilibrium of system (4.9) is Q 0 = (X 0 , 0). It is globally asymptotically stable if and only if the following conditions are satisfied.
where Γ is the constant attraction region of system (4.9).
. System (4.1) can be rewritten as Clearly, G(X, 0) = 0 and the disease-free equilibrium of system (4.1) is Thus, X 0 = (K a , Λ h µ h ) is globally asymptotically stable. Then condition (H 1 ) is satisfied. Let Thus, B ah of system (4.1) is globally asymptotically stable when R 0 < 1.
According to the above discussions, we can reduce system (4.1) to (4.10) We shall use Li-Muldowney's geometry method [24] to study the global stability of C ′ ah (S * h , E * h , I * h ) in system (4.10) . Let | · | denote a vector norm in R n and also denote the induced matrix norm in R n×n , the space of all n × n matrices. For matrix A in R n×n , the Lozinskǐmeasure or the logarithmic norm of A with respect to | · | (see [25]) is Let y(t) be a solution of linear differential equatioṅ where A(t) is m × m matrix-valued continuous function. Then, we have Let B be an n × n matrix. The second additive compound matrix of B, denoted by B [2] , is an Consider the following autonomous systeṁ where f : Ω → R n , Ω ⊂ R n is an open set and simply connected and f ∈ C 1 (Ω). Let x(t, x 0 ) be the solution of system (4.11) such that x(0, x 0 ) = x 0 . Let x * be an equilibrium of system (4.11), i.e., f (x * ) = 0. A set K is said to be absorbing in Ω for system (4.11) if x(t, K 1 ) ⊂ K for each compact set K 1 ⊂ Ω and sufficiently large t. Assume the following assumptions hold: (H 3 ) System (4.11) has a unique equilibrium point x * in Ω.
. Definē and the matrix Q f is obtained by replacing each entry q ij of Q by its derivative in the direction of f , (q ij )f , and J [2] is the second additive compound matrix of the Jacobian matrix J of system (4.11) . The following lemma in Li and Muldowney [24] will be used here.

Lemma 4.1 ( [24]). Assume that Ω is simply connected and assumptions (H 3 ) and
(H 4 ) hold. Then, the unique equilibrium x * of system (4.11) is globally asymptotically stable in Ω if there exist a function Q and a Lozinskǐ measure µ such that q 2 < 0.
From the above statement, we now state our main result. Proof. The Jacobian matrix of system (4.10) is where By Theorems 4 and 5, we have (1) If R 0 > 1, then system (4.1) has a unique endemic equilibrium which is locally asymptotically stable in D 2 . Assumption (H 3 ) holds.
(2) If R 0 > 1, then A ah and B ah are unstable. The instability of A ah ∈ ∂D 2 and B ah ∈ ∂D 2 imply the uniform persistence [16], i.e., there exists a constant m > 0 such that Due to D 2 is bounded, there must exist a compact set in the interior of D 2 which is absorbing for system (4.10). So, (H 3 ) is verified. From system (4.10), we get and

Take any vector
Let µ(·) be the Lozinskǐ measure with this vector norm, which can be obtained according to [28], then where |B 12 |,|B 21 | are the matrix norms with respect to l 1 vector form Then we have Hence, we have From the second and third equations of system (4.10), we obtaiṅ Then Since holds, then Along each solution (S h , E h , I h ) of system (4.10), we have This means thatq According to the analysis above-mentioned, if R 0 > 1 and (4.13) hold, C ′ ah (S * h , E * h , I * h ) of system (4.10) is globally asymptotically stable. Based on the above-mentioned analysis, then we can obtain the following theorem.  1) is globally asymptotically stable when condition (4.13) is not satisfied (see Figure 1(b)).

Discussions and conclusions
In order to investigate the influence of saturation inhibition, psychological effect and medical resources on the spread of avian influenza, we proposed an SI-SEIR avian influenza epidemic model with nonlinear incidence rate and nonlinear recovery rate functions and presented the detailed analysis. Firstly, the result of the avian sub-model (3.1) is given. It shows that the disease-free equilibrium A a is always unstable. Moreover, when the basic reproduction number R 0 < 1 (R 0 > 1), the disease-free equilibrium B a (the endemic equilibrium C a ) is globally asymptotically stable in D 1 ; Secondly, the dynamic behavior of the whole avian influenza model is analyzed. From the theoretical analysis, the disease-free equilibrium A ah is unstable; when the basic reproduction number R 0 < 1, the disease-free equilibrium B ah is globally asymptotically stable in D 2 (see Figure 1(a)), i.e., the disease will be eradicted; when R 0 > 1 and (4.13) holds, the endemic equilibrium C is globally asymptotically stable in D 2 (see Figure 1(b)), that is to say, the avia inluenza will persist for a long time; Finally, numerical simulations (Figures 1-5) are also given and some conclusions can be drawn as follows.
(1) The saturation effect inhibition coefficient α, the psychological effect coefficient c, and the bed-to-population ratio b do not change the stability of system (2.1).
(2) When the inhibitory effect α, the ratio of hospital beds-population b and the psychological effect c increase together or seperately, the peak value of infected human will drop, and then the final size of the number of infected individual is relatively reduced (see .
In order to better embody the epidemiology of avian flu, it will be very interesting to induce the periodic incidence rate and the periodic recovery rate in avian influenza epidemic model, we leave it in the future work.