Modelling and Stability Analysis of the Dynamics of Measles with Application to Ethiopian Data

1. Introduction

Measles is an acute viral illness caused by the pathogen $\mathit{Morbillivirus}$ whose only reservoir is the human host. Despite an extensive vaccination against measles, it is still endemic in many countries and is the main cause of morbidity and mortality in developing regions [1,2]. Clinical symptoms of measles include: cough, runny nose, red eyes, sore throat, fever and a widespread skin rash. If the virus is transmitted to the susceptible ones, they become exposed and pass the latent period within the first 6 to 9 days of exposure, after that the infectious period follows and are able to transmit the disease for 6 to 7 days [3,4]. An individual once infected with measles becomes lifelong immune after running its entire course.

Mathematical modelling has been broadly used as a tool to understand the mechanisms by which infectious diseases spread, to decide on how to control the spread and to minimize expenses in controlling disease outbreaks [5,6,7,8,9]. The predictive nature of these epidemic models which aids in decision making is one of the key strengths of such models. One of the challenges in epidemic modelling is the analyses of the global stability of the disease-free and endemic equilibria. So far, the most powerful method for global stability analyses of infectious diseases is the direct Lyapunov and the geometric methods [10,11,12,13,14]. The authors [10,11] studied the global stability analyses of $SIS$, $SIR$, and $SIRS$ epidemic models with constant recruitment, disease-induced death and standard incidence rate using Lyapunov functions. However, constructing a Lyapunov function to establish the global stability of an $SEIR$ system with constant recruitment, standard incidence, and disease-induced death is still difficult. Moreover, an epidemic model must be validated to check if it is to be actually used [7,15,16,17], which contains comparing the real data with the model solutions in order to evaluate if the model proposed correspond to reality. In most cases, the first model solutions fail to agree with the observed real data, probably due to wrong choice of the initial parameters and lack of adequate model assumptions. Even though, the first model solutions fail to agree with the observed real data, it can be adjusted through a process known as model calibration. Model calibration is the process whereby parameter values that promote a good agreement between the model solutions and observed real data are estimated.

Quite recently, studying the dynamics of the measles disease in the context of mathematics is an important issue [4,18,19,20,21,22,23]. Some authors such as [20,22,23] have studied the dynamics using constant and continuous controls. They have studied the global stability analyses of the disease-free and endemic equilibria. Other authors, for instance, [19] investigated the global stability analysis of the disease dynamics but failed to predict the disease dynamics in the context of a particular situation. There are also others such as [21] who proposed a mathematical model of measles and predicted the disease situation in Bangladesh using double vaccination doses. However, the majority of previous works failed to assume standard incidence in the rate of disease transmission and include disease-induced death rate in the study of the dynamics of the disease. Moreover, few researchers have addressed the estimation of the parameters of the proposed model using real data. Therefore, it is important to study the global dynamics of the measles model by assuming the standard incidence in the disease transmission and estimating the parameters of the proposed model to predict the future situation of a particular region.

This work aims to propose a deterministic compartmental mathematical model for the dynamics of measles disease and calibrate the model to real data of measles outbreaks in Ethiopia. Particularly, we study the long-term global stability analyses of the disease-free and endemic equilibria using the Lyapunov stability method. We also investigate the bifurcation analysis. Next, the model is extended by assuming a fraction of the individuals recruited due to birth or immigration are vaccinated and then model predictions are done to compare model solutions and the real data using some plausible parameters as initial values. Then, a rigorous process of model calibration is done through the solution of an inverse problem to find the best parameter estimates. Finally, the local sensitivity analysis of the basic reproduction number to the estimated parameter values is done to check the robustness of the system proposed.

The work has seven sections: The measles model is formulated in Section 2. In Section 3, we study the mathematical analysis: The basic reproduction number, existence of equilibria, the local and global stability analysis of the equilibria and following this the bifurcation analysis is done in Section 4. The system parameters are estimated in Section 5. In Section 6, numerical simulation of the system is performed. Finally, we conclude the work in Section 7.

2. Mathematical Model Formulation

To formulate the mathematical model of measles, the population is divided into four mutually exclusive classes with each class representing the health condition of an individual. Susceptible individuals (S), individuals carrying the pathogen transmitting the disease but show no clinical symptoms of measles (E), individuals infected with measles, showing the symptoms and can spread the disease (I). Finally, the recovered from measles (R), members who are not infected. The total human population is denoted by $N=S+E+I+R$.

Some of the main assumptions made in the formulation of the deterministic system (1) are: The population is homogeneously mixed. All susceptible individuals have the same probability of being infected. The incidence is assumed to be standard incidence from the assumption that the rate of contact is constant. An individual infected of measles will either die or recover. There exist life long immunity for the measles disease survivors and hence the model proposed is a $SEIR$. Diseases induced death rate exist for the members of the infected group. All parameters are non-negative.

The class of susceptible individuals is recruited by immigration or birth at per capita rate $\Lambda$ and decreased following infection with measles at a rate $\frac{\beta I}{S+E+I+R}$. These population is also decreased by the natural death at a rate $\mu$ which is assumed to be the same for all classes. Hence, the first equation of system (1) is proposed. The population of humans exposed to measles (E) is generated at the rate $\frac{\beta I}{S+E+I+R}$. It is also decreased by the rate of progression to infected class at a rate of $ϵ$ and the natural death rate. Hence, the second equation of system (1) is proposed. The time elapsed in the exposed class is called the latent (or incubation) period. The population of humans infected with measles (I) is increased at the rate of $ϵ$. These individuals recover naturally at a rate of $\gamma$. It is further diminished by the disease-induced death rate at a rate of d and the natural death rate. Hence, the third equation of system (1) is proposed. The population of measles recovered individuals (R) is generated at a rate of $\gamma$. It is also decreased by the natural death rate. Hence, the fourth equation of system (1) is developed.

The system model dynamics for the measles is given by the following deterministic system of nonlinear differential equations and the parameter descriptions are given in Table 1, respectively.

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =\Lambda -\left(\frac{\beta I}{S+E+I+R}+\mu \right)S,\hfill \\ \hfill \frac{dE}{dt}& =\frac{\beta IS}{S+E+I+R}-(\mu +ϵ)E,\hfill \\ \hfill \frac{dI}{dt}& =ϵE-\left(\mu +d+\gamma \right)I,\hfill \\ \hfill \frac{dR}{dt}& =\gamma I-\mu R,\hfill \\ \hfill S\left(0\right)& >0,E\left(0\right)\ge 0,I\left(0\right)\ge 0,R\left(0\right)\ge 0.\hfill \end{array}$$

(1)

3. Model Analysis

In this section, we examine the basic qualitative features which are useful in the subsequent sections.

3.1. Nonnegativity of Solutions and the Feasible Region

Lemma 1. 

The solutions of system (1) with nonnegative initial data remain positive for all $t>0$.

