Global Stability Analysis of SEIR Model with Holling Type II Incidence Function

A deterministic model for the transmission dynamics of a communicable disease is developed and rigorously analysed. The model, consisting of five mutually exclusive compartments representing the human dynamics, has a globally asymptotically stable disease-free equilibrium (DFE) whenever a certain epidemiological threshold, known as the basic reproduction number (ℛ 0), is less than unity; in such a case the endemic equilibrium does not exist. On the other hand, when the reproduction number is greater than unity, it is shown, using nonlinear Lyapunov function of Goh-Volterra type, in conjunction with the LaSalle's invariance principle, that the unique endemic equilibrium of the model is globally asymptotically stable under certain conditions. Furthermore, the disease is shown to be uniformly persistent whenever ℛ 0 > 1.


Introduction
Mathematical models have been widely used to gain insight into the spread and control of emerging and reemerging disease. The dynamics of these models is usually determined by a threshold quantity known as the basic reproduction number (denoted by R 0 ), which is defined as the number of secondary cases generated by an infected individual in a completely susceptible population [1][2][3][4][5]. Characteristically, when R 0 is less than unity, a small influx of infected individuals will not generate large outbreaks, and the disease dies out in time. On the other hand, when R 0 exceeds unity, the disease will persist. A basic epidemic model supports at least two equilibria (a disease-free equilibria and endemic equilibria); R 0 plays important role in the study of equilibria of a model. Several models found in the literature [2,[4][5][6][7][8][9][10][11][12][13][14][15] have been used to show that when R 0 crosses the threshold, R 0 = 1, a transcritical bifurcation takes place. That is, asymptotic local stability is transferred from the disease-free state to the new (emerging) endemic (positive) equilibria. In some cases, it can be shown that the transfer of asymptotic stability is independent of initial conditions; that is, it is global; see, for instance, [6,[8][9][10][11][12]16]. Establishing global properties of a dynamical system using Lyapunov function is generally a nontrivial problem. This is owing to the fact that there are no systematic methods for constructing Lyapunov function for infectious disease models with standard incidence rate [17]. The most successful approach to the problem is the direct Lyapunov method which involves the use of quadratic function of the form ω(x 1 , x 2 , . . . , x n ) = n i=1 (c i /2)(x i − x * i ) or by using nonlinear Lyapunov function of Goh-Volterra type of the form ω(x 1 , x 2 , . . . , ). However, other methods used in establishing global properties of some epidemic models include Dulac's criterion to eliminate the existence of the periodic solution and prove the global stability by the Poincaré Bendixson theorem [18] and those reported in Kamgang and Sallet [19] and Qiao et al. [20].
Let S(t), I(t), and N(t) denote the number of susceptible individuals, infectious individuals, and the total size of the population at time t, respectively. Further, let β(N) be the average number of contacts that is sufficient to transmit infection (effective contact rate). Then, the force of infection, given by β(N)I/N, represents the average number of contacts a susceptible individual makes with infectious individuals per unit time. If β(N) = βN (i.e., the contact rate depends on the total population, N), then the incidence function g 1 (I) = βI is called mass action incidence. If β(N) = β (a constant), then the incidence function g 2 (I) = βI/N is called standard incidence [4]. These two functions are widely used in modeling the transmission dynamics of the human diseases [1,21]. Another widely used incidence function is the Holling type II incidence function, given by g 3 (I) = βI/(1 + ωI), with ω > 0, [22][23][24][25][26][27][28].
The nonlinear incidence function of type g 3 (I) was first introduced by Capasso and Serio [22], in their study of cholera epidemic. The main justification for using such a functional form of the incidence function stems from the fact that the number of effective contacts between infective individuals and susceptible individuals may saturate at high infective levels due to crowding of infective individuals or due to the preventive measures (and behavioral changes) taken by the susceptible individuals in response to the severity of the disease [23,[25][26][27].
A number of mathematical models have been developed in the literature to gain insights into the transmission dynamics of diseases with subpopulation (compartments) [4,[6][7][8][9][10][11][12][13][14][29][30][31]. The choice of which compartment to include in a model depends on the characteristic of the particular disease being modelled and the purpose of the model [4].
The classical SIR or SIRS model assumed that the disease incubation is negligible so that each susceptible individual (in the S class) once infected becomes infectious (and move to I class) and later recovers (to move to R class) where they acquire permanent or temporary immunity [32]. However, more general models than SIR or SIRS models assumed that susceptible individuals, once infected, first go through the latent period (in the E class) before becoming infectious; the resulting models are of SEIR or SEIRS type, depending on whether recovered individuals acquired permanent or temporary immunity [33,34].
The model considered in this study is based on SEIR (S: susceptible, E: exposed, I: infected, R: recovered) where recovered individuals acquire permanent immunity, so that they will not become infected again. This is owing to the fact that for many viral diseases such as measles, smallpox, rubella, HIV/AIDS, and influenza recovered individuals confer lifelong immunity [29,33]. Huang and Takeuchi [29] study the classical SIR, SIS, SEIR, and SEI models with time delay and a general incidence function. Safi et. al. [13] consider the effect of periodic fluctuations on the transmission dynamics of a communicable disease using SEIRS model, subject to quarantine and isolation; the authors show that adding periodicity to the autonomous model does not alter the threshold dynamics of the model with respect to the control of the disease in the population. Zhang and Ma [35] considered the global dynamics of an SEIR model with saturating contact rate. Korobeinikov [36] established global asymptotic dynamics of SEIR and SIR models with several parallel infectious stages. Li et al. [33] analyzed the global dynamics of an SEIR model with vertical transmission and bilinear incidence. Li and Jin [34] consider global stability of an SEIR epidemic model with infectious force in latent, infected, and immune period. It should be stated, however, that the aforementioned three studies [33,34,36] considered mass action (bilinear) incidence to model the infection. This study focuses on the mathematical modeling of the transmission dynamics of an arbitrary disease with educated (counsel) and uneducated infectious stages.
This study has three important differences from those reported in [33,34,36]. The first is established essential qualitative features of an SEIR model with Holling type II incidence function. The second is that a public health education and counselling are offered to infected individuals (I u ). Recent studies further reinforced the widely held belief that one key strategy for preventing and controlling the spread of communicable disease such as HIV (especially in resource-poor nations) is to provide HIV-related public health education and counselling (such as sexual education and awareness of the risk and life-threatening consequences of HIV/AIDS) which would, hopefully, lead to reduction in risky sexual behavior and safer lifestyle within the community [7]. In addition to the above extensions, rigorous qualitative analysis will be provided for the resulting SEIR model. In particular, the paper gives special emphasis on the global asymptotic stability of the disease-free equilibrium and endemic equilibrium.
The paper is organized as follows. The model is formulated and analysed for its basic dynamical features in Section 2. The stability analysis of the model is presented in Sections 3 and 4.

Model Formulation
The total population at time t, denoted by N(t), is subdivided into five compartments of susceptible (S(t)), exposed (those who have been infected but are not yet infectious) (E(t)), uneducated infected individuals (I u (t)), educated infected individuals (I e (t)), and recovered (R(t)) individuals, so that The susceptible population is increased by the recruitment of individuals into the population (assumed susceptible), at a rate Π. Susceptible individuals may acquire infection, following effective contact with infected individuals (in the I u or I e class) at a rate λ(t), where In (2), β is the effective contact rate (contact capable of leading to infection), while the modification parameter 0 < η < 1 accounts for the assumed reduction in disease transmission by educated infected individuals in comparison to uneducated infected individuals in the I u class. The population of susceptible individuals is further decreased by natural death (at a rate μ). Thus, the rate of change of the susceptible population is given by The population of exposed individuals is generated by the infection of susceptible individuals (at the rate λ(t)). This Computational and Mathematical Methods in Medicine 3 population is decreased by development of disease symptoms (at a rate κ) and natural death (at a rate μ), so that The population of uneducated infected individuals is generated at the rate κ. It is decreased by natural recovery (at a rate γ 1 ), education (at a rate σ), natural death (at the rate μ), and disease-induced death (at a rate δ 1 ). This gives The population of educated infected individuals is generated by the education of infected individuals (at the rate σ). This population is decreased by recovery (at a rate γ 2 ), natural death (at the rate μ), and disease-induced death (at a rate δ 2 < δ 1 ). It is assumed that the disease-induced mortality rate of educated infected individuals is low in comparison with uneducated infected individuals. Hence, the rate of change of this population is given by Finally, the population of recovered individuals is generated by the recovery of uneducated and educated infected individuals (at the rates γ 1 and γ 2 , resp.). It is decreased by natural death (at the rate μ), so that Thus, the model for the transmission dynamics of an infectious disease in the presence of educated (counsel) infected individuals is given by the following nonlinear system of differential equations: The model (8) extends some SEIR models in the literature such as those in [33,34,36] by (i) replacing the mass action incidence function with Holling type II incidence function, (ii) splitting the compartment of infected individuals into educated (counsel) and uneducated (noncounsel) infected individuals thereby allowing time-varying infection rate. The epidemiological implication of this assumption is that educated infected individuals transmit the disease at a reduced rate (0 < η < 1) in comparison to the uneducated infected individuals (due to behavioral changes). For instance, In Zambia, the decline in HIV incidence since early 1990s is attributed to behavioral changes [37]. Public health education campaigns have also been successfully implemented in numerous countries and communities, such as Uganda, Thailand, Zambia, and the US gay community [38,39]. Between 1991 and 1998, HIV prevalence dramatically declined in Uganda from 21 to 9.8% (with a corresponding reduction in nonregular sexual partners by 65% coupled with greater levels of awareness about HIV/AIDS [38]. The Ugandan programme fostered community mobilization towards change in risky behavior, without increasing stigma [40]. Further, unlike in the aforementioned modelling studies, detailed rigorous mathematical analysis of the model (8) will be provided.

Basic Properties.
Since the model (8) which can be rewritten as Hence, so that Similarly, it can be shown that E > 0, I u > 0, I e > 0 and R > 0, for all time t > 0.

Computational and Mathematical Methods in Medicine
The previous result can also be established using the method in of [41,Appendix A].
We claim the following result.

Lemma 2. The closed set
is positively invariant.
Proof. Adding all the equations of the model (8) gives, Since Since the region D is positively invariant, it is sufficient to consider the dynamics of the flow generated by the model (8) in D, where the usual existence, uniqueness, continuation results hold for the system [10].

Local Stability of Disease-Free Equilibrium (DFE)
The DFE of the model (8) is given by The local stability of E 0 will be explored using the next generation operator method [43,44]. Using the notation in [44], the nonnegative matrix, F, of the new infection terms and the M-matrix, V , of the transition terms associated with the model (8) are given, respectively, by It follows that the control reproduction number [4,21], denoted by R 0 = ρ(FV −1 ), where ρ is the spectral radius, is given by where Using in [44,Theorem 2], the following result is established.
The quantity R 0 measures the average number of new infections generated by a single infected individual in a population. Lemma 3 implies that the disease can be eliminated from the community (when R 0 < 1) if the initial sizes of the subpopulations of the model are in the basin of attraction of the DFE (E 0 ). To ensure that disease elimination is independent of the initial sizes of sub-populations, it is necessary to show that the DFE is globally asymptotically stable (GAS) if R 0 < 1. This is explored below. (8), given by (15), is GAS in D whenever R 0 ≤ 1.

Theorem 4. The DFE of the model
Proof. Consider the following Lyapunov function: with Lyapunov derivative (where a dot represents differentiation with respect to time) given bẏ Since all the parameters and variables of the model (8)  {(E, I u , I e ) = (0, 0, 0)}. Thus, it follows, by the LaSalle's invariance principle [18,45], that (E, I u , I e ) −→ (0, 0, 0) as t −→ ∞. (21) Since lim sup t → ∞ I u = 0 and lim sup t → ∞ I e = 0 (from (21)), it follows that, for sufficiently small * > 0, there exist constants M 1 > 0 and M 2 > 0 such that lim sup t → ∞ I u ≤ * for all t > M 1 and lim sup t → ∞ I e ≤ * for all t > M 2 . Hence, it follows from the last equation of the model (8) that, for t > max{M 1 , Thus, by comparison theorem [42], so that, by letting * → 0, Similarly (by using lim inf t → ∞ I u = 0 and lim inf t → ∞ I e = 0), it can be shown that Thus, it follows from (24) and (25) that Hence, Similarly, it can be shown that Thus, by combining (21), (27), and (28), it follows that every solution of the equations of the model (8), with initial conditions in D, approaches E 0 as t → ∞ (for R 0 < 1).
The previous result shows that the disease can be eliminated from the community if the associated reproduction number of the model is less than unity.

Existence and Stability for Endemic Equilibrium Point
In this section, the possible existence and stability of endemic (positive) equilibria of the model (8) (i.e., equilibria where at least one of the infected components of the model is nonzero) will be explored.

Existence of Endemic Equilibrium Point (EEP).
First of all, the persistence of the disease in the population will be investigated below. The model system (8) is said to be uniformly persistent if there exists a constant c such that any solution (S(t), E(t), I u (t), I e (t), R(t)) satisfies ( [31,46]) provided that (S(0), E(0), I u (0), I e (0), R(0)) ∈ D • (the interior of the region D).
Proof. The theorem can be proved by using the approach used to prove Proposition 3.3 of [32], by applying a uniform persistence result in [46] and noting that the DFE of the model (8) is unstable whenever R 0 > 1 (Lemma 3).
When R 0 > 1, it follows (from Theorem (3)) that model (8) is uniformly persistent; by using in [47,Theorem 2.8.6] and in [42,Theorem D.3] it follows that the model (8) has at least one endemic equilibrium in D • . Thus the following result is established.
The uniqueness of the endemic equilibrium will be investigated at the end of the next subsection.

Global Stability of Endemic Equilibrium for Special Case.
Here, the global asymptotic stability property of the endemic equilibrium of the model (8) is given for a special case when educated (counsel) infected individuals do not transmit infection (η = 0). The model (8), with η = 0, then reduces to: where The reproduction number of the model (30) is given by We claim the following result.

6
Computational and Mathematical Methods in Medicine Theorem 7. The endemic equilibrium of the reduced model, given by (30), is GAS in D • if R 0r > 1. Thus it is unique.
Proof. Consider the reduced model, given by (30). Let R 0r > 1, so that the associated endemic equilibrium exists. Further, consider the following nonlinear Lyapunov function: with Lyapunov derivativė It can be shown from (30) that, at endemic steady state, Using the relations (35) in (34) giveṡ Adding and subtracting βS * * f (I * * u ) and βS * * f 2 (I * * u )I u / f (I u )I * * u in (36) giveṡ