Global Dynamics of an SEIRS Compartmental Measles Model with Interrupted Vaccination

Dominic Otoo, Justice A. Kessie, Elvis K. Donkoh, Eric Okyere, Williams Kumi

Department of Mathematics & Statistics, University of Energy & Natural Resources, Sunyani, Ghana.

DOI: 10.4236/am.2019.107042                

Abstract

Measles is a reemerging disease that has a devastating impact, especially among children under 5. In this paper, an SEIRS model is developed to investigate a possible outbreak among the population of children under 5 in the Sunyani Municipality. We consider waning immunity or loss of immunity among those who were vaccinated, which leads to secondary attacks among some in the population. Using Routh-Hurwitz criterion, Matrix Theoretic and Goh-Volterra Lyapunov functions, the stability of the model was investigated around the equilibria. We have computed the threshold parameter, R₀, using the Next Generation Matrix method. The disease-free equilibrium is globally stable whenever R₀ ≤1 and unstable otherwise. The endemic equilibrium is globally stable when R₀ >1.

Keywords

Global Stability, Lyapunov Function, Matrix Theoretic Method, Next Generation Matrix, SEIRS

1. Introduction

Measles, a recurrent virus infection has a short term outbreak but its impact is devastating especially among children under five. Severe measles results in pneumonia, which is the disease with high mortality rate in Ghana. Found in [1], the graphs show the mortality during an outbreak of measles. In this research paper, the authors looked at a case study of the dynamics of measles disease in the case of interrupted vaccination in Sunyani, a thriving metropolis in Ghana, a country where only 57% deliveries are handled by qualified health personnel [2]. In the Sunyani city and its immediate environs, new deliveries are commonly done in the hospital, clinics and maternity homes but significantly some cases of new births are successfully done without the presence of qualified trained health practitioner, usually in their homes. The health program in Ghana makes it necessary for regular attendant mothers with newborns to vaccinate their children against the six childhood killer diseases including measles in any health facility. The first dose of MMR (mumps, measles, rubella) vaccine is given at 12 - 13 months of age. The second dose is usually taken when the child is 3 years and 4 months old. This can be taken in later years in the management of a disease outbreak. The irregular attendants do not successfully vaccinate their children against measles. New deliveries done outside a health facility make it difficult to track these children for vaccination. There are some children who missed the first and second doses, and there are some still who also missed the second dose because of irregularity in attending a health facility and for immigration, too [3]. There is also among the population a strong believe that MMR vaccines are leaving children with autism. Thereby, there are many who do not have a strong immunity against measles and many have their immunity wane in time. Their development of a strong immunity is interrupted so several children are susceptible to the disease and several still may have secondary infection. We seek to investigate the nature of a possible outbreak of measles in a population with dynamics such as this. Many researchers have studied the global dynamics of infectious diseases in general (see [4] [5]) while others also studied the dynamics of specific infectious diseases such as measles, chickenpox, mumps and rubella (see [6] [7]). In this paper, we also study the dynamics of measles with a relapse in a varying population where birth, death and immigration dynamics are considered. [8] [9] [10] studied the global stability of disease models and inculcated immigration and/or births and death dynamics into the population system. Again much work has already been done (see [11] [12]) defining what measles is and the symptoms and the effect of chronic infection of measles. It is mentioned that measles is severe among children under five (5) years.