Proof. 

Let $T=sup\left\{t>0:S>0,E>0,I>0,R>0\right\}\in [0,t]$. The first equation of system (1) can be written as

$$\begin{array}{c}\hfill \frac{dS}{dt}=\Lambda -\left(\lambda +\mu \right)S,\mathrm{where}\lambda =\frac{\beta I}{S+E+I+R}.\end{array}$$

(2)

Equation (2) is an initial value problem (IVP) with variable coefficients. It has a unique positive solution $S\left(t\right)$ given by

$$\begin{array}{c}\hfill S\left(t\right)=S\left(0\right)exp\left\{-A\left(t\right)\right\}+exp\left\{-A\left(t\right)\right\}{\int }_0^t\Lambda exp\left\{A\left(s\right)\right\}ds,\end{array}$$

where the function $A\left(t\right)={\int }_0^t\left(\lambda \left(s\right)+\mu \right)ds$ is a particular antiderivative of the function $\left(\lambda \right(s)+\mu )$.

Since T is the maximum of all the time $t\in [0,t]$, we have that $A\left(T\right)={\int }_0^T(\lambda \left(s\right)+\mu )ds$. So that the solution $S\left(T\right)$ is

$$\begin{array}{c}\hfill S\left(T\right)=S\left(0\right)exp\left\{-A\left(T\right)\right\}+exp\left\{-A\left(T\right)\right\}{\int }_0^T\Lambda exp\left\{A\left(s\right)\right\}ds>0.\end{array}$$

Likewise, it is easy to show that all other state variables of the system model remain positive for all $t>0$. □

Lemma 2. 

Every nonnegative solution is bounded in $L^1$-norm by $max\left\{N\left(0\right),\frac{\Lambda }{\mu }\right\}$.

Proof. 

The total dynamics of the system model obtained by adding all the equations of system (1) is:

$$\begin{array}{cc}\hfill \frac{dN}{dt}& =\Lambda -\mu N-dI.\hfill \end{array}$$

(3)

Consequently, the total population $N\left(t\right)$ may vary in time. From 3, we also have $\frac{dN}{dt}\le \Lambda -\mu N$. It is easy to show that

$$0\le N\left(t\right)\le \left(N\left(0\right)-\frac{\Lambda }{\mu }\right)e^{-\mu t}+\frac{\Lambda }{\mu }.$$

(4)

Taking $t\to \infty$, we observe that $0<N\left(t\right)\to \frac{\Lambda }{\mu }$.

The $L^1$-norm of each non-negative solution is N and satisfies $N^{\prime }\le \Lambda -\mu N$ where $N^{\prime }$ represents the derivative of N. Consider a comparison differential equation $N_1^{\prime }=\Lambda -\mu N_1$ where $N^{\prime }\le N_1^{\prime }$. Consider the solutions to the equation $N_1^{\prime }=\Lambda -\mu N_1$. If $N_1\left(0\right)\le \frac{\Lambda }{\mu }$ then ${lim}_{t\to \infty }{N}_1\left(t\right)=\frac{\Lambda }{\mu }$ and $\frac{\Lambda }{\mu }$ is the upper bound. If $N_1\left(0\right)>\frac{\Lambda }{\mu }$, then the solution will decrease to $\frac{\Lambda }{\mu }$ as $t\to \infty$ and $N_1\left(0\right)=N_{10}$ is the upper bound of $N_1$. Since $N^{\prime }\le N_1^{\prime }$ the claim follows for $N\left(t\right)$[24]. □

The biological feasible region for the system model (1) is given by

$$\Omega =\left\{\left(S,E,I,R\right)\in {\mathbb{R}}_{+}^4|0\le S+E+I+R\le \frac{\Lambda }{\mu }\right\}.$$

(5)

In the next subsection we prove the positively invariant and attracting behavior of 5 with respect to system (1).

3.2. Positively Invariant

Theorem 1. 

The nonnegative orthant ${\mathbb{R}}_{+}^4$ is positively invariant for the system model (1).

Proof. 

To prove this Theorem, we first write the system (1) in the form $Y^{\prime }=\mathbf{M}Y+\mathbf{K}$, where

$$\begin{array}{cc}\hfill \mathbf{M}& =\left[\begin{array}{cccc}-(\lambda +\mu )& 0& 0& 0\\ \lambda & -(\mu +ϵ)& 0& 0\\ 0& ϵ& -(\mu +d+\gamma )& 0\\ 0& 0& \gamma & -\mu \end{array}\right]\mathrm{and}\mathbf{K}=\left[\begin{array}{c}\Lambda \\ 0\\ 0\\ 0\end{array}\right].\hfill \end{array}$$

Clearly the column vector $\mathbf{K}\ge 0$ and the matrix $\mathbf{M}$ has the properties of Metzler matrix(all the off diagonal entries of $\mathbf{M}$ is nonnegative). Following the results in [24], the system (1) is positively invariant in ${\mathbb{R}}_{+}^4$. □

3.3. Equilibria and the Basic Reproduction Number

System (1) has a disease-free equilibrium $E_0=(\frac{\Lambda }{\mu },0,0,0)$. In the study of disease dynamics mathematically, the basic reproduction number $R_0$ is the most important parameter to determine disease transmissibility. It determines the spread or die out of the infection with time [13]. Based on the next generation matrix approach developed by [25], the basic reproduction number of system (1) is given as:

$$R_0=\frac{\beta ϵ}{(\mu +ϵ)(\mu +d+\gamma )}.$$

It could be easily noticed that the basic reproduction number $R_0$ is independent of the fraction $\frac{\Lambda }{\mu }$. The endemic equilibrium is given as $E_1=(S^{*},E^{*},I^{*},R^{*})$ with coordinates:

$$\begin{array}{cc}\hfill S^{*}& =\frac{\Lambda \left(\mu \left(\mu +\gamma +d\right)+ϵ\left(\mu +\gamma \right)\right)}{\mu ϵ(\beta -d)},E^{*}=\frac{\Lambda (\mu +\gamma +d)(R_0-1)}{ϵ(\beta -d)},\hfill \\ \hfill I^{*}& =\frac{\Lambda (R_0-1)}{\beta -d},R^{*}=\frac{\gamma \Lambda (R_0-1)}{\mu (\beta -d)}.\hfill \end{array}$$

(6)

It can be easily observed from Equation (6) that the susceptible ($S^{*}$) is feasible only if $\beta >d$. So in what follows we will assume that the effective contact rate ($\beta$) is always greater than the disease induced death rate (d).

3.4. Local Stability Analysis of Disease-Free Equilibrium

Theorem 2. 

The disease-free equilibrium $E_0$ of the system (1) is locally asymptotically stable when $R_0<1$, and unstable if $R_0>1$.

Proof. . The proof involves the evaluation of the Jacobian matrix of the system (1) at $E_0$, which is given by

$$J\left(E_0\right)=\left[\begin{array}{cccc}-\mu & 0& -\beta & 0\\ 0& -(\mu +ϵ)& 0& 0\\ 0& ϵ& -(\mu +d+\gamma )& 0\\ 0& 0& \gamma & -\mu \end{array}\right].$$

The two negative eigenvalues are ${\lambda }_1={\lambda }_2=-\mu$ and the other eigenvalues are found from the submatrix:

$$A=\left[\begin{array}{cc}-(\mu +ϵ)& 0\\ ϵ& -(\mu +d+\gamma )\end{array}\right].$$

Correspondingly, the other eigenvalues are the solutions of the quadratic equation

$${\lambda }^2-\mathrm{trace}\left(A\right)\lambda +det\left(A\right)=0$$

where

$$\mathrm{trace}\left(A\right)=-(2\mu +ϵ+d+\gamma )anddet\left(A\right)=(\mu +ϵ)(\mu +d+\gamma ).$$

Since $\mathrm{trace}\left(\mathrm{A}\right)<0$ and $det\left(A\right)=(\mu +ϵ)(\mu +d+\gamma )>0$ then the quadratic equation, ${\lambda }^2-\mathrm{trace}\left(\mathrm{A}\right)\lambda +det\left(A\right)=0$ has negative eigenvalues [26]. Consequently, the system (1) is locally asymptotically stable at $E_0$ if $R_0<1$. □

3.5. Global Stability of Disease-Free Equilibrium

Theorem 3. 

The disease-free equilibrium, $E_0$, of system (1) is globally asymptotically stable in $\Omega$ if $R_0\le 1$.

Proof. 

To prove this, we define the Lyapunov function $L:\left\{(S,E,I,R)\in \Omega :S>0\right\}\to \mathbb{R}$ by

$$\begin{array}{cc}\hfill L& =ϵE+(\mu +ϵ)I.\hfill \end{array}$$

Differentiating L with respect to time in the solutions of system (1) we get

$$\begin{array}{cc}\hfill \frac{dL}{dt}& =ϵE^{\prime }+(\mu +ϵ)I^{\prime }\hfill \\ \hfill & =ϵ\left(\frac{\beta IS}{S+E+I+R}-(\mu +ϵ)E\right)+(\mu +ϵ)\left(ϵE-(\mu +d+\gamma )I\right)\hfill \\ \hfill & =ϵ\left(\frac{\beta IS}{S+E+I+R}\right)-(\mu +ϵ)(\mu +d+\gamma )I\hfill \\ \hfill & =\frac{I}{S+E+I+R}\left(ϵ\beta S-(\mu +ϵ)(\mu +d+\gamma )(S+E+I+R)\right)\hfill \\ \hfill & =\frac{I}{S+E+I+R}\left((\mu +ϵ)(\mu +d+\gamma )R_{0}S-(\mu +ϵ)(\mu +d+\gamma )(S+E+I+R)\right)\hfill \\ \hfill & =\frac{I(\mu +ϵ)(\mu +d+\gamma )}{S+E+I+R}\left(R_{0}S-(S+E+I+R)\right)\hfill \\ \hfill & =-\frac{I(\mu +ϵ)(\mu +d+\gamma )}{S+E+I+R}\left(-R_{0}S+(S+E+I+R)\right)\hfill \\ \hfill & =-\frac{I(\mu +ϵ)(\mu +d+\gamma )}{S+E+I+R}\left((1-R_0)S+(E+I+R)\right).\hfill \end{array}$$

Then, if $R_0\le 1$, then $\frac{dL}{dt}\le 0$. Furthermore, $\frac{dL}{dt}=0$ if and only if $I=0$ or $E=I=R=0$ and $R_0=1$. Hence, L is a Lyapunov function on $\Omega$. Thus, $(E,I,R)\to (0,0,0)$ as $t\to \infty$. Using $E=I=R=0$ in the first equation of system (1) we obtain $S\to \frac{\Lambda }{\mu }$ as $t\to \infty$. Therefore, the largest compact invariant set in $\left\{(S,E,I,R)\in \Omega :\frac{dL}{dt}(E,I)=0\right\}$ is the singleton $\left\{E_0\right\}$. It follows from the LaSalle’s-Lyapunov theorem [27] that $E_0$ is globally asymptotically stable in $\Omega$. □

3.6. Local Stability of Endemic Equilibrium

Theorem 4. 

The endemic equilibrium is locally asymptotically stable if $R_0>1$.

Proof. 

Linearizing the system (1) around the endemic equilibrium $E_1$, gives the following Jacobian matrix

$$\begin{array}{cc}\hfill J\left(E_1\right)& =\left(\begin{array}{cccc}-x_1& \frac{-\beta I^{*}{S}^{*}}{{(S^{*}+E^{*}+I^{*}+R^{*})}^2}& \frac{-\beta S^{*}(S^{*}+E^{*}+R^{*})}{{(S^{*}+E^{*}+I^{*}+R^{*})}^2}& \frac{-\beta I^{*}{S}^{*}}{{(S^{*}+E^{*}+I^{*}+R^{*})}^2}\\ \frac{\beta I^{*}(E^{*}+I^{*}+R^{*})}{{(S^{*}+E^{*}+I^{*}+R^{*})}^2}& -x_2& \frac{\beta S^{*}(S^{*}+E^{*}+R^{*})}{{(S^{*}+E^{*}+I^{*}+R^{*})}^2}& \frac{-\beta I^{*}{S}^{*}}{{(S^{*}+E^{*}+I^{*}+R^{*})}^2}\\ 0& ϵ& -x_3& 0\\ 0& 0& \gamma & -x_4\end{array}\right)\hfill \end{array}$$

where $x_1=\frac{\beta I^{*}(E^{*}+I^{*}+R^{*})}{{(S^{*}+E^{*}+I^{*}+R^{*})}^2}+\mu >0$, $x_2=(\mu +ϵ)$, $x_3=(\mu +d+\gamma )$, $x_4=\mu$. The characteristic equation of the matrix $J\left(E_1\right)$ evaluated at $E_1$ is

$$\begin{array}{c}\hfill {\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0,\end{array}$$

where

$$\begin{array}{cc}\hfill & a_1=x_1+x_2+x_3+x_4>0,\hfill \\ \hfill & a_2=x_3{x}_4+(x_1+x_2)(x_3+x_4)+x_1{x}_2>0,\hfill \\ \hfill & a_3=x_3{x}_4(x_1+x_2)+x_1{x}_2(x_3+x_4)>0,\hfill \\ \hfill & a_4=x_1{x}_2{x}_3{x}_4>0.\hfill \end{array}$$

This directly follows from the condition $R_0>1$ for $\beta -d>0$.

Next remains to check if $D_i>0$, $i=1,2,3,4$ holds, where

$D_1=a_1$, $D_2=\left|\begin{array}{cc}{a}_1& a_3\\ a_0& a_2\end{array}\right|$, $D_3=\left|\begin{array}{ccc}{a}_1& a_3& a_5\\ a_0& a_2& a_4\\ 0& a_1& a_3\end{array}\right|$, $D_4=\left|\begin{array}{cccc}{a}_1& a_3& a_5& a_7\\ a_0& a_2& a_4& a_6\\ 0& a_1& a_3& a_5\\ 0& a_0& a_2& a_4\end{array}\right|$.

Using matrix manipulation we get

$$\begin{array}{cc}\hfill D_1& =a_1=x_1+x_2+x_3+x_4>0,\hfill \\ \hfill D_2& =\left|\begin{array}{cc}{a}_1& a_3\\ a_0& a_2\end{array}\right|=\left|\begin{array}{cc}{x}_1+x_2+x_3+x_4& x_3{x}_4(x_1+x_2)+x_1{x}_2(x_3+x_4)\\ 1& x_3{x}_4+(x_1+x_2)(x_3+x_4)+x_1{x}_2\end{array}\right|\hfill \\ \hfill & =x_3{x}_4(x_3+x_4)+{(x_1+x_2)}^2(x_3+x_4)+(x_1+x_2){(x_3+x_4)}^2+x_1{x}_2(x_1+x_2)>0,\hfill \\ \hfill D_3& =\left|\begin{array}{ccc}{a}_1& a_3& 0\\ a_0& a_2& a_4\\ 0& a_1& a_3\end{array}\right|=\left|\begin{array}{ccc}{x}_1+x_2+x_3+x_4& x_3{x}_4(x_1+x_2)+x_1{x}_2(x_3+x_4)& 0\\ 1& x_3{x}_4+(x_1+x_2)(x_3+x_4)+x_1{x}_2& x_1{x}_2{x}_3{x}_4\\ 0& x_1+x_2+x_3+x_4& x_3{x}_4(x_1+x_2)+x_1{x}_2(x_3+x_4)\end{array}\right|\hfill \\ \hfill & =a_1(a_2{a}_3-a_1{a}_4)-a_3^2\hfill \\ \hfill & =(x_1^2+x_2^2)(x_1+x_2)(x_3+x_4)x_3{x}_4+{(x_1+x_2)}^2{x}_1{x}_2{(x_3+x_4)}^2+(x_1+x_2){\left(x_1{x}_2\right)}^2(x_3+x_4)\hfill \\ \hfill & +{\left(x_3{x}_4\right)}^2(x_1+x_2)(x_3+x_4)+{(x_1+x_2)}^2{(x_3+x_4)}^2{x}_3{x}_4+(x_1+x_2)x_1{x}_2{(x_3+x_4)}^3>0,\hfill \\ \hfill D_4& =\left|\begin{array}{cccc}{a}_1& a_3& 0& 0\\ a_0& a_2& a_4& 0\\ 0& a_1& a_3& 0\\ 0& 1& a_2& a_4\end{array}\right|=a_1\left|\begin{array}{ccc}{a}_2& a_4& 0\\ a_1& a_3& 0\\ 1& a_2& a_4\end{array}\right|-a_3\left|\begin{array}{ccc}{a}_0& a_4& 0\\ 0& a_3& 0\\ 0& a_2& a_4\end{array}\right|\hfill \\ \hfill & =a_4\left(a_1(a_2{a}_3-a_1{a}_4)-a_3^2\right)\hfill \\ \hfill & =a_4{D}_3>0.\hfill \end{array}$$

Since $D_1$, $D_2$, $D_3$, $D_4$ are all positive, the Routh-Hurwitz criterion [28] for local stability holds. Therefore, the endemic equilibrium is locally asymptotically stable when it exists. □

3.6.1. Global Stability of Endemic Equilibrium

Theorem 5. 

If $R_0>1$ and $d=0$, the endemic equilibrium $E_1$ is globally asymptotically stable in the interior of Ω.

Proof. 

Define $L:\left\{\left(S,E,I,R\right)\in \Omega :S,E,I,R>0\right\}\to \mathbb{R}$ by

$$\begin{array}{cc}\hfill L(S,E,I,R)& =(S-S^{*})+(E-E^{*})+(I-I^{*})+(R-R^{*})-(S^{*}+E^{*}+I^{*}+R^{*})ln\frac{(S+E+I+R)}{(S^{*}+E^{*}+I^{*}+R^{*})}.\hfill \end{array}$$

This function is defined, continuous and positive definite for all $S,E,I,R>0$. It can be verified that the function $L(S,E,I,R)$ takes the value $L(S,E,I,R)=0$ at the equilibrium $E_1$, and thus, the global minimum of $L(S,E,I,R)$ occurs at the endemic equilibrium $E_1=(S^{*},E^{*},I^{*},R^{*})$. At the endemic equilibrium system (1) satisfies

$$\Lambda =\mu (S^{*}+E^{*}+I^{*}+R^{*}),\frac{\beta S^{*}{I}^{*}}{S^{*}+E^{*}+I^{*}+R^{*}}=(\mu +ϵ)E^{*},ϵE^{*}=(\mu +\gamma )I^{*},\gamma I^{*}=\mu R^{*}.$$

(15)

Computing the time derivative of $L(S,E,I,R)$ along the solutions of system (1), we obtain

$$\begin{array}{cc}\hfill \frac{dL}{dt}& =\frac{d}{dt}(S+E+I+R)-\frac{S^{*}+E^{*}+I^{*}+R^{*}}{S+E+I+R}\frac{d}{dt}(S+E+I+R)\hfill \\ \hfill & =\left(\frac{(S-S^{*})+(E-E^{*})+(I-I^{*})+(R-R^{*})}{S+E+I+R}\right)\frac{d}{dt}(S+E+I+R).\hfill \end{array}$$

(16)

Substituting Equation (3) into Equation (16), we get

$$\begin{array}{c}\hfill \frac{dL}{dt}=\left(\frac{(S-S^{*})+(E-E^{*})+(I-I^{*})+(R-R^{*})}{S+E+I+R}\right)\left(\Lambda -\mu \left(S+E+I+R\right)\right).\end{array}$$

Using $\Lambda =\mu (S^{*}+E^{*}+I^{*}+R^{*})$ from Equation (15), we obtain

$$\begin{array}{cc}\hfill \frac{dL}{dt}& =\left(\frac{(S-S^{*})+(E-E^{*})+(I-I^{*})+(R-R^{*})}{S+E+I+R}\right)\left(\mu (S^{*}+E^{*}+I^{*}+R^{*})-\mu \left(S+E+I+R\right)\right)\hfill \\ \hfill & =-\mu \frac{{\left\{(S-S^{*})+(E-E^{*})+(I-I^{*})+(R-R^{*})\right\}}^2}{S+E+I+R}.\hfill \end{array}$$

(17)

It is clear that $\frac{dL}{dt}\le 0$ always holds. Note that, $\frac{dL}{dt}=0$ if and only if $S=S^{*}$, $E=E^{*}$, $I=I^{*}$, and $R=R^{*}$. Hence, the largest invariant compact set in $\{(S,E,I,R)\in \Omega :\frac{dL}{dt}=0\}$ is the singleton $\left\{E_1\right\}$, where $E_1$ is the endemic equilibrium. By the LaSalle’s invariant principle [29], it can be concluded that $E_1$ is globally asymptotically stable in the interior of $\Omega$. □

4. Bifurcation Analysis

Theorem 6. 

The system (1) exhibits a forward bifurcation at $R_0=1$.

Proof. 

Setting $x_1=S$, $x_2=E$, $x_3=I$, and $x_4=R$ we write system (1) as follows

$$\begin{array}{cc}\hfill \frac{d{x}_1}{dt}& =\Lambda -\left(\frac{\beta x_3}{x_1+x_2+x_3+x_4}+\mu \right)x_1:=f_1,\hfill \\ \hfill \frac{d{x}_2}{dt}& =\frac{\beta x_3{x}_1}{x_1+x_2+x_3+x_4}-(\mu +ϵ)x_2:=f_2,\hfill \\ \hfill \frac{d{x}_3}{dt}& =ϵx_2-\left(\mu +d+\gamma \right)x_3:=f_3,\hfill \\ \hfill \frac{d{x}_4}{dt}& =\gamma x_3-\mu x_4:=f_4.\hfill \end{array}$$

(18)

Fix $R_0=1$ and let

$${\beta }^{*}:=\frac{(\mu +ϵ)(\mu +d+\gamma )}{ϵ}$$

be the bifurcation parameter. Jacobian of the linearized of system (18) around the disease-free equilibrium $E_0$ when $\beta ={\beta }^{*}$ is

$$J\left(E_0\right)=\left(\begin{array}{cccc}-\mu & 0& -{\beta }^{*}& 0\\ 0& -(\mu +ϵ)& {\beta }^{*}& 0\\ 0& ϵ& -(\mu +d+\gamma )& 0\\ 0& 0& \gamma & -\mu \end{array}\right).$$

The eigenvalues of the characteristic polynomial are ${\lambda }_1={\lambda }_2=-\mu$, ${\lambda }_3=-(2\mu +d+\gamma +ϵ)$, and ${\lambda }_4=0$. We can observe that the three eigenvalues are real and negative and one is 0 and simple.

To study the bifurcation analysis, we denote the right eigenvector corresponding to the zero eigenvalue ${\lambda }_4=0$ by $\mathbf{w}$=${(w_1,w_2,w_3,w_4)}^T$. Computing $J\left(E_0\right)\mathbf{w}=\mathbf{0}\mathbf{w}$ gives

$$\begin{array}{c}\left(\begin{array}{cccc}-\mu & 0& -{\beta }^{*}& 0\\ 0& -(\mu +ϵ)& {\beta }^{*}& 0\\ 0& ϵ& -(\mu +d+\gamma )& 0\\ 0& 0& \gamma & -\mu \end{array}\right)\left(\begin{array}{c}{w}_1\\ w_2\\ w_3\\ w_4\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right)\hfill \end{array}.$$

Direct calculation gives the right eigenvector

$$\begin{array}{cc}\hfill \mathbf{w}& ={\left(-\frac{\gamma (\mu +d+\gamma )}{ϵ\mu },\frac{\mu (\mu +d+\gamma )}{ϵ},\mu ,\gamma \right)}^T.\hfill \end{array}$$

Next, we compute the left eigenvector $\mathbf{v}$=$(v_1,v_2,v_3,v_4)$ associated with ${\lambda }_4=0$ by computing $vJ\left(E_0\right)=0$ which is

$$\begin{array}{c}(v_1,v_2,v_3,v_4)\left(\begin{array}{cccc}-\mu & 0& -{\beta }^{*}& 0\\ 0& -(\mu +ϵ)& {\beta }^{*}& 0\\ 0& ϵ& -(\mu +d+\gamma )& 0\\ 0& 0& \gamma & -\mu \end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right)\hfill \end{array}.$$

It follows that

$$\begin{array}{c}\hfill \left\{\begin{array}{c}-\mu v_1=0,\hfill \\ -(\mu +ϵ)v_2+ϵv_3=0,\hfill \\ -{\beta }^{*}{v}_1+{\beta }^{*}{v}_2-(\mu +d+\gamma )v_3+\gamma v_4=0,\hfill \\ -\mu v_4=0.\hfill \end{array}\right.\end{array}$$

From the first and fourth equations, we have that $v_1=v_4=0$. Since $\mathbf{v}.\mathbf{w}$=1 one solution of the left eigenvector is $\mathbf{v}={\left(0,ϵ,\mu +ϵ,0\right)}^T$ by taking the expression of ${\beta }^{*}$ into account.

Based on theoretical results in [30], we have to compute the bifurcation coefficients $\mathbf{a}$ and $\mathbf{b}$ given by

$$\begin{array}{cc}\hfill \mathbf{a}& =\sum _{k,i,j=1}^4{v}_k{w}_i{w}_j\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}(E_0,{\beta }^{*}),\hfill \\ \hfill \mathbf{b}& =\sum _{k,i=1}^4{v}_k{w}_i\frac{{\partial }^2{f}_k}{\partial x_i\partial \beta }(E_0,{\beta }^{*}).\hfill \end{array}$$

Since the first and fourth components of $\mathbf{v}$ are zero, we do not need the derivatives of $f_1$ and $f_4$. From the derivatives of $f_2$ and $f_3$, the only nonzero values are:

$$\begin{array}{c}\hfill \frac{{\partial }^2{f}_2}{\partial x_3^2}(E_0,{\beta }^{*})=-\frac{2\mu (\mu +ϵ)(\mu +d+\gamma )}{\Lambda ϵ}\mathrm{and}\frac{{\partial }^2{f}_2}{\partial x_3\partial \beta }(E_0,{\beta }^{*})=1.\end{array}$$

Which directly follows that,

$$\begin{array}{cc}\hfill \mathbf{a}& =v_2{w}_3^2\frac{{\partial }^2{f}_2}{\partial x_3^2}(E_0,{\beta }^{*})=-\frac{2{\mu }^3(\mu +ϵ)(\mu +d+\gamma )}{\Lambda },\hfill \\ \hfill \mathbf{b}& =v_2{w}_3\frac{{\partial }^2{f}_2}{\partial x_3\partial \beta }(E_0,{\beta }^{*})=\mu ϵ.\hfill \end{array}$$

It is clear that $\mathbf{a}<0$ and $\mathbf{b}>0$ for all positive parameter values. Thus, the system (1) exhibits forward bifurcation at $R_0=1$. □

5. Model Calibration with Vaccination

Recently great progress has been done towards measles elimination through constant vaccination against measles globally. Therefore for any mathematical model regarding measles epidemic to be realistic it must include vaccination. In this section, we additionally assume that a fraction of the individuals recruited due to birth or immigration are vaccinated at a rate of $\theta$ and hence go to the recovered class. Vaccination and natural recovery from measles provide the same level of protection. The proportion of those vaccinated is represented by $\theta (0\le \theta \le 1)$ where $\theta =0$ represents vaccination is not applied and $\theta =1$ indicates all individuals recruited due to birth or immigration are vaccinated. The evolution of the system with the incorporation of vaccination is governed by the following nonlinear ordinary differential equations:

$$\begin{array}{cc}\hfill \frac{dS}{dt}& =\Lambda (1-\theta )-\left(\frac{\beta I}{S+E+I+R}+\mu \right)S,\hfill \\ \hfill \frac{dE}{dt}& =\frac{\beta IS}{S+E+I+R}-(\mu +ϵ)E,\hfill \\ \hfill \frac{dI}{dt}& =ϵE-\left(\mu +d+\gamma \right)I,\hfill \\ \hfill \frac{dR}{dt}& =\theta \Lambda +\gamma I-\mu R.\hfill \end{array}$$

(19)

The vaccination reproduction number with the implementation of vaccination is given as

$$R_{0v}=\frac{\beta ϵ(1-\theta )}{(\mu +ϵ)(\mu +d+\gamma )}=(1-\theta )R_0.$$

(20)

It can be easily observed that $R_{0v}\le R_0$ and the equality holds when $\theta =0$. This shows that the vaccine will minimize the $R_0$. For measles eradication, we have obtained that $R_0\le 1$ (Theorem 3) and it follows from Equation (20) that the critical vaccination for measles eradication is $\theta \ge 1-\frac{1}{R_0}={\theta }_c$.

5.1. Data

Ethiopia is one of the countries with high measles mortality burden [18]. Measles is reported weekly by the Ministry of Health. According to the Ethiopian Weekly Epidemiological Bulletin [31] report, the national surveillance data summary of measles disease in 2015 is presented in Figure 1. The real data points are represented by the red dot histogram.

5.2. Nominal Parameters

In order to simulate the proposed measles model in the context of Ethiopian data it is necessary to set all the parameter values. Some of the parameter values are obtain from national census data bases. Among these are the natural birth and death rates and the total population. The natural mortality rate for Ethiopia is $\mu =1/(64\times 52){\left[weeks\right]}^{-1}$[32]. Some others, such as the parameter rate of progress from exposed to infectious class ($ϵ$) and recovery rate ($\gamma$) are also found from measles infected patients’ clinical studies. Considering the assumption that the latent period for measles ends in $6-9$ days and the infectious period as 5 days [33] we choose $ϵ=7.5/7{\left[weeks\right]}^{-1}$ and $\gamma =5/7{\left[weeks\right]}^{-1}$.

5.3. Forward Problem

The system (19) can be given in vector form as:

$$\frac{d\mathbf{x}\left(t\right)}{dt}=\mathbf{f}(\mathbf{x}\left(t\right),\mathbf{p}),\mathbf{x}\left(t_0\right)={\mathbf{x}}_0,$$

(21)

where $\mathbf{x}\left(t\right)=(S,E,I,R)\in {\mathbb{R}}^4$ represents the state variables and ${\mathbf{x}}_0=(S\left(0\right),E\left(0\right),I\left(0\right),R\left(0\right))\in {\mathbb{R}}^4$ is the initial condition at the initial time $t_0$, the vector $\mathbf{p}=(p_1,p_2,p_3,p_4,p_5,p_6)=(\Lambda ,\beta ,\mu ,ϵ,\gamma ,d)\in {\mathbb{R}}^6$ represents the model parameters and $\mathbf{f}:\mathbf{U}\subset {\mathbb{R}}^4\times {\mathbb{R}}^6\to {\mathbb{R}}^4$ is a nonlinear map representing the evolution law of the system defined on the open set:

$$\begin{array}{c}\hfill \mathbf{U}=\left\{\left(x\left(t\right),\mathbf{p}\right)\in {\mathbb{R}}^4\times {\mathbb{R}}^6\mid x_i>0,fori=1,\cdots ,4and{p}_i>0,fori=1,\cdots ,6\right\}.\end{array}$$

(22)

The solution to the forward problem is obtained by numerical integration using a Runge-Kutta fourth and fifth orders ($\mathbf{o}\mathbf{d}\mathbf{e}\mathbf{45}$) method. To simulate the forward problem we assume the initial set of parameters as ${\mathbf{p}}_0=(500,2.59605,1/(64\times 52),7.5/7,5/7,0.168)$. Initially infected humans are assumed as the number of confirmed measles cases in Ethiopia in 2015 which is $I\left(0\right)=250$ individuals. The exposed and infected groups are also equal, $E\left(0\right)=I\left(0\right)$. Moreover, initially recovered individuals are assumed to be $R\left(0\right)=1800$. The susceptible population is also assumed to be $S\left(0\right)=21000$. Accordingly, the initial state variables at the initial time $t_0$ are given as ${\mathbf{x}}_0=(21000,250,250,1800)$. It is also assumed that the proportion of vaccinated in the year 2015 is $\theta =0.65$.

Following the above initial assumptions, Figure 1 depicts the $SEIR$ model solution compared with the real data reported weekly. Apart from the numerical variation between the model solution and the real data, there exists an agreement in the pattern of the model simulation and the real data. This characteristic pattern agreement and numerical disparity proposes that the model solution may vary from the real data due to the use of improper initial parameter values or improper initial state variables in the numerical simulation. Therefore, it is important that model calibration is done to fit the model solution to the observed real data.

Figure 1. Simulation of the model solution compared to the real data using the initial parameter and state values. Observed real data points are represented by the red dot histogram and the best fitting solution curve for the system model is shown by the blue curve.

Figure 1. Simulation of the model solution compared to the real data using the initial parameter and state values. Observed real data points are represented by the red dot histogram and the best fitting solution curve for the system model is shown by the blue curve.

Model calibration is a numerical simulation for the search of parameter values that make the model simulation fit well with the real data. It is discussed in the following subsection.

5.4. Inverse Problem

In most cases, it is impractical the forward problem to reproduce the real data. This could be due to the error in data collection or data registration and error on model formulation caused by wrong model assumptions. Therefore, model calibration consists of seeking suitable set of parameters that, to a certain degree, makes the model solution as close as to the real data.

For the model calibration, the system contains $\mathit{p}=(p_1,p_2,p_3,p_4,p_5,p_6)$ unknown parameters to be estimated based on real data. For the purpose of comparison between the real data observations and model solutions, a discrete time set with $N=52$ weeks presented in Figure 1 is used. The N time discrete real data points are represented by $\mathbf{Z}=(z_1,z_2,...,z_N)$ and the time discrete model solution $\mathit{y}\left(t_k\right)=(y_1,y_2,\cdots ,y_N)$ is generated for the 4 outputs $y_1,y_2,y_3,y_4$ corresponding to N data points.

$$\begin{array}{cc}\hfill & \mathit{y}\left(t_k\right)={[y_1\left(t_k\right),y_2\left(t_k\right),y_3\left(t_k\right),y_4\left(t_k\right)]}^{T}and\hfill \\ \hfill & \mathit{y}={[y_1\left(t_1\right),\cdots ,y_1\left(t_N\right),\cdots ,y_4\left(t_1\right),\cdots ,y_4\left(t_N\right)]}^T(length4·N).\hfill \end{array}$$

The model error (or the the residual error) ${\mathit{e}}_k$ which is the difference between a column vector model solution $\mathit{y}$ and a column vector real data points $\mathit{z}$ is weighted in a quadratic criterion ${\mathit{J}}_N$:

$$\begin{array}{cc}\hfill & {\mathit{e}}_k=\mathit{z}\left(t_k\right)-\widehat{\mathit{y}}(t_k,\widehat{\mathit{p}}),k=1,\cdots ,N,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & {\mathit{J}}_N=\left|\right|\mathit{z}-{\mathit{y}\left|\right|}^2=\sum _{i=1}^M{\mathit{e}}_k^T\mathit{W}{\mathit{e}}_k.\hfill \end{array}$$

(24)

In formal terms, given a real data observation vector z and a model solution vector $\mathit{y}$, the calibration aims at finding a vector of parameters $\widehat{\mathit{p}}$ such that

$$\widehat{\mathit{p}}=arg\underset{\widehat{\mathit{p}}\ge \mathbf{0}}{min}{\mathit{J}}_N\left(\widehat{\mathit{p}}\right),$$

(25)

for a function Equation (16). Here, $\widehat{\mathit{p}}$ is the optimal estimated parameter and $\widehat{\mathit{y}}$ is the model solution corresponding to the optimal $\widehat{\mathit{p}}$ value. $\mathit{W}$ is a $4·N\times 4·N$ positive definite symmetric weighting matrix to be solved using the weighted least squares algorithm. For the optimal estimated parameter values the functional ${\mathit{J}}_N$ reaches a minimum value [34].

For the implementation of the calibration, we used the least squares method using the Trust-Region-Reflective method algorithm lsqnonlin from the MATLAB’s optimization Toolbox. The parameters are estimated with lower and upper bounds. Global minimum convergence is not guaranteed in this case. However, to find the best optimal parameter estimates with minimum functional value, the algorithm is started with different initial parameter and state values. The estimated parameter values are presented in Table 2 and the estimated basic reproduction number of measles for Ethiopia in the year 2015 is found out to be $R_0=1.3973$.

Figure 2a depicts the proposed system model fitted to real data. We found that there are only very few points that are poorly fitted with the model solution as compared. The fit of a system model to a real data is also measured by its residuals which is defined as the difference between the real data and the value predicted by the model. The residual plots seem to be random proving that the model is well calibrated (Figure 2b).

5.5. Local Sensitivity Analysis

Studying the sensitivity of the basic reproduction number as model parameters change is important in system dynamics. It is used to determine the robustness of system predictions to parameter values. A highly sensitive parameter should be carefully estimated, because a small variation in that parameter could lead to large quantitative changes. Less sensitive parameter, does not require as much effort to estimate, since a small change in that parameter will not result in a big influence in the disease dynamics.

Definition 3. 

The normalized forward local sensitivity index of the basic reproduction number $R_{0v}$ that depends differentiably on a parameter $\omega$ is defined by

$${\Pi }_{\omega }^{R_{0v}}=\frac{\partial R_{0v}}{\partial \omega }\frac{\omega }{|R_{0v}|}.$$

(26)

Notice that ${\Pi }_{\omega }^{R_{0v}}$ has a maximum value of magnitude 1. ${\Pi }_{\omega }^{R_{0v}}=1$ implies an increase (decrease) of $\omega$ by $y\%$ increases (decreases) $R_{0v}$ by $y\%$. On the other hand, ${\Pi }_{\omega }^{R_{0v}}=-1$ indicates an increase (decrease) of $\omega$ by $y\%$ decreases (increases) $R_{0v}$ by $y\%$. Following Equation (26), the normalized forward local sensitivity index of $R_{0v}$ with respect to $\beta ,\mu ,ϵ,\gamma ,d,\theta$ is given by

$$\begin{array}{cc}\hfill {\Pi }_{\beta }^{R_{0v}}& =1,{\Pi }_{\mu }^{R_{0v}}=-\frac{\mu (2\mu +ϵ+\gamma +d)}{(\mu +ϵ)(\mu +\gamma +d)},{\Pi }_{ϵ}^{R_{0v}}=\frac{\mu }{\mu +ϵ},\hfill \\ \hfill {\Pi }_{\gamma }^{R_{0v}}& =-\frac{\gamma }{\mu +\gamma +d},{\Pi }_d^{R_{0v}}=-\frac{d}{\mu +\gamma +d},{\Pi }_{\theta }^{R_{0v}}=-\frac{\theta }{1-\theta }.\hfill \end{array}$$

It can be easily observed from Table 3 that the local sensitivity indices of the parameters $\beta$ and $ϵ$ are positive and the remaining are negative. Positive sensitivity indices have positive result on the increase of the basic reproduction number and negative sensitivity indices have reverse effect on the basic reproduction number. The effective contact rate of measles $\beta$ and the rate of vaccination $\theta$ are the most sensitive parameters and significantly impact the disease dynamics. More importantly, decreasing the effective contact rate of measles $\beta$ and increasing the vaccination rate $\theta$ decrease the basic reproduction number. On the other hand, the natural death rate of humans and the rate of progress from the exposed to the infectious class are less sensitive to the basic reproduction number.

6. Numerical Simulations

6.1. The System without Vaccination

In this subsection, numerical simulations of the forward bifurcation of the system (1) and the mesh and contour plots for the basic reproduction number are discussed.

Forward Bifurcation

It is shown in Equation (6) that the disease dies out if $R_0\le 1$ and is endemic if $R_0>1$ for $\beta >d$. These results further prove that the model exhibits forward bifurcation. The global stability of the disease free and endemic equilibria for $R_0\le 1$ and $R_0>1$, receptively is also depicted graphically in Figure 3. The dots colored in red show that the disease free equilibrium is unstable if the reproduction number is greater than unity. On the other hand, the line segment colored in blue shows that the disease free equilibrium is globally asymptotically stable if $R_0\le 1$ and the endemic equilibrium is globally stable if $R_0>1$. $R_0=1$ is the threshold value.

Mesh and Contour Plots

We draw the mesh and contour plots for the basic reproduction number dependence on the effective transmission rate $\beta$ and the natural recovery rate $\gamma$. The mesh and contour plots show that the basic reproduction number increases significantly for higher values of $\beta$ and lower values of $\gamma$ (Figure 4a,b). Thus, to control $R_0$ policymakers are recommended to work on reducing $\beta$ and increasing $\gamma$. The variation of the basic reproduction number as the parameters effective transmission rate $\beta$ and the natural recovery rate $\gamma$ covary is depicted in Figure 4.

We also draw the mesh and contour plots for the basic reproduction number dependence on the effective transmission rate $\beta$ and the rate of progression to the infectious class $ϵ$. The mesh and contour plots show that the basic reproduction number increases significantly for higher values of $\beta$, but $ϵ$ has no significant value (Figure 5a,b). In this case, to control $R_0$ policymakers are recommended to work on reducing $\beta$. The variation of the basic reproduction number as the parameters effective transmission rate $\beta$ and the rate of progression to the infectious $ϵ$ covary is depicted in Figure 5.

6.2. The System with Vaccination

In this subsection, numerical simulations of the proposed system (19) are done. We also validate the importance of the threshold parameter $R_0$ to show disease elimination and endemic cases. The simulations are done using the built-in MATLAB function ode45.

6.2.1. the Disease Extinction Case

The measles disease extinction (Figure 6). The starting points where the solution trajectories start are colored in magenta and the end point of the solution trajectories is colored in red. The direction of the solution trajectories is shown by a forward arrow. In this case the disease-free equilibrium $E_0=$(1690709,0,0,32123479) is globally stable. All solution values tend to move to the disease-free equilibrium what ever the initial value is. The disease free equilibrium is the plane containing the susceptible and recovered classes. From this graphical result it could be concluded that the disease dies out from the community. The calculated value of the reproduction number is $R_0=0.1996$.

6.2.2. the Disease Endemic Case

The measles disease endemic case (Figure 7). The starting points where the solution trajectories start are colored in magenta and the end point of the solution trajectories is colored in red. The direction of the solution trajectories is shown by a forward arrow. In this case the endemic equilibrium $E^{*}=$(2250705,642,324,5601600) is globally stable. All solution values tend to move to the endemic equilibrium. From this graphical result it could be concluded that the disease is endemic. The calculated value of the reproduction number is $R_0=1.3973$.

6.2.3. the Effects of Vaccination and Transmission Rates

The effects of vaccination and transmission rates are discussed here. This is important to identify which one would be the most valid parameter for controlling disease spreading since both are the most sensitive parameters. Figure 8 shows the time series of the infected when the most sensitive parameters are changed. As the rate of vaccination increases from 0.65 to 0.75 it can be observed that the disease can be eradicated from the community (Figure 8a). Compared to the impact of the vaccination rate, lowering the effective contact rate from 7.10724 to 6.10724 does not have a significant impact on the eradication of the disease eradication (creffig:bettavary). The epidemiological implication of this is that vaccination of the susceptible population is the best means for controlling the disease.

It is also worth to note that the $SEIR$ model assumption is not like the $SIR$ in which it cannot be concluded the disease dies out if the number of infected individuals is minimized. The time series of the exposed and infected individuals presented by Figure 9 can prove this hypothesis. It can be observed that the number of individuals in the exposed class is enormous as compared to infected ones. Even though the infected individuals seemingly to decrease which may cause ultimately the disease to die out, the disease may remain endemic due to the transfer of the exposed individuals to the infected class.

7. Discussions and Conclusions

A deterministic $SEIR$ model to describe the dynamics of the measles outbreak in Ethiopia is analyzed and calibrated using weekly real data. The long-term dynamics of the system are analyzed. It is shown that the system is biologically feasible. It is also established that the $R_0$ is a sharp threshold value and it does completely determine the global stability of the disease-free equilibrium for $R_0<1$. And therefore the system exhibits forward bifurcation. Furthermore, the global stability analysis of the endemic equilibrium is analyzed for $R_0>1$ using a suitable Lyapunov function by assuming the disease-induced death rate $d=0$. Based on our mathematical analysis, to eliminate the measles disease from a community, policymakers may work on reducing the basic reproduction number to less than unity. This result is different from the previous results on measles dynamics [4,20], where reducing the basic reproduction number to less than one is not enough for measles eradication. It is important to note that it is still challenging to propose an appropriate Lyapunov function to study the stability of the endemic equilibrium when $d>0$. Therefore, further study using the Lyapunov function theory is needed to determine the global stability of the endemic equilibrium of an SEIR with standard incidence and disease-induced death rate. However, the better way to deal with so far is to use the geometric approach proposed in the works of [12].

The calibration process is done through the solution of an inverse problem using the least squares method using the Trust-Region-Reflective method algorithm lsqnonlin from MATLAB’s optimization Toolbox. Nominal initial parameter and range values are selected from the related literature studying measles disease dynamics. Estimated parameter values with 95% confidence intervals are presented. Furthermore, unlike the previous work by the authors [21] who studied the global sensitivity analysis of the basic reproduction number, local sensitivity analysis of the basic reproduction number to the estimated parameters is done to investigate the robustness of the model as parameters change. It is found that a decrease in the contact rate and an increase in the vaccination rate would cause the disease to die out.

The $SEIR$ model is not like the $SIR$ in that we cannot say the epidemic is in control if the number of infected individuals drops. Definitely, in the $SEIR$ models, it may come about that temporarily the number of people in the infected class is small though the exposed class is still enormous. Consequently, the epidemic may appear died out but will then be out of control again when the people in the exposed class transfer to the infected class and contaminate other people. Therefore, an effective strategy essentially takes into consideration the time necessary for the exposed class to transfer to infected and will minimize to zero new infected cases throughout the period.

Like many epidemic models, the system studied has drawbacks. We calibrated the proposed model to weekly reported cases of measles to estimate the parameters of the model. Estimation of parameters is a tedious work because a large part of the infectious process is missed in most cases [35]. One of the drawbacks is that we could obtain an optimal parameter values that result in a local minimum. Therefore, to solve this problem the recommendation is to simulate the system for several sets of initial parameter values and take the optimal parameters with the smallest error as the best parameter estimate.

Even though we used a comparatively simple model of measles to demonstrate our methodology, the method that we use can be functional on a more complex models by incorporating additional information into the model. For instance, age structure would be very much suitable when considering potential vaccination strategies for measles [5,36]. Furthermore, parameter uncertainty will always exist during the estimation of the parameters and will always influence model predictions. Consequently, an age structured measles model that incorporates uncertainty analysis in the parameter estimation is essential to be investigated.