Global stability analysis of an SVEIR epidemic model with general incidence rate

In this paper, a susceptible-vaccinated-exposed-infectious-recovered (SVEIR) epidemic model for an infectious disease that spreads in the host population through horizontal transmission is investigated, assuming that the horizontal transmission is governed by an unspecified function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f(S,I)$\end{document}f(S,I). The role that temporary immunity (vaccinated-induced) and treatment of infected people play in the spread of disease, is incorporated in the model. The basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R}_{0}$\end{document}R0 is found, under certain conditions on the incidence rate and treatment function. It is shown that the model exhibits two equilibria, namely, the disease-free equilibrium and the endemic equilibrium. By constructing a suitable Lyapunov function, it is observed that the global asymptotic stability of the disease-free equilibrium depends on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R}_{0}$\end{document}R0 as well as on the treatment rate. If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R}_{0}>1$\end{document}R0>1, then the endemic equilibrium is globally asymptotically stable with the help of the Li and Muldowney geometric approach applied to four dimensional systems. Numerical simulations are also presented to illustrate our main results.


Introduction
Mathematical modeling enjoys popularity in both preventing and controlling infectious diseases such as severe acute respiratory syndrome (SARS) [1], human immunodeficiency virus infection/acquired immune deficiency syndrome (HIV/AIDS) [2], H5N1 (avian flu) [3] and H1N1 (swine flu) [4]. In recent years, a lot of efforts have been made to develop realistic diseases and further study the asymptotic behavior of such epidemic models [5]. In the field of studying epidemic model behavior, one of the most important parts is to analyze steady states together with their stability [6]. In general, there are two distinct techniques named Lyapunov's direct method and Li-Muldowney's geometric approach to give sufficient conditions of global stability for the equilibrium states (see, for example, [7][8][9][10][11][12][13][14]). We would like to mention some related work concerned with the existence of positive solutions for the discrete fractional boundary value problem [15], the sensitivity analysis for optimal control problems governed by nonlinear evolution inclusions [16] and the nonexistence of global in time solution of the mixed problem for the nonlinear evolution equation with memory generalizing the Voigt-Kelvin rheological model [17].
It is well known that the rate of incidence plays the main part in modeling infectious diseases. The rise and fall of epidemics can be influenced by some factors, such as density of population and life style [18,19]. Many researchers have adopted different nonlinear incidence rates in their works. For more details, we refer the reader to [8][9][10][11][12][13][14][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] and the references therein. When it comes to control of a disease, it is generally known that the spread of many diseases can be prevented by vaccinating. When massive vaccination is impossible, the second stage of defensive mechanism could be medical treatment. Individuals need to bear in mind that the treatment is an indispensable way to take precautions for some diseases (for instance, measles, phthisis and influenza). In recent years, many treatment functions have been introduced by several authors to study some epidemic models under different conditions (see, for instance, [12,14,27,31,[35][36][37][38]).
Recently, Dénes and Röst [27] investigated the following SI model: where a population of constant size (assumed to be equal to 1) is divided into three compartments: susceptible (denoted by S), infected (denoted by I) and recovered (denoted by R). The transmission of the infection is governed by the incidence rate f (S, I) and μ is birth rate as well as the death rate of the susceptible class. The nonlinear function g(I) denotes the sum of the death rate and the recovery rate for the infected individuals satisfying g(0) = 0 and g(I) > 0 for I > 0. Using a Dulac function approach, which aims at eliminating the existence of the periodic solution and proving the global stability by the Poincaré-Bendixson theorem (see [39], p. 54), they obtained the global stability of the disease-free equilibrium and the endemic equilibrium for system (1.1). Very recently, Upadhyay et al. [12] considered the following e-epidemic model: All the parameters in model (1.2) are positive and are defined as follows: S, E, I, R and V represent the number of susceptible nodes, exposed nodes, infectious nodes, recovered nodes and vaccinated nodes at time t, respectively; A is the recruitment rate of new nodes, c is the half saturation constant for susceptible nodes S, α is the contact rate or the rate of transfer of virus from an infectious node to the susceptible node, η is rate at which the vaccinated nodes lose their immunity and join the susceptible class, β is the maximal treatment capacity of a network, δ 0 is the natural crashing rate of nodes all classes, a is the half saturation constant for an infected node I, μ is the vaccination rate coefficient, δ 3 is the virus induced crashing rate and δ 1 , δ 2 are the state transition rates.
Using a Lyapunov function and a geometric approach, they obtained the global stability of virus-free equilibrium and endemic equilibrium for system (1.2). As pointed out by Liu and Yang [11], due to the high similarity between computer virus and biological virus, it is acceptable to establish dynamical models describing biological virus among a population by modifying an e-epidemic model. Thus, it is interesting and important to extend model (1.2) to study the biological virus in the infectious disease.
Inspired by these research results above, in this paper, we consider the following system with five compartments: where S(t), E(t), I(t), R(t), V (t) are the number of susceptible population, exposed population, infective population, recovered population, vaccinated population, respectively. The two-variable function f (S, I) represents incidence rate and the nonlinear function g(I) denotes the removal rate of infective individuals because of the treatment of infective. The initial conditions for system (1.3) are as follows: denotes the total number of high-risk human population at time t. The model parameters of system (1.3) are described as follows: A: the rate at which new individuals (including newborns and immigrants) enter the susceptible population, δ 0 : natural death rate of population all classes, η: the rate at which the vaccinated population lose their immunity and join the susceptible class, μ: vaccination rate coefficient, δ 1 : the rate at which exposed population become infective, δ 2 : natural recovery rate of infective population, δ 3 : disease-related death rate of infective population. Model (1.3) involves certain assumptions which consist of the following: (i) The new individuals enter the population with a constant rate and all the new individuals are susceptible. (ii) Susceptible individuals move to exposed class by adequate contact with infective individuals and after some time (i.e., latency period), they become infectious and move to infectious class. (iii) The infectious individuals are assumed to leave the infectious class as a result of natural death and disease-related death as well as recovery of infected individuals. (iv) After recovery the individuals become immunized and hence they are no longer susceptible to it.
(v) It is assumed that a fraction of susceptible individuals get vaccinated and join the vaccinated class. A part of vaccinated individuals may lose their immunity and rejoin the susceptible class. It is easy to see that system (1.3) includes (1.1) and (1.2) as special cases and so model (1.3) provides a uniform setting for the computer virus and biological virus studies. Following the classical assumptions [27,40], it is reasonable to suppose that the transmission of the infection is governed by an incidence rate f (S, I) in model (1.3). Moreover, as pointed out by Wang [31], the recovery rate is naturally dependent on the number of infected individuals provided the health care resources are constrained and so it is natural to use the nonlinear function g(I) as the treatment function in model (1.3).
The main purpose of this paper is to derive the expression for the basic reproduction number and further show the global stability of disease-free as well as endemic equilibria by the aid of Lyapunov function and Li-Muldowney geometric approach applied to four dimensional systems. This paper is organized as follows. In Sect. 2, some elementary assumptions on the functions f and g will be given, and the basic reproduction number R 0 is provided. Also the equilibrium points are discussed. The global stability of disease-free equilibrium and endemic equilibrium are analyzed in Sects. 3 and 4, respectively. All our important analytical findings are numerically verified with the help of Mathlab in Sect. 5. Finally, a brief conclusion is given in Sect. 6.

Basic reproduction number and equilibrium
To define the basic reproduction number R 0 and indicate the existence of equilibrium, we give some hypotheses.
(H1) f : (1) It is easy to check that the classes of f (S, I) satisfying (H1) include incidence rates such as for β, a, b > 0 and 0 ≤ q ≤ 1.
(2) It is straightforward to show that the classes of g(I) satisfying (H2) include removal rates such as (3) By hypothesis (H2), we know that Φ(I) = g(I) I is a monotone decreasing function on I > 0. (4) The assumption (H3) is equivalent to the following inequality: which can be found in [27]. (5) By the assumptions, it is easy to find that system (1.3) always has a disease-free equilibrium We shall assume that (H1), (H2) and (H3) hold in the rest of this paper. Now we define the basic reproduction number R 0 for model (1.3) as follows: In order to find the positive equilibria of model (1.3), set and V = μS m 4 . Substituting the above equalities into the second equation in (2.1), one has Then it is easy to see that the positive equilibrium points of system (2.1) are given by zeros of F in the interval (0, Aδ 1 m 3 m 2 ]. We denote G(I) = m 3 I + g(I) -Aδ 1 m 2 for convenience. Then it is easy to see Therefore, we conclude that G(I) has at least one root namedĨ in the interval (0, Aδ 1 m 3 m 2 ]. That is, m 3Ĩ + g(Ĩ) = Aδ 1 m 2 and so , 0) = 0 and so If R 0 > 1, then system (2.1) has a positive equilibrium point In the following, we show that P * is the unique positive equilibrium point of system (2.1). For any positive equilibrium point P * , by (2.2) and hypotheses (H1) and (H2), we have Since m 1 m 4 = (δ 0 + μ)(δ 0 + η) > μη and hypotheses (H1) and (H2) hold, we have Since g(0) = 0 and g is differentiable on R + , there exists ξ ∈ (0, I * ) such that g(I * ) I * = g (ξ ). By using hypotheses (H2) and (H3), one has Thus, it follows from (2.3), (2.4) and (2.5) that F (I * ) < 0. Suppose that there exists another positive equilibrium point P 1 = (S 1 , E 1 , I 1 , R 1 , V 1 ). Then F (I 1 ) ≥ 0 due to the property of continuous function. This is a contradiction. Therefore, system (2.1) has a unique endemic equilibrium P * when R 0 > 1. It can be stated as follows.
3) has a disease-free equilibrium P 0 as follows: which exists for all parameter values. For R 0 > 1, the endemic equilibrium P * admits the unique positive equilibrium point for system (1.3).
Remark 2.2 From the proof of the existence of endemic equilibrium P * , it is not difficult to arrive at such a conclusion that the nonlinear treatment function g(I) has an upper bound Aδ 1 m 2 , which is reasonable for limited medical resources in our daily life. Proof Summing up the five equations in system (1.3) and denoting

