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 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 R 0 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 R 0 as well as on the treatment rate. If R 0 > 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-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–14, 20–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–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 and μ is birth rate as well as the death rate of the susceptible class. The nonlinear function denotes the sum of the death rate and the recovery rate for the infected individuals satisfying and for . 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: with initial conditions: , , , and . 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, 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, is the virus induced crashing rate and , 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 , , , , are the number of susceptible population, exposed population, infective population, recovered population, vaccinated population, respectively. The two-variable function represents incidence rate and the nonlinear function denotes the removal rate of infective individuals because of the treatment of infective. The initial conditions for system (1.3) are as follows: Clearly, denotes the total number of high-risk human population at time t.

The model parameters of system (1.3) are described as follows:

the rate at which new individuals (including newborns and immigrants) enter the susceptible population,

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,

the rate at which exposed population become infective,

natural recovery rate of infective population,

disease-related death rate of infective population.

Model (1.3) involves certain assumptions which consist of the following:

The new individuals enter the population with a constant rate and all the new individuals are susceptible.

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.

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.

After recovery the individuals become immunized and hence they are no longer susceptible to it.

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 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 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 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 and indicate the existence of equilibrium, we give some hypotheses.

is differentiable such that

for all ;

for all ;

for all and ;

for all ;

for all .

is differentiable such that , and for .

holds for all , where for .

Remark 2.1

It is easy to check that the classes of satisfying (H1) include incidence rates such as for and .

It is straightforward to show that the classes of satisfying (H2) include removal rates such as for , , .

By hypothesis (H2), we know that is a monotone decreasing function on .

The assumption (H3) is equivalent to the following inequality: which can be found in [27].

By the assumptions, it is easy to find that system (1.3) always has a disease-free equilibrium point , where

We shall assume that (H1), (H2) and (H3) hold in the rest of this paper.

Now we define the basic reproduction number for model (1.3) as follows: where

In order to find the positive equilibria of model (1.3), set It follows that and .

Substituting the above equalities into the second equation in (2.1), one has Let Then it is easy to see that the positive equilibrium points of system (2.1) are given by zeros of F in the interval .

We denote for convenience. Then it is easy to see Therefore, we conclude that has at least one root named Ĩ in the interval . That is, and so Since , we know that and so If , then system (2.1) has a positive equilibrium point , where

In the following, we show that is the unique positive equilibrium point of system (2.1). For any positive equilibrium point , by (2.2) and hypotheses (H1) and (H2), we have

Since and hypotheses (H1) and (H2) hold, we have Since and g is differentiable on , there exists such that . By using hypotheses (H2) and (H3), one has Thus, it follows from (2.3), (2.4) and (2.5) that .

Suppose that there exists another positive equilibrium point . Then due to the property of continuous function. This is a contradiction. Therefore, system (2.1) has a unique endemic equilibrium when . It can be stated as follows.

Theorem 2.1

System (1.3) has a disease-free equilibrium as follows: which exists for all parameter values. For , the endemic equilibrium admits the unique positive equilibrium point for system (1.3).

Remark 2.2

From the proof of the existence of endemic equilibrium , it is not difficult to arrive at such a conclusion that the nonlinear treatment function has an upper bound , which is reasonable for limited medical resources in our daily life.

Proposition 2.1

The set is a positively invariant and attracting region for the disease transmission model given by system (1.3) with initial conditions (1.4).

Proof

Summing up the five equations in system (1.3) and denoting we get i.e., Now integrating both sides of the above inequality and using the theory of differential inequality [41], we obtain Clearly, as . If , then . Thus, the set Ω is positive-invariant, that is, all initial solutions belong to Ω remain in Ω for all . □

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 for system (1.3).

Theorem 3.1

If , then the disease-free equilibrium of system (1.3) is globally asymptotically stable in the feasible region Ω. If , then is unstable.

The Jacobian matrix of system (1.3) at is Obviously, is an eigenvalue of . The other eigenvalues of are determined by the equations and respectively. If , then one eigenvalue is positive. Thus, is unstable.

When , to prove the disease-free equilibrium is globally asymptotically stable, we consider the Lyapunov function . The derivative of along system (1.3) satisfies Furthermore, iff . Thus, the largest compact invariant set in , when , is the singleton . By the LaSalle invariance principle theorem ([42], p. 30), the disease-free equilibrium is globally asymptotically stable if . 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 . First, we show the local stability of the endemic equilibrium of system (1.3) around the endemic equilibrium .

Theorem 4.1

If , then the endemic equilibrium exists and is locally asymptotically stable if the following conditions hold:

;

;

and . Here all the parameters , , , , are defined in (4.1). The values of h and z equal to and , respectively.

The Jacobian matrix of system (1.3) at is given by where Clearly, one of the roots of is negative, i.e. . The remaining roots can be determined from the following equation: where Using assumptions (i) and (ii), by a direct calculation, we have for . It follows from the Routh–Hurwitz criteria [43] that all the eigenvalues associated to have negative real parts iff for and .

Now, if (iii) holds. This ends the proof. □

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 satisfies , 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 , which is an open set, simply connected and . Suppose that is an equilibrium point of (4.3), i.e. . Therefore, is said to be globally stable in D if it is locally stable and all trajectories in D converge to .

Let be a matrix-valued function of order that is on D. We also consider the matrix A which is defined as where the matrix is 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 as follows:

We will apply the following theorem according to [44].

Lemma 4.1

If is a compact absorbing subset in the interior of D, and there exist and a Lozinskiĭ measure for all , then every omega limit point of system (4.2) in the interior of D is an equilibrium in .

Theorem 2.1 states that implies the existence and uniqueness of the endemic equilibrium . Further, we know that the disease-free equilibrium is unstable when . The instability of , together with , which implies the uniform persistence of the state variables (see [45]). Thus, there exists a constant such that any solution with 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 .

According to [47], the Lozinskiĭ measure in Lemma 4.1 can be evaluated as: where is the right-hand derivative. The endemic equilibrium is locally asymptotically stable, provided . Hence, to get the global asymptotic stability, according to Lemma 4.1, the trick of the proof is to find a norm such that for all x in the interior of Θ.

Starting with the Jacobian matrix J of (4.2), the second additive compound matrix is given by Hence, the second additive compound matrix of J is given as follows: where Now we consider the following matrix: Then we obtain the matrix , where is the derivative of Q in the direction of the vector field f. More accurately, we have Hence, in view of the fact that we obtain where As in [8], we consider the following norm on : where , with components , , and is defined as: and is defined as In the next, we will use the following inequalities: and Furthermore, we assume that We will use the inequalities mentioned above to get the estimates on .

Theorem 4.2

For , system (4.2) admits an unique endemic equilibrium which is globally asymptotically stable in the interior of Θ, provided that inequality (4.6) is satisfied and that for some positive constant ω, where

The basic idea of the proof is to obtain the estimate of the right derivate 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 and , it follows from (4.8) that

Case 2: Thus, we have Using and , 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: where 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 which is a particular case of system (1.3) by letting and , where m, n, γ, a are positive and . The other parameters in (5.1) have the same biological meanings as in model (1.3).

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.

Table 1

Parameter for Figs. 1–5

Parameter	Values
A	2
δ_0	0.2
m	0.2
n	2
η	0.2
μ	0.4
δ_1	0.8
δ_2	0.5
δ_3	0.55
γ	0.3
a	2

Figure 1

Time series plot of the susceptible population for with various initial conditions, parameter values are given in Table 1

Figure 5

Time series plot of the vaccinated population for with various initial conditions, parameter values are given in Table 1

Time series plot of the exposed population for with various initial conditions, parameter values are given in Table 1

Time series plot of the infective population for with various initial conditions, parameter values are given in Table 1

Time series plot of the recovered population for with various initial conditions, parameter values are given in Table 1

Next, we consider the case when 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 2

Parameter for Figs. 6–10

Parameter	Values
A	6
δ_0	0.5
m	0.8
n	2
η	0.2
μ	0.4
δ_1	0.8
δ_2	0.5
δ_3	0.55
γ	0.3
a	2

Figure 6

Time series plot of the susceptible population for with various initial conditions, parameter values are given in Table 2

Figure 10

Time series plot of the vaccinated population for with various initial conditions, parameter values are given in Table 2

Time series plot of the exposed population for with various initial conditions, parameter values are given in Table 2

Time series plot of the infective population for with various initial conditions, parameter values are given in Table 2

Time series plot of the recovered population for with various initial conditions, parameter values are given in Table 2

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 when When , 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.