Mathematical models are strong instruments used extensively to study the spread and control of infectious diseases. One important measure that determines the dynamics of disease models is the threshold parameter known as basic reproductive number $R_0$ . This parameter measures the number of infectives generated by a single infectious individual introduced into the susceptible population. In measuring the foregoing, researchers in recent literature use the NGM approach. When $R_0\le 1$ (usually when the disease-free equilibrium exists) the introduction of infectious individual can not generate a large outbreak. On the other hand, when $R_0>1$ , (where the endemic equilibrium exists), the disease will persist in the population. [4] [13] have mentioned in their works that establishing global properties of dynamical systems using Lyapunov functions is not a trivial problem because there is no systematic procedure in constructing Lyapunov functions for infectious diseases. But one of the most successful procedures used to construct Lyapunov functions by many researchers is the quadratic function of the form ${\mathcal{D}}_i={\displaystyle {\sum }_{i=1}^n{c}_i{\left(x-x^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}$ (see [14]) or the nonlinear Goh-Volterra function ${\mathcal{D}}_i={\displaystyle {\sum }_{i=1}^n{c}_i\left(x-x^{*}-x^{*}\ln \frac{x}{x^{*}}\right)}$ (see [15] [16]).

In this paper, the matrix theoretic method, a type of Lyapunov function, is used to establish the stability in the disease free system. Later in this study, the Goh-Volterra method is used to establish the stability of the endemic equilibrium. [4] [17] used this matrix theoretic method to study the global dynamics of several specific disease models.

2. The Model

We consider a susceptible-exposed-infectious-recovered-susceptible (SEIRS) model with relapse since there is a vaccination parameter that provides a near-permanent immunity against the disease with non-negative initial conditions. E represents the number of latent(exposed) individuals, $\sigma >0$ represents the rate at which the exposed individuals become infectious (i.e.) $1/\sigma$ represents the average latent period. $\lambda >0$ represents the rate at which the susceptible individual become immune and gain near-permanent immunity and $\gamma >0$ also represents the rate at which the infectious individuals go through the successful treatment process of the disease. Near-permanent immunity because there exist a rate $ϵ>0$ of the recovered individuals falling back into the susceptible compartment. It is considered that this $ϵR$ individual have a waning immunity over time so they fall back susceptible to the infectious disease. Already, [18] studied the “Global Stability Results in an SVIR Epidemic Model with Immunity Loss Rate Depending on the Vaccine-Age” where it is assumed that the waning rate of vaccine-induced immunity is depending on the vaccine-age. $\alpha$ represents a vital dynamics-birth and immigration rates. The total number of new infections at a time t is given by $\beta SI/N$ . It is also considered that the natural immunity of some exposed individuals ${\delta }_E$ moves them into the recovered population. While $\mu$ is the natural death rate for the susceptible and recovered compartments, ${\mu }_E=\mu +k_E$ and $\tau =\mu +k_I$ represent the death rates for the exposed and infected classes respectively. $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right)$ is a growing population since births and immigration rates are added at any time t. The schematic diagram as shown in Figure 1 is drawn or developed to depict the interactions in the thriving metropolis, Sunyani of Ghana. The schematic diagram above in Figure 1 has the following as its system of equations;

$\begin{array}{l}{S}^{\prime }=\alpha N+ϵR-\left(\mu +\beta I/N+\lambda \right)S\\ E^{\prime }=\beta SI/N-\left({\mu }_E+\sigma +{\delta }_E\right)E\\ I^{\prime }=\sigma E-\left(\tau +\gamma \right)I\\ R^{\prime }=\gamma I-\left(\mu +ϵ\right)R+\lambda S+{\delta }_{E}E\end{array}$ (1)

3. Model Transformation & Analysis

The total population size $N_t$ can be determined by $N=S\left(t\right)+E\left(t\right)+I\left(t\right)+R(t)$

Figure 1. Schematic diagram of the SEIRS model.

or from the differential equation $N^{\prime }=\left(\alpha -\mu \right)N-k_{E}E-k_{I}I$ , which is derived by adding the equations in (1). Let $s=S/N$ , $e=E/N$ , $i=I/N$ and $r=R/N$ denote the fractions of the classes S, E, I and R in the population, respectively. Then for $s^{\prime }=\left(N\cdot S^{\prime }-S\cdot N^{\prime }\right)/N^2$ , $e^{\prime }=\left(N\cdot E^{\prime }-E\cdot N^{\prime }\right)/N^2$ , $i^{\prime }=\left(N\cdot I^{\prime }-I\cdot N^{\prime }\right)/N^2$ and $r^{\prime }=\left(N\cdot R^{\prime }-R\cdot N^{\prime }\right)/N^2$ respectively, the following assignments are appropriate;

$\begin{array}{l}{s}^{\prime }=\frac{N\cdot \left(\alpha N+ϵR-\left(\mu +\beta I/N+\lambda \right)S\right)-S\cdot \left(\left(\alpha -\mu \right)N-k_{E}E-k_{I}I\right)}{N^2}\\ =\alpha +ϵr-\left(\alpha +\lambda \right)s-\left(\beta -k_I\right)si+k_{E}se\\ e^{\prime }=\frac{N\cdot \left(\beta SI/N-\left({\mu }_E+\sigma +{\delta }_E\right)E\right)-E\cdot \left(\left(\alpha -\mu \right)N-k_{E}E-k_{I}I\right)}{N^2}\\ =\beta si-\left(\alpha +\sigma +k_E+{\delta }_E\right)e+k_E{e}^2+k_{I}ie\end{array}$

$\begin{array}{l}{i}^{\prime }=\frac{N\cdot \left(\sigma E-\left(\tau +\gamma \right)I\right)-I\cdot \left(\left(\alpha -\mu \right)N-k_{E}E-k_{I}I\right)}{N^2}\\ =\sigma e-\left(\alpha +k_I\right)i+k_{E}ei+k_I{i}^2\\ r^{\prime }=\frac{N\cdot \left(\gamma I-\left(\mu +ϵ\right)R+\lambda S+{\delta }_{E}E\right)-R\cdot \left(\left(\alpha -\mu \right)N-k_{E}E-k_{I}I\right)}{N^2}\\ =\gamma i-\left(\alpha +ϵ\right)r+\lambda s+{\delta }_{E}e+k_{E}er+k_{I}ir\end{array}$ (2)

From the transformation above, the following equations are derived

$\begin{array}{l}{s}^{\prime }=\alpha +ϵr-\left(\alpha +\lambda \right)s-\left(\beta -k_I\right)si+k_{E}se\\ e^{\prime }=\beta si-\left(\alpha +\sigma +k_E+{\delta }_E\right)e+k_E{e}^2+k_{I}ie\\ i^{\prime }=\sigma e-\left(\alpha +k_I\right)i+k_{E}ei+k_I{i}^2\\ r^{\prime }=\gamma i-\left(\alpha +ϵ\right)r+\lambda s+{\delta }_{E}e+k_{E}er+k_{I}ir\end{array}$ (3)

subject to the restriction that $s+e+i+r=1$ . Let us hereon refer to (3) as the full model. It is also worthy to note that the total population $N\left(t\right)$ does not appear in Equation (2). We can attribute that to the homogeneity in the system (1). Also, since r appears in the first equation of the Equation (3), we substitute $r=1-s-e-i$ into the first equation to get:

$\begin{array}{l}{s}^{\prime }=\alpha +ϵ-\left(\alpha +\lambda +ϵ\right)s-\left(\beta -k_I\right)si+k_{E}se-ϵe-ϵi\\ e^{\prime }=\beta si-\left(\alpha +\sigma +k_E+{\delta }_E\right)e+k_E{e}^2+k_{I}ie\\ i^{\prime }=\sigma e-\left(\alpha +k_I\right)i+k_{E}ei+k_I{i}^2\end{array}$ (4)

so that we can focus our attention on the subsystem (4).

$\begin{array}{l}{s}^{\prime }=\left(\alpha +ϵ\right)-\left(\alpha +\lambda +ϵ\right)s-\left(\beta -k_I\right)si+k_{E}se-ϵe-ϵi\\ e^{\prime }=\beta si-\left(\alpha +\sigma +k_E+{\delta }_E\right)e+k_E{e}^2+k_{I}ie\\ i^{\prime }=\sigma e-\left(\alpha +k_I\right)i+k_{E}ei+k_I{i}^2\end{array}$ (5)

The transformation has revealed some dynamics and interactions embedded in our model which is not intended to be measured in this study. Let us rewrite system (5) by making these substitutions; $h_1=\alpha +ϵ$ ; $h_2=\alpha +\lambda +ϵ$ , $h_3=\alpha +\sigma +k_E+{\delta }_E$ ; and $h_4=\alpha +k_I$ to obtain the system of equations:

$\begin{array}{l}{s}^{\prime }=h_1-h_{2}s-\left(\beta -k_I\right)si\\ e^{\prime }=\beta si-h_{3}e\\ i^{\prime }=\sigma e-h_{4}i\end{array}$ (6)

For biological considerations we study the system (6) in the region

$\Omega =\left\{\left(s,e,i\right)\in {\Re }_{+}^3|0\le s+e+i\le 1\right\}$

It can be shown that $\Omega$ is positively invariant. In the latter part of this study, we show that the transformed model (3) is similar to the new model (6) by numerical simulation. [19] also did this transformation in their study and used the full model, but for the purpose of this study, with the established equality of the models (3) and model (6) by numerical computation we turn our attention to the latter model.

4. Equilibria & Stability

It can be seen that Equation (6) has a disease-free equilibrium is $DFE=\left(\frac{h_1}{h_2},0,0\right)$ . The Jacobian matrix evaluated at the DFE is given thus,

$DFE=\left(\begin{array}{ccc}-h_2& 0& \frac{h_1\left(k_I-\beta \right)}{h_2}\\ 0& -h_3& \frac{\beta h_1}{h_2}\\ 0& \sigma & -h_4\end{array}\right)$

The characteristic polynomial of the DFE system $F\left(\lambda \right)$ is given as,

$\begin{array}{c}F\left(\lambda \right)={\lambda }^3+\left(h_2+h_3+h_4\right){\lambda }^2+\left(h_2{h}_3+h_2{h}_4+h_3{h}_4-\frac{\beta \sigma h_1}{h_2}\right)\lambda \\ +h_2{h}_3{h}_4-\beta \sigma h_1\end{array}$ (7)

From the characteristic polynomial above,

$b_1=h_2+h_3+h_4$ (8)

$b_2=h_2{h}_3+h_2{h}_4+h_3{h}_4-\frac{\beta \sigma h_1}{h_2}$ (9)

$b_3=h_2{h}_3{h}_4-\beta \sigma h_1$ (10)

By the Routh-Hurwitz’s criterion, for the characteristic equation to have negative eigenvalues, then $b_1>0$ , $b_3>0$ and $b_1{b}_2>b_3$ . The Routh-Hurwitz’s criterion is a tool used to establish the negativity or otherwise of the roots of a characteristic equation. Let us first introduce our threshold parameter $R_0$ .

4.1. Basic Reproductive Number, R₀

Following the theory made explicit in [4] on $R_0$ above, we use the features of the model under study to carve out an appropriate measure for the $R_0$ . Let $\mathcal{F}$ , $\mathcal{V}$ F, V and $f\left(x\mathrm{,}y\right)$ be defined as in [4]. From the system (6),

$\mathcal{F}=\left[\begin{array}{c}\left(\beta -k_I\right)si\\ 0\end{array}\right]\text{and}\mathcal{V}=\left[\begin{array}{c}{h}_{3}e\\ -\sigma e+h_{4}i\end{array}\right]$ (11)

From (11), we proceed to find F and V as defined in [4] above.

$F=\left[\begin{array}{cc}0& \left(\beta -k_I\right)s\\ 0& 0\end{array}\right]\text{and}V=\left[\begin{array}{cc}{h}_3& 0\\ -\sigma & h_4\end{array}\right]$ (12)

It can be easily shown that

$V^{-1}=\frac{1}{h_3{h}_4}\left[\begin{array}{cc}{h}_4& 0\\ \sigma & h_3\end{array}\right]$ (13)

Consequently,

$\begin{array}{l}{R}_0=\rho \left(F{V}^{-1}\right)=\frac{1}{h_3{h}_4}\left[\begin{array}{cc}0& \left(\beta -k_I\right)s\\ 0& 0\end{array}\right]\times \left[\begin{array}{cc}{h}_3& 0\\ -\sigma & h_4\end{array}\right]\\ =\frac{\left(\beta -k_I\right)\sigma s_{DFE}}{h_3{h}_4},\text{but}{s}_{DFE}=\frac{h_1}{h_2}\\ R_0=\frac{\left(\beta -k_I\right)h_1\sigma }{h_2{h}_3{h}_4}\end{array}$ (14)

Another construction of the measure $R_0$ mostly considered by researchers long before the introduction of the new generation matrix (NGM) method was introduced is

$R_0=\frac{DFE}{EE}$

where $EE$ is the endemic equilibrium.

$\begin{array}{l}{R}_0=\frac{DFE}{EE}=\frac{\beta \sigma h_1}{h_2{h}_3{h}_4}\hfill \end{array}$ (15)

We remark that, the NGM’s measure of $R_0$ (as shown in 14) is equivalent to this measure (15) if and only if $k_I=0$ . But for the purpose of this study, we choose (14) as our measure for $R_0$ .

Theorem 1. Denote $R_0=\frac{\left(\beta -k_I\right)h_1\sigma }{h_2{h}_3{h}_4}$ . When $R_0\le 1$ , then the system (6) has only a $DFE=\left(\frac{h_1}{h_2},0,0\right)$ on the set $\Omega$ . When $R_0>1$ then besides the DFE, there exist a unique endemic equilibrium EE where the disease is persistent.

4.2. Routh-Hurwitz’s Criterion

Theorem 2 (Routh-Hurwitz’s Criterion) Consider the characteristic equation

$\left|\lambda I-A\right|={\lambda }^n+b_1{\lambda }^{\left(n-1\right)}+\cdots +b_{\left(n-1\right)}\lambda +b_n=0$ (16)

determining the n eigenvalues $\lambda$ of a real $n\times n$ square matrix A, where I is the identity matrix. Then the eigenvalues $\lambda$ all have negative real parts if ${\Delta }_1>0$ , ${\Delta }_2>0$ , $\cdots$ , ${\Delta }_n>0$ , where

${\Delta }_k=\left|\begin{array}{cccccccc}{b}_1& 1& 0& 0& 0& 0& \cdots & 0\\ b_3& b_2& b_1& 1& 0& 0& \cdots & 0\\ b_5& b_4& b_3& b_2& b_1& 1& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ b_{2k-1}& b_{2k-2}& b_{2k-3}& b_{2k-4}& b_{2k-5}& b_{2k-6}& \cdots & b_k\end{array}\right|$

The above theorem (2) is reproduced from [20] and used in the article [21].

Theorem 3. The disease-free equilibrium DFE is locally stable as $R_0\le 1$ and unstable as $R_0>1$ .

Proof. It is easy to show from (7) that $b_1>0$ and $b_2>0$ whenever $R_0\le 1$ . Now for $b_1{b}_2>b_3$ ,

$\begin{array}{l}\Rightarrow b_1{b}_2>b_3\\ \Rightarrow \left(h_2+h_3+h_4\right)\left(h_2{h}_3+h_3{h}_4-\frac{\beta \sigma h_1}{h_2}\right)>h_2{h}_3{h}_4-\beta \sigma h_1\\ \Rightarrow h_2^2{h}_3+h_2{h}_3{h}_4-\beta \sigma h_1+h_2{h}_3^2+h_3^2{h}_4-\frac{\beta \sigma h_1{h}_3}{h_2}+h_2{h}_3{h}_4+h_3{h}_4^2-\frac{\beta \sigma h_1{h}_4}{h_2}\\ >h_2{h}_3{h}_4-\beta \sigma h_1\end{array}$

And since $R_0\le 1$ ,

$h_2^2{h}_3+h_2{h}_3^2+h_3^2{h}_4+h_3{h}_4^2+h_2{h}_3{h}_4>\frac{\beta \sigma h_1}{h_2}\left(h_3+h_4\right)$

therefore $b_1{b}_2>b_3$ . The DF system has negative eigenvalues, hence the DF system is locally stable when $R_0\le 1$ . Again, when $R_0>1$ , then the condition 2 of the Routh-Hurwitz’s criterion is violated; thus $b_2$ may become less than 0 and therefore the system becomes unstable.

4.3. Global Stability of the DFE

The matrix theoretic method is used to prove the global stability of the DFE.

A Matrix-Theoretic Method

The matrix-theoretic method is used to prove the sharp threshold statement in theorem (5). It is a systematic method, and it is presented to guide the construction of a Lyapunov function. Taking the same path as [4] [22] [23], let us set

$f\left(x,y\right):=\left(F-V\right)x-\mathcal{F}\left(x,y\right)+\mathcal{V}\left(x,y\right)$ (17)

Then the equation for the disease compartment can be written as

$x^{\prime }=\left(F-V\right)x-f\left(x,y\right)$ (18)

Let ${\psi }^{\text{T}}\le 0$ be the left eigenvector of the non-negative matrix $V^{-1}F$ corresponding to the eigenvalue $\rho \left(V^{-1}F\right)=\rho \left(F{V}^{-1}\right)=R_0$ . The following result provides a general method to construct a Lyapunov function for [1.1]. [17] [24] used this Lyapunov function involving the Perron eigenvectors to study the global dynamics of several specific disease models while [4] used it to consider a general case for infectious diseases. In this paper, this method reproduces to establish the global stability of our system.

Theorem 4 Let F, V and $f\left(x\mathrm{,}y\right)$ be defined as in (1.2) and (2.1) in [4] respectively. If $f\left(x\mathrm{,}y\right)\ge 0$ , in the $\Omega \subset {\Re }_{+}^{n+m}$ , $F\ge 0$ , $V^{-1}\ge 0$ and $R_0\le 1$ , then the function $\mathcal{D}={\psi }^{\text{T}}{V}^{-1}x$ is a Lyapunov function for the system (1.1) (again in [4]) on $\Omega$ .

Proof. The proof as followed in [4] gives

$\begin{array}{l}{\mathcal{D}}^{\prime }={\mathcal{D}}^{\prime }|=\psi V^{-1}{x}^{\prime }=\psi V^{-1}\left(F-V\right)x-\psi V^{-1}f\left(x,y\right)\\ =\left(R_0-1\right){\psi }^{\text{T}}x-{\psi }^{\text{T}}{V}^{-1}\left(x,y\right)\end{array}$ (19)

Since ${\psi }^{\text{T}}\ge 0$ , $V^{-1}\ge 0$ , and $f\left(x\mathrm{,}y\right)\ge 0$ in the region $\Omega$ , the last term is non-positive. If $R_0\le 1$ , then ${\mathcal{D}}^1\le 0$ in $\Omega$ and thus $\mathcal{D}$ is a Lyapunov function for the system [1.1] as defined in [4].

[4] has proven that the Lyapunov function used to prove the global stability of the DFE in $\Omega$ can also be extended to establish a uniform persistence and thus establish the existence of an EE in ${\Re }_{+}^{n+m}$ . Find the theorem 2.2 in [4] and the proof thereof.

Theorem 5. The DFE is globally stable.

Proof. Taking the left side expansion of (19), the Lyapunov function $\mathcal{D}$ as defined in theorem (4) therefore is

$\mathcal{D}={\psi }^{\text{T}}{V}^{-1}x=\left(\frac{\sigma }{h_4^2}+\frac{\sigma }{h_3{h}_4}\right)e+\frac{1}{h_3}i$ (20)

${\mathcal{D}}^{\prime }\mathrm{<mo>\; =\; </mo>}\left(R_0-1\right)\left(\frac{\sigma }{h_4}e+i\right)-\left(\frac{2\sigma }{h_3{h}_4}+\frac{1}{h_4}\right)f\left(x\mathrm{,}y\right)$

it is easy to prove that

$f\left(x,y\right)=0$ (21)

${\mathcal{D}}^{\prime }=\left(R_0-1\right)\left(\frac{\sigma }{h_4}e+i\right)$ (22)

It is obvious from this equation that when $R_0\le 1$ then ${\mathcal{D}}^{\prime }<0$ and therefore the DFE is globally asymptotically stable and unstable when $R_0>1$ .

4.4. Existence and Stability of the Endemic Equilibrium

Referring to theorem (1), there exists a unique EE. In this section, we establish the stability of the EE using Routh-Hurwitz’s criterion for the proof of local stability and Lyapunov functions for the global stability.

4.4.1. Local Stability of Endemic Equilibrium

Theorem 6. The EE is locally stable.

Proof. The Jacobian of the EE is given thus

$J_{EE}=\left(\begin{array}{ccc}-\frac{\left(\beta -k_I\right)\left(\beta h_1\sigma -h_2{h}_3{h}_4\right)}{\left(\beta -k_I\right)h_3{h}_4}-h_2& 0& -\frac{h_3{h}_4\left(\beta -k_I\right)}{\beta \sigma }\\ \frac{\beta \left(\beta h_1\sigma -h_2{h}_3{h}_4\right)}{\left(\beta -k_I\right)h_3{h}_4}-h_2& -h_3& \frac{h_3{h}_4}{\sigma }\\ 0& \sigma & -h_4\end{array}\right)$

Then the characteristic polynomial of the EE would be

$\begin{array}{c}F\left(\lambda \right)=-\frac{\beta h_1{h}_3{h}_4{k}_I\sigma }{h_3{h}_4{k}_I-\beta h_3{h}_4}-\frac{\lambda \beta h_1{h}_4{k}_I\sigma }{h_3{h}_4{k}_I-\beta h_3{h}_4}-\frac{\lambda \beta h_1{h}_3{k}_I\sigma }{h_3{h}_4{k}_I-\beta h_3{h}_4}-\frac{{\lambda }^2\beta h_1{k}_I\sigma }{h_3{h}_4{k}_I-\beta h_3{h}_4}\\ +\frac{{\beta }^2{h}_1{h}_3{h}_4\sigma }{h_3{h}_4{k}_I-\beta h_3{h}_4}+\frac{\lambda {\beta }^2{h}_1{h}_4\sigma }{h_3{h}_4{k}_I-\beta h_3{h}_4}+\frac{\lambda {\beta }^2{h}_1{h}_3\sigma }{h_3{h}_4{k}_I-\beta h_3{h}_4}+\frac{{\lambda }^2{\beta }^2{h}_1\sigma }{h_3{h}_4{k}_I-\beta h_3{h}_4}\\ +\frac{h_2{h}_3^2{h}_4^2{k}_I}{h_3{h}_4{k}_I-\beta h_3{h}_4}+\frac{\lambda h_2{h}_3{h}_4^2{k}_I}{h_3{h}_4{k}_I-\beta h_3{h}_4}+\frac{\lambda h_2{h}_3^2{h}_4{k}_I}{h_3{h}_4{k}_I-\beta h_3{h}_4}+\frac{{\lambda }^2{h}_2{h}_3{h}_4{k}_I}{h_3{h}_4{k}_I-\beta h_3{h}_4}\\ -\frac{\beta h_2{h}_3^2{h}_4^2}{h_3{h}_4{k}_I-\beta h_3{h}_4}-\frac{\lambda \beta h_2{h}_3{h}_4^2}{h_3{h}_4{k}_I-\beta h_3{h}_4}-\frac{\lambda \beta h_2{h}_3^2{h}_4}{h_3{h}_4{k}_I-\beta h_3{h}_4}-\frac{{\lambda }^2\beta h_2{h}_3{h}_4}{h_3{h}_4{k}_I-\beta h_3{h}_4}\\ -\lambda h_2{h}_4-{\lambda }^2{h}_4-\lambda h_2{h}_3-{\lambda }^2{h}_3-{\lambda }^2{h}_2-{\lambda }^3\end{array}$ (23)

which is of the form $F\left(\lambda \right)={\lambda }^3+b_1{\lambda }^2+b_2\lambda +b_3$ where $b_1$ , $b_2$ and $b_3$ are simplified from Equation (23) as

$b_1=h_3+h_4+\frac{\beta \sigma h_1}{h_3{h}_4}$

$b_2=\beta \sigma h_1+\frac{h_1}{h_3}\sigma$

$b_3=\beta \sigma h_1-h_2{h}_3{h}_4$

when $F\left(\lambda \right)=0$ . It can be seen that $b_1$ , $b_2$ and $b_3$ are all positive. For $b_1{b}_2>b_3$ ,

$\begin{array}{c}{b}_1{b}_2=\left(h_3+h_4+\frac{\beta \sigma h_1}{h_3{h}_4}\right)\left(\beta \sigma h_1+\frac{h_1}{h_3}\sigma \right)>\beta \sigma h_1-h_2{h}_3{h}_4=b_3\\ =\beta \sigma h_1{h}_3+h_1\sigma +\beta \sigma h_{1}h4+\frac{h_1{h}_4}{h_3}\sigma +\frac{{\left(\beta \sigma h_1\right)}^2}{h_3{h}_4}+\frac{\beta {\left(\sigma h_1\right)}^2}{h_3^2{h}_4}\\ >\beta \sigma h_1-h_2{h}_3{h}_4=b_3\\ =\left(h_3+h_4\right)\beta \sigma h_1+h_1\sigma +\frac{h_1{h}_4}{h_3}\sigma +\frac{{\left(\beta \sigma h_1\right)}^2}{h_3{h}_4}+\frac{\beta {\left(\sigma h_1\right)}^2}{h_3^2{h}_4}\\ >\beta \sigma h_1-h_2{h}_3{h}_4=b_3\end{array}$ (24)

Now $h_1<h_3+h_4$ and by extension $\left(h_3+h_4\right)\beta \sigma h_1>\beta \sigma h_1$ . The rest of the terms in the left hand side of this inequality are non-negative and therefore obviously greater than $-h_2{h}_3{h}_4$ . Therefore the EE is locally stable.

4.4.2. Global Stability Analysis of the Endemic Equilibrium

A general form of Lyapunov functions coined from the first integral form of the Goh-Volterra system which is often used in the literature of mathematical biology is used to prove the global stability of the EE. This function takes the form

$\begin{array}{l}{\mathcal{L}}_i={\displaystyle {\sum }_{i=1}^n{c}_i\left(x_i-x_i^{*}-x_i^{*}\ln \frac{x_i}{x_i^{*}}\right)}\hfill \end{array}$ (25)

where x is the variables and $c_i$ are carefully selected constants. This criterion has been used many times in establishing the stability or otherwise of many disease models and also present in [4, 13, 15].

Theorem 7. The EE is globally stable.

Proof. Let $s>s^{*}$ and $i>i^{*}$ ,

${\mathcal{L}}_1=s-s^{*}-s^{*}\ln \frac{s}{s^{*}}$ (26)

${{\mathcal{L}}^{\prime }}_1=-\left(\frac{s^{*}-s}{s}\right)s^{\prime }$ (27)

$=-\left(\frac{s^{*}-s}{s}\right)\left(h_1-h_{2}s-\left(\beta -k_I\right)si\right),$ (28)

and the equilibrium relation for $h_1=h_2{s}^{*}+\left(b-k_I\right)s^{*}{i}^{*}$

$\le -\left(\frac{s^{\mathrm{<mo>\; *\; </mo>}}-s}{s}\right)\left(h_2\left(s^{\mathrm{<mo>\; *\; </mo>}}-s\right)+\left(\beta -k_I\right)\left(s^{\mathrm{<mo>\; *\; </mo>}}{i}^{\mathrm{<mo>\; *\; </mo>}}-si\right)\right)$ (29)

${{\mathcal{L}}^{\prime }}_1\le 0$ , whenever $s>s^{*}$ and $i>i^{*}$ in the region ${\Re }_{+}$ (30)

In line (26), we define the Goh-Volterra function for the first variable s and differentiate it in the second line (27). In the following line, $s^{\prime }$ is substituted for its relation in (6) and evaluated at the equilibrium relation for $h_1$ in line (29). Consider that $s>s^{*}$ and $i>i^{*}$ , (which is a necessary condition for (26) to hold). Then the inequality in line (29) is justified. Therefore, ${\mathcal{L}}_1$ defined above is a Lyapunov function.

Again, let $e>e^{*}$ and ${\mathcal{L}}_2=e-e^{*}-e^{*}\ln \frac{e}{e^{*}}$

${\mathcal{L}}_2=e-e^{*}-e^{*}\ln \frac{e}{e^{*}}$ (31)

${{\mathcal{L}}^{\prime }}_2=-\left(\frac{e^{\mathrm{<mo>\; *\; </mo>}}-e}{e}\right)e^{\prime }$ (32)

${{\mathcal{L}}^{\prime }}_2=-\left(\frac{e^{\mathrm{<mo>\; *\; </mo>}}-e}{e}\right)\left(\beta si-h_{3}e\right)$ (33)

$=-\left(e^{*}-e\right)\left(\frac{\beta si}{e}-h_3\right)$ (34)

$\le -\left(e^{*}-e\right)\left(\frac{\beta s^{*}{i}^{*}}{e^{*}}-h_3\right),$ (35)

then from the system under study, (6), $\frac{i^{*}}{e^{*}}=\frac{s}{h_4}$ and $s^{*}=\frac{h_3{h}_4}{bs}$

$\le -\left(e^{*}-e\right)\left(\frac{\beta h_3{h}_4\sigma }{\beta h_4\sigma }-h_3\right)=-\left(e^{*}-e\right)\times 0=0$ (36)

${{\mathcal{L}}^{\prime }}_2\le 0$ (37)

The next Lyapunov function for the e variable is given in line (31), differentiated and evaluated in lines (32) and (33) respectively. The equilibrium relation for $\frac{\beta s^{\mathrm{<mo>\; *\; </mo>}}{i}^{\mathrm{<mo>\; *\; </mo>}}}{e^{\mathrm{<mo>\; *\; </mo>}}}$ is introduced in lines (35) and (36). Again, the inequality is obviously justified whenever $e>e^{*}$ . The function defined as ${\mathcal{L}}_2=e-e^{*}-e^{*}\ln \frac{e}{e^{*}}$ is successful Lyapunov function.

Lastly, let us set $i>i^{*}$ and ${\mathcal{L}}_3=i-i^{*}-i^{*}\ln \frac{i}{i^{*}}$ ,

${\mathcal{L}}_3=i-i^{*}-i^{*}\ln \frac{i}{i^{*}}$ (38)

${{\mathcal{L}}^{\prime }}_3=-\left(\frac{i^{*}-i}{i}\right)i^{\prime }$ (39)

${{\mathcal{L}}^{\prime }}_3=-\left(\frac{i^{*}-i}{i}\right)\left(\sigma e-h_{4}i\right)$ (40)

$=-\left(i^{*}-i\right)\left(\sigma \frac{e}{i}-h_4\right)$ (41)

$\le -\left(i^{*}-i\right)\left(\sigma \frac{e^{*}}{i^{*}}-h_4\right),\frac{e^{*}}{i^{*}}=\frac{h_4}{s}$ (42)

$\le -\left(i^{*}-i\right)\left(\sigma \frac{h_4}{\sigma }-h_4\right)=0$ (43)

${{\mathcal{L}}^{\prime }}_3\le 0$ (44)

The Lyapunov function for the last variable understudy is introduced in line (38) and differentiated in the following line (39). The equilibrium relation for $e^{\mathrm{<mo>\; *\; </mo>}}/i^{\mathrm{<mo>\; *\; </mo>}}$ is substituted in lines (42) and (43) thereby justifying the inequality.

Therefore $\mathcal{L}$ defined as

$\mathcal{L}={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(x_i-x_i^{*}-x_i^{*}\ln \frac{x_i}{x_i^{*}}\right)={\mathcal{L}}_1+{\mathcal{L}}_2+{\mathcal{L}}_3$

is a Lyapunov function for the system (6). Arbitrary constants $c_i$ can be chosen from ${\Re }_{+}$ and any linear combination of $\mathcal{L}$ would be a Lyapunov function for the system. Hence the proof.

5. Numerical Simulations

In this section, we simulate our model to demonstrate the theoretic results found by this study. We only talk about susceptibility and recovery over time since they explain the rest of the compartments. Note that from the schematic diagram (figure [1]), individuals in the exposed compartments are recruited from the susceptible and the infectious individuals are also recruited from the exposed compartment. We wish you to see that we have shown by this numerical simulation that the first system (1) which was transformed to give the full system (3) is mathematically the same as the reduced model (6). Parameters used for this simulation are shown in Table 1. With all other parameters held constant, we vary vaccination from 0.2 to 0.8 with a step size of 0.2 to see when the disease would die out of the system.

From Figure 2, we can see that the higher the vaccination cover in the population, the lesser the susceptibility to the measles disease in the population. Some significant fraction of the population remains susceptible to measles at a 20% yearly vaccination cover after ten-year span. This means then that in the wake of an outbreak there would still be some that would be impacted by the disease. At a 40% vaccination cover yearly, the susceptible individuals reduce sharply. The same can be said of a vaccination cover of 60% and 80% vaccination cover. It can be drawn from here that when vaccination cover is aimed at 80%, the number of susceptible individuals approaches zero.

In Figure 3, it can be seen that awareness of the disease is already in the population and control steps for treatment and vaccination do not allow for an outbreak to get worse before it becomes better. It can be seen that no matter the

Table 1. Parameter and variable definition & values.

Figure 2. Plot of the susceptible individuals with time.

Figure 3. Plot of the exposed individuals with time.

vaccination cover in the population, the disease will persist in the population for a while and die out. When the time is extended, it would be seen that the disease will completely die out of the system as the exposed individuals approaches zero. It makes sense because, from the schematic diagram (4.1) of the model, the exposed individuals are assumed to go through the infectious stage or attained a natural cure against the disease and move into the recovery stage. Again, we see that the full model and the reduced model are the same.

In Figure 4 also, with time, measles would no longer be seen in the population. At any level of vaccination cover, the disease will die out of the population. Only that, at a higher vaccination rate, more susceptible individuals do not have to get the disease but will gain a strong immunity against the disease.

In Figure 5, the higher the vaccination rate, the more susceptible individuals gain immunity against the disease so that more people fall into the safe haven of being recovered. But it can be seen at 60% vaccination cover is comparable to an 80% vaccination cover. It would therefore be a good practice if vaccination in a population whose system can be compared to that which is studied in this paper vaccinate 80% of its people, newborns and new immigrants so that in any case of an outbreak of measles, the disease will be eradicated in the shortest possible time.

Figure 4. Plot of the infectious individuals with time.

Figure 5. Plot of the recovered individuals with time.

6. Conclusion

In this paper, the global dynamics of measles has been studied in a varying population system such as could be compared to the thriving metropolis of Sunyani in Ghana. A matrix theoretic method and Goh-Volterra types of Lyapunov function were used to establish the stability of the disease-free and endemic equilibria respectively. The basic reproductive number $R_0$ is defined using the New Generation Matrix method and proved as the threshold parameter. The DFE is globally stable whenever $R_0\le 1$ and unstable otherwise. The EE is globally stable when $R_0>1$ . Also, the transformation of the original system (1) uncovered some unnecessary measures embedded in the model which was removed to give the approximate system (6). The latter system was used to replace the former and towards the end of this paper, we showed by numerical simulation that approximate system (3) is the same as (6). The numerical simulation also showed the theoretic dynamics discussed in this paper. It showed that more people will gain immunity against the measles and by extension, an outbreak of measles would not impact the community if vaccination cover is 80% annually.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Earn, D.J.D. (2008) A Light Introduction to Modelling Recurrent Epidemics. In: Mathematical Epidemiology, Springer, Berlin, 317. https://doi.org/10.1007/978-3-540-78911-6_1
[2]	Ghana Health Service and UNICEF Encourage Mothers to Deliver with Help from Skilled Birth Attendants UNICEF. https://www.unicef.org/infobycountry/ghana_62444.html
[3]	NHS Public Health Functions Agreement 2015-16, Service Specification No. 10 Measles, Mumps and Rubella (mmr) Immunisation Programme. https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file /383184/1516_No10_Measles_Mumps_and_Rubella__MMR__Immunisation_Programme _FINAL.pdf
[4]	Shuai, Z. and van den Driessche, P. (2013) Global Stability of Infectious Disease Models Using Lyapunov Functions. SIAM Journal on Applied Mathematics, 73, 1513-1532. https://doi.org/10.1137/120876642
[5]	Smith, H.L., Wang, L. and Li, M.Y. (2001) Global Dynamics of an SEIR Epidemic Model with Vertical Transmission. SIAM Journal on Applied Mathematics, 62, 58-69. https://doi.org/10.1137/S0036139999359860
[6]	Abdelhadi, A. and Hassan, L. (2013) Optimal Control Strategy for SEIR with Latent Period and a Saturated Incidence Rate. ISRN Applied Mathematics, 2013, Article ID: 706848. https://doi.org/10.1155/2013/706848
[7]	Onyejekwe, O.O. and Kebede, E.Z. (2015) Epidemiological Modeling of Measles Infection with Optimal Control of Vaccination and Supportive Treatment. Applied and Computational Mathematics, 4, 264-274. https://doi.org/10.11648/j.acm.20150404.15
[8]	Henshaw, S. and McCluskey, C.C. (2015) Global Stability of a Vaccination Model with Immigration. Electronic Journal of Differential Equations, 92, 110.
[9]	Wang, W., Xin, J. and Zhang, F. (2010) Persistence of an SEIR Model with Immigration Dependent on the Prevalence of Infection. Discrete Dynamics in Nature and Society, 2010, Article ID: 727168. https://doi.org/10.1155/2010/727168
[10]	Li, J. and Ma, Z. (2002) Qualitative Analyses of SIS Epidemic Model with Vaccination and Varying Total Population Size. Mathematical and Computer Modelling, 35, 1235-1243. https://doi.org/10.1016/S0895-7177(02)00082-1
[11]	Clements, C.J. and Cutts, F.T. (1995) The Epidemiology of Measles: Thirty Years of Vaccination. In: Measles Virus, Springer, Berlin, Heidelberg, 13-33. https://doi.org/10.1007/978-3-642-78621-1_2
[12]	Centers of Disease Control (CDC) (1989) Measles Prevention. MMWR Supplements, 38, 1.
[13]	Safi, M.A. and Garba, S.M. (2012) Global Stability Analysis of SEIR Model with Holling Type II Incidence Function. Computational and Mathematical Methods in Medicine, 2012, Article ID: 826052. https://doi.org/10.1155/2012/826052
[14]	Lahrouz, A., Omari, L. and Kiouach, D. (2011) Global Analysis of a Deterministic and Stochastic Nonlinear SIRS Epidemic Model. Nonlinear Analysis: Modelling and Control, 16, 59-76.
[15]	Melesse, D.Y. and Gumel, A.B. (2010) Global Asymptotic Properties of an SEIRS Model with Multiple Infectious Stages. Journal of Mathematical Analysis and Applications, 366, 202-217. https://doi.org/10.1016/j.jmaa.2009.12.041
[16]	Adewale, S.O., Podder, C.N. and Gumel, A.B. (2009) Mathematical Analysis of a TB Transmission Model with DOTS. Canadian Applied Mathematics Quarterly, 17, 1-36.
[17]	Shuai, Z. and Van den Driessche, P. (2011) Global Dynamics of Cholera Models with Differential Infectivity. Mathematical Biosciences, 234, 118-126. https://doi.org/10.1016/j.mbs.2011.09.003
[18]	Peralta, R., Vargas-De-Len, C. and Miramontes, P. (2015) Global Stability Results in a SVIR Epidemic Model with Immunity Loss Rate Depending on the Vaccine-Age. Abstract and Applied Analysis, 2015, Article ID: 341854. https://doi.org/10.1155/2015/341854
[19]	Li, J. and Cui, N. (2013) Dynamic Analysis of an SEIR Model with Distinct Incidence for Exposed and Infectives. The Scientific World Journal, 2013, Article ID: 871393. https://doi.org/10.1155/2013/871393
[20]	Weisstein, E.W. (2000) Routh-Hurwitz Theorem. From MathWorld—A Wolfram Web Resource. http://mathworld.wolfram.com/Routh-HurwitzTheorem.html
[21]	Li, M.Y., Graef, J.R., Wang, L. and Karsai, J. (1999) Global Dynamics of a SEIR Model with Varying Total Population Size. Mathematical Biosciences, 160, 191-213. https://doi.org/10.1016/S0025-5564(99)00030-9
[22]	Castillo-Chavez, C., Feng, Z. and Huang, W. (2002) On the Computation of RO and Its Role on Global Stability. In: Castillo-Chavez, P.C., Blower, S., Driessche, P., Kirschner, D. and Yakubu, A.-A., Eds., Mathematical Approaches for Emerging and Reemerging Infectious Diseases: An Introduction, Springer, Berlin, 229. https://doi.org/10.1007/978-1-4757-3667-0_13
[23]	Van den Driessche, P. and Watmough, J. (2008) Further Notes on the Basic Reproduction Number. In: Mathematical Epidemiology, Springer, Berlin, 159-178.
[24]	Guo, H., Li, M. and Shuai, Z. (2008) A Graph-Theoretic Approach to the Method of Global Lyapunov Functions. Proceedings of the American Mathematical Society, 136, 2793-2802. https://doi.org/10.1090/S0002-9939-08-09341-6