Proposition 2.1 The set
we get Now integrating both sides of the above inequality and using the theory of differential inequality [41], we obtain Thus, the set Ω is positiveinvariant, that is, all initial solutions belong to Ω remain in Ω for all t > 0.

Global stability of the disease-free equilibrium by means of Lyapunov function
In this section, we investigate the global stability of the disease-free equilibrium P 0 for system (1.3).
A(m 3 +g (0)) , then the disease-free equilibrium P 0 of system (1.3) is globally asymptotically stable in the feasible region Ω. If R 0 > 1, then P 0 is unstable.
Proof The Jacobian matrix of system (1.3) at P 0 is Obviously, λ 1 = -δ 0 is an eigenvalue of J(P 0 ). The other eigenvalues of J(P 0 ) are determined by the equations A(m 3 +g (0)) , is the singleton P 0 . By the LaSalle invariance principle theorem ( [42], p. 30), the disease-free equilibrium P 0 is globally asymptotically stable (0)) . This completes the proof.

Global stability of the endemic equilibrium by means of geometric approach
In this section, we analyze the stability of the endemic equilibrium P * . First, we show the local stability of the endemic equilibrium of system (1.3) around the endemic equilibrium P * .
Here all the parameters a 11 , a 13 , a 21 , a 33 , a 43 are defined in (4.1). The values of h and z equal to B 1 and B 3 , respectively.
Proof The Jacobian matrix of system (1.3) at P * is given by  Using assumptions (i) and (ii), by a direct calculation, we have B i > 0 for i = 1, 2, 3, 4. It follows from the Routh-Hurwitz criteria [43] that all the eigenvalues associated to J(P * ) have negative real parts iff B i > 0 for i = 1, 2, 3, 4 and To find the global stability of system (1.3), it is necessary to reduce system (1.3) first. Since recovered class R does not have any effect on the dynamics of S, V , E and I class, we shall investigate the following system:  The solutions of (4.2) corresponding to nonnegative initial values remain nonnegative for all time. Moreover, we observe that the total population size of (4.2) denoted by X(t) satisfiesẊ = Aδ 0 Xδ 2 Iδ 3 Ig(I), so that we can study the model in the region: Here we follow the approach used in [8] for a SVEIR model of SARS epidemic spread. Let us consider the following autonomous dynamical system: where f : D → R n , D ⊂ R n which is an open set, simply connected and f ∈ C 1 (D). Suppose that x * is an equilibrium point of (4.3), i.e. f (x * ) = 0. Therefore, x * is said to be globally stable in D if it is locally stable and all trajectories in D converge to x * . Let Q(x) be a matrix-valued function of order n 2 × n 2 that is C 1 on D. We also consider the matrix A which is defined as and the matrix M is the second additive compound matrix of the Jacobian matrix J. Further the Lozinskiȋ measureμ of A with respect to a vector norm · can be defined in R n 2 as follows: We will apply the following theorem according to [44]. Theorem 2.1 states that R 0 > 1 implies the existence and uniqueness of the endemic equilibrium P * . Further, we know that the disease-free equilibrium P 0 is unstable when R 0 > 1. The instability of P 0 , together with P 0 ∈ ∂Θ, which implies the uniform persistence of the state variables (see [45]). Thus, there exists a constant a > 0 such that any solution (S(t), E(t), I(t), V (t)) with (S(0), E(0), I(0), V (0)) in the orbit of system (4.2) satisfies The uniform persistence of system (4.2), incorporating the boundedness of Θ, suggests that the compact absorbing set in the interior of Θ; see [46]. Hence, Lemma 4.1 may be applied, with D = Θ.
According to [47], the Lozinskiȋ measure in Lemma 4.1 can be evaluated as: where D + is the right-hand derivative. The endemic equilibrium is locally asymptotically stable, provided R 0 > 1. Hence, to get the global asymptotic stability, according to Lemma 4.1, the trick of the proof is to find a norm · such thatμ(A) < 0 for all x in the interior of Θ.
Furthermore, we assume that We will use the inequalities mentioned above to get the estimates on D + z .
Proof The basic idea of the proof is to obtain the estimate of the right derivate D + z of the norm (4.5). For this purpose, we need to discuss sixteen cases according to the different orthants and the definition of the norm (4.5) within each orthant. Case 1: This shows that By using |z 5 | < U 2 < |z 1 | and |z 2 | + |z 3 | < |z 1 |, it follows from (4.8) that Case 2: U 1 > U 2 and z 1 , z 2 , z 3 > 0 with |z 1 | < |z 2 | + |z 3 |. Then, z = |z 2 | + |z 3 |. (4.9) Thus, we have Using |z 6 | < U 2 < |z 2 | + |z 3 | and |z 1 | < |z 2 | + |z 3 |, in view of (4.9), one has The discussion for the other fourteen cases are similar to the ones discussed in [7] and so we omit it here. Thus, we can get the following estimate: Now the global stability follows from Lemma 4.1.
Remark 4.1 As pointed out by Buonomo and Lacitignola [7], in some real situations, different choices of the matrix Q and of the vector norm · may lead to better sufficient conditions than those we presented here, in the sense that the assumptions on the parameters may be weakened. Thus, it is worth to note that sufficient conditions (4.6) and (4.7) in Theorem 4.2 are derived from the application of the method and numerical simulations suggest that they may be not necessary (see Example 5.1).

Numerical simulations
The aim of this section is to give a numerical example to illustrate our main results.
Example 5.1 Consider the system We first consider the case when by using the parameter values given in Table 1. Using these parameter values, for different initial conditions the dynamics of model (5.1) is presented in Figs. 1-5. It shows that system (5.1) has a disease-free equilibrium and it is globally asymptotically stable. This numerical verification supports the result stated in Theorem 3.1.
Next, we consider the case when R 0 = 2.211436 > 1 using the parameter values given in Table 2. Using these parameter values, for different initial conditions the dynamics of model (5.1) is presented in Figs. 6-10. It shows that system (5.1) has an endemic equilibrium and it is globally asymptotically stable with different initial values, which supports our analytical results stated in Theorem 4.2.   Table 1 Figure 2 Time series plot of the exposed population for R 0 = 0.571429 < 0.910714 with various initial conditions, parameter values are given in Table 1 Figure 3 Time series plot of the infective population for R 0 = 0.571429 < 0.910714 with various initial conditions, parameter values are given in Table 1 Figure 4 Time series plot of the recovered population for R 0 = 0.571429 < 0.910714 with various initial conditions, parameter values are given in Table 1 Figure 5 Time series plot of the vaccinated population for R 0 = 0.571429 < 0.910714 with various initial conditions, parameter values are given in Table 1 Table 2 Figure 7 Time series plot of the exposed population for R 0 = 2.211436 > 1 with various initial conditions, parameter values are given in Table 2 Figure 8 Time series plot of the infective population for R 0 = 2.211436 > 1 with various initial conditions, parameter values are given in Table 2 Figure 9 Time series plot of the recovered population for R 0 = 2.211436 > 1 with various initial conditions, parameter values are given in Table 2 Figure 10 Time series plot of the vaccinated population for R 0 = 2.211436 > 1 with various initial conditions, parameter values are given in Table 2 6 Conclusions In this paper, we have considered an SVEIR epidemic model with general nonlinear incidence rate. In model (1.3), we have divided the total population into five compartments, namely susceptible, exposed, infective, recovered, vaccinated population and investigated the dynamical behavior of this model. Here, we have found that is a basic reproduction number of system (1.3), which helps us to determine the dynamical behavior of the system. We have showed that system (1.3) to be globally asymptotically stable at disease-free equilibrium P 0 when 0)) .
When R 0 > 1, the endemic equilibrium stable both locally and globally has been derived and analyzed under some conditions. The important mathematical findings for the dynamical behavior of model (1.3) have also numerically been verified for a special case of model (1.3). We would like to point out that the model considered in this paper is not a case study and so it is difficult to choose parameter values from quantitative estimation. We have used hypothetical sets of parameters to verify our analytical results. It is worth to mention that the results presented in this paper improve and extend some related results in [9,10,12,13]. Finally, we remark that there are quite a few spaces to deserve further investigation. For example, we can continue the research in this line considering the vaccination rate μ in our model (1.3) as a continuous function, and, later, a discontinuous function. On the other hand, as is well known, epidemiological models which incorporate the control strategies can be useful to both control the spread of disease and minimize the intervention costs. For our model, it is natural to consider vaccination rate coefficient as a control to reduce the disease burden. Thus, it is important and interesting to prove the existence of optimal control, characterize the optimal control, prove the uniqueness of optimal control, compute the optimal control numerically and investigate how the optimal control depends on various parameters in the models. We will devote to these questions our future work.