Mathematical analysis for a new nonlinear measles epidemiological system using real incidence data from Pakistan.

Modeling of infectious diseases is essential to comprehend dynamic behavior for the transmission of an epidemic. This research study consists of a newly proposed mathematical system for transmission dynamics of the measles epidemic. The measles system is based upon mass action principle wherein human population is divided into five mutually disjoint compartments: susceptible S(t)-vaccinated V(t)-exposed E(t)-infectious I(t)-recovered R(t). Using real measles cases reported from January 2019 to October 2019 in Pakistan, the system has been validated. Two unique equilibria called measles-free and endemic (measles-present) are shown to be locally asymptotically stable for basic reproductive number R 0 < 1 and R 0 > 1 , respectively. While using Lyapunov functions, the equilibria are found to be globally asymptotically stable under the former conditions on R 0 . However, backward bifurcation shows coexistence of stable endemic equilibrium with a stable measles-free equilibrium for R 0 < 1 . A strategy for measles control based on herd immunity is presented. The forward sensitivity indices for R 0 are also computed with respect to the estimated and fitted biological parameters. Finally, numerical simulations exhibit dynamical behavior of the measles system under influence of its parameters which further suggest improvement in both the vaccine efficacy and its coverage rate for substantial reduction in the measles epidemic. © Società Italiana di Fisica and Springer-Verlag GmbH Germany, part of Springer Nature 2020.

Introduction

Measles is a highly contagious respiratory disease caused by a virus in the Paramyxoviridae family [1, 2]. Clinical symptoms include high fever, cough, conjunctivitis, rhinitis, Koplik’s spots and maculopapular rash. The incubation period for measles is 10–14 days, and the infected individuals usually recover in three weeks of illness without undergoing any complications. However, people suffering from malnutrition or vitamin A deficiency are prone to diarrhea, pneumonia, ear infection, blindness and inflammation of brain [3].

Despite being vaccine preventable, measles continues to pose a serious concern for global health management. The disease has been a primary cause of morbidity and mortality among young children under five years of age. The world has faced measles epidemic several times. California faced measles epidemic between 1988 and 1990 with over 16,000 cases and more than 70 deaths reported [4]. In 2018, Madagascar was affected by measles outbreak which infected 50,000 people and resulted in about 300 deaths, majority of them being children [5]. According to World Health Organization (WHO), measles caused more than 140,000 deaths in 2018. Although vaccination has resulted in a 73% drop in measles deaths worldwide between 2000 and 2018, measles is still prevalent in the developing countries in Asia and Africa [6]. The majority of measles-related deaths occur in countries with poor health infrastructures and low per capita incomes.

Pakistan is among the highly measles burdened countries in WHO’s Eastern Mediterranean Region [7]. There are recurrent measles outbreaks in the country every 8–10 years. In 2016, there were 2845 confirmed measles cases in Pakistan. This number surged to 6791 in 2017 and 33,007 in 2018. These figures account for about 44%, 20% and 51% of the total number of cases reported in the respective years in the Eastern Mediterranean Region comprising 22 countries. Around 130 children died from the disease in 2017, while the number rose to over 300 in 2018 [8].

Immunization is regarded as one of the most cost-effective and successful public health interventions. The WHO recommends two doses of measles vaccine for all children. The first dose given to infants at nine months provides 85% immunity, while a second dose at the age of twelve months imparts 95% immunity to the disease. A Demographic and Health Survey conducted in Pakistan during 2017–2018 indicated the nationwide coverage of the first and second dose of measles vaccine at 73% and 67%, respectively. The survey illustrated the significant variation in the estimates of vaccine coverage among different provinces and federal areas in the country with Sindh (61%, 60%), Punjab (85%, 82%), Khyber Pakhtunkhwa (63%, 50%), Baluchistan (33%, 34%), Azad Kashmir (83%, 75%), Gilgit Baltistan (66%, 62%) and Federally Administered Tribal Areas (35%, 21%). These figures are well below the WHO recommended coverage of 95% for both doses of the vaccine [9].

The epidemic models help to describe the mechanism of disease spread and evaluate strategies for the disease control. In recent decades, there has been a growing interest in the use of deterministic compartmental models to study the dynamics of measles and finding ways for its control and prevention. For example, in [10], the authors have taken into account the impact of asymptomatic individuals on measles dynamics. In [11], the authors have shown that wider distancing between measles-infected and non-infected people proves effective in controlling the disease spread. Smith et al. [12] and Peter et al. [14] examined the role of vaccination on measles dynamics. Garba et al. [13] designed a deterministic model to assess the effect of vaccination and treatment on measles transmission. The effect of quarantine and treatment on measles spread is studied in [15]. Other significant contributions can be found in [16-19]. There are a number of case studies found in the literature related to mathematical study of measles, using deterministic models, focusing different regions of the world, for example, Madagascar [20], London [21], China [22, 23], Kenya [24], Cape Coast [25], Ontario [26], Italy [27], Taiwan [28], Senegal [29] and Afghanistan [30].

The objective of present study is to find, via mathematical modeling, a public health strategy based on using vaccine for efficient control of measles in Pakistan. In particular, we aim to analyze the effect of vaccine efficacy and its coverage in preventing the disease spread in the country. Our motivation derives from a few studies [31-34] in the literature focused on deterministic modeling of measles disease in Pakistan. Each of these studies is based on a four-compartmental SEIR (S-susceptible, E-exposed, I-infectious and R-recovered) model, and none investigates the role of vaccine efficiency and its coverage rate on the disease control. The model in this study is an extension of SEIR model that includes a separate compartment V for the vaccinated class. The SVEIR model is based on the assumption of continuous vaccination. The findings of present study may assist government and public health authorities in formulating strategic vaccination plans to deal with the immunization gaps and thus prevent measles outbreaks.

This paper is organized as follows. In Sect. 2, model is formulated and estimates are obtained for the model parameters. Model equilibria are obtained in Sect. 3 along with analysis of backward bifurcation, local and global stability. Section 4 discusses herd immunity, while a discussion on sensitivity analysis is carried out in Sect. 5. In Sect. 6, numerical simulations are presented to study the effects of various model parameters on the dynamics of measles infection. Conclusion and future research directions are given in Sects. 7 and 8, respectively.

Model description

We formulate a deterministic mathematical model comprising five ordinary differential equations. The total population is divided into five compartments that denote the sub-populations: susceptible (S), vaccinated (V), asymptomatic or exposed (E), symptomatic or infectious (I) and recovered (R). A flow diagram for the model is given in Fig. 1.

Fig. 1

Flow diagram for measles model specified in (2.1)

The equations describing the model are:with force of infection . In (2.1), denotes the recruitment rate and keeps the total population N a constant, the effective contact rate, the vaccination coverage rate, the vaccine efficiency, the natural mortality rate, the rate of developing clinical symptoms and the recovery rate . In this study, the vaccine is assumed to be imperfect, i.e., it does not provide a 100% prevention of infection. Thus, vaccinated individuals also become infected via contact with symptomatic individuals. Note that, ( means perfect vaccine, while represents a vaccine that offers no protection at all). The initial conditions of the model (2.1) are of the formIt can be easily shown that the solution of model (2.1) subject to the initial conditions (2.2) exists and is nonnegative for all . Further,is the positively invariant region for the model (2.1).

Parameters estimation and curve fitting

One of the most important steps to be taken during model validation is the use of real data (if available) which assists to get values of some unknown biological parameters used in the epidemiological model under study. In this connection, real measles incidence cases as given in Table 1 are used for validation of the proposed measles model and also to obtain best fitted values of some unknown biological parameters that occur in the model. For the model in present research analysis, there are seven parameters among which four are to be fitted, whereas remaining three are estimated such as the natural mortality rate of a Pakistani is 66.5 years (1.253133e−03 per month) according to WHO data (year-2018) and the population of Pakistan in 2018 is 207862518 and in this way, the recruitment rate is estimated to be . Further, it is also known from [35] that the measles vaccine is about 97% effective; therefore, the vaccine efficacy is estimated to be 0.97. In addition to these estimated values, values of other parameters are mentioned in Table 2 where the parameters (contact rate), (recovery rate), (rate of developing clinical symptoms) and (vaccination coverage rate) are obtained through parameter estimation technique under lsqcurvefit routine via MATLAB software. The simulation results obtained for the measles incidence cases by fitting the proposed model (2.1) with the real statistics of the first 10 months of 2019 are shown in Fig. 2 along with the respective residuals as depicted in Fig. 3. Figure 2 presents a reasonably good fit thereby including reality to the predictions obtained from the proposed measles model (2.1). The associated average relative error of the fit using the formula is used to measure goodness of the fit which is further confirmed by reasonably small relative error’s value (1.4685e−01).

Table 1

Real statistics of measles infected cases from January to October 2019 in Pakistan [36]

Jan	Feb	Mar	Apr	May	June	July	Aug	Sep	Oct
237	252	397	399	276	168	70	28	23	19

Table 2

Biological parameters used in the proposed epidemic model of measles

Parameter	Description	Value	Source
Π	Recruitment rate of susceptible humans	260,479	Estimated
μ	Natural mortality rate	1.253133e−03	Estimated
τ	Efficacy of vaccine	0.97	Estimated
β	Measles contact rate	1.60056e−07	Fitted
δ	Recovery rate	9.3408	fitted
α	Rate of developing clinical symptoms	9.2373e−01	Fitted
ξ	Vaccination coverage rate	5.8306e−01	Fitted

Fig. 2

Fitting the proposed measles model to the real statistical data using parameters from Table 2

Fig. 3

Residuals for the proposed model using parameters from Table 2

Equilibria and stability

Disease-free and endemic equilibrium points

Given that the total population N is a constant, it is possible to obtain disease-free equilibrium by equating to zero the right side of equations in system (2.1) asNext, we compute the basic reproductive ratio using only the two equations corresponding to compartments E and I from system (2.1) using the next-generation matrix method [37]. The next-generation matrix is the product of matrices and where the matrix represents the rate of infection transmission in these compartments, and the matrix describes all other transfers across the compartments. The matrices and are given asThe basic reproductive ratio, defined as the spectral radius of the matrix , is obtained as

It is easy to prove that if , in addition to the disease-free equilibrium point , system (2.1) also has an endemic equilibrium point , withwhere

Backward bifurcation analysis

In this section, we analyze the model (2.1) for backward bifurcation [38, 39]. The phenomenon occurs in models having multiple endemic equilibria for . Evaluating force of infection at the endemic equilibrium yields following quadratic equation:withHence, the following theorem is established:

Theorem 3.1

The model (2.1) has

a unique endemic equilibrium state if ,

a unique endemic equilibrium if and or ,

two endemic equilibria if and ,

no endemic equilibrium otherwise.

Clearly, and or according to or , respectively. Note that case (i) of Theorem 3.1 suggests the existence of a unique endemic equilibrium for . Further, case (iii) of the theorem describes the condition for occurrence of backward bifurcation, i.e., coexistence of a disease-free and an endemic equilibrium point. To determine whether the model (2.1) has this phenomenon, the discriminant of Eq. (3.5) is equated to zero and the resulting equation is solved for critical value of , denoted as , given bywhere . Therefore, backward bifurcation occurs for .

Stability of equilibrium points

Local stability

Theorem 3.2

The disease-free equilibrium is locally asymptotically stable if and unstable if .

Proof

The Jacobian matrix of the system (2.1) at isIt is easy to verify that three of the eigenvalues of are (recall that and remaining parameters are nonnegative). The remaining two eigenvalues and can be obtained from the equationwhich gives a quadratic equationwith and as its roots satisfying the following:The two conditions stated in (3.7) imply that both and have negative real parts provided . This proves stability of . If , one of and has positive real part. In this case, is unstable.

Theorem 3.3

The endemic equilibrium is locally asymptotically stable if .

The Jacobian matrix of system (2.1) at isThe characteristic equation for is given bywhereFrom (3.8), an eigenvalue of is . Further, it is easy to verify from the four equations in (3.9) thatHence, by Routh–Hurwitz stability criterion, the endemic equilibrium point is locally asymptotically stable for .

Global stability

Theorem 3.4

The disease-free equilibrium is globally asymptotically stable in if .

Consider the Lyapunov functionThe derivative of this function isClearly, if with only for . The largest compact invariant set in is the singleton set . Hence, by LaSalle’s invariance principle [40], is globally asymptotically stable in .

Theorem 3.5

The endemic equilibrium is globally asymptotically stable in for .

Consider the following functionTaking derivative of this function, we obtainIt may be noted thatSubstituting the equilibrium relations from (3.12) in the above expression for , we obtainUsing the inequality of arithmetic and geometric means, one can show thatFrom (3.13) and (3.14), it follows that for . Furthermore, the equality holds if and only if . Substituting in the first equation of (2.1), we get and thus . Hence, by LaSalle’s invariance principle [40], as . Therefore, is globally asymptotically stable in for .

Herd immunity

Not everyone in a given population needs to be immunized in order to eliminate the disease. The fraction of individuals with immunity in the population required to prevent an epidemic is called herd immunity. Let denote the fraction of population which is vaccinated at (the disease-free equilibrium). Then,In the absence of vaccination, i.e., when , the basic reproductive number given in (3.2) reduces toHence, we can writeIt must be noted that . The equality holds only when (i.e., ) or . This implies that the vaccine, even not 100% efficient, will certainly diminish the basic reproductive number of the disease. As is a necessary and sufficient condition for the elimination of measles (Theorems 3.2 and 3.4 ), it follows from Eq. (4.1) thatis also a necessary and sufficient condition for measles control. Here, denotes herd immunity. Combining Theorems 3.2 and 3.4, we get the following result:

Corollary 4.1

Measles can be eliminated from the population if .

Table 3 illustrates the threshold values for the herd immunity subject to different values for and . For instance, if and , this implies that at least 95.7% of the population must be vaccinated in order to be a measles-free community. This highlights the equally important roles of vaccine efficiency and the vaccinated proportion of population in the disease control. It further suggests that both these parameters must be significantly high for the disease to be eliminated.

Table 3

Percent values of corresponding to various values of and

τ	R¯_0
	12	13	14	15	16	17	18
0.95	96.5	97.2	97.7	98.2	98.7	99.1	99.4
0.97	94.5	95.2	95.7	96.2	96.6	97.0	97.4
0.99	92.6	93.2	93.8	94.3	94.7	95.1	95.4

Sensitivity analysis of

As the initial transmission and persistence of a disease are both dependent on the basic reproductive number , we perform a sensitivity analysis to identify the model parameters that have the greatest and the least impact on . There are various techniques available to carry out the required sensitivity analysis including the sensitivity heat map method [41], scatter plots [42], Latin hypercube sampling–partial rank correlation coefficient [43], the Morris [44] and Sobol’ [45] methods and the normalized forward sensitivity index technique [46]. In the present study, this is done by computing sensitivity indices of known as normalized forward sensitivity indices using the following formula:where denotes the model parameter. It has been observed that the basic reproduction number is highly sensitive to the efficacy of vaccine that requires the most attention in order to stop the spread of measles epidemic. On the other hand, the parameter (rate of developing clinical symptoms) is observed to be least effective parameter toward the basic reproduction number during forward sensitivity analysis. The sensitivity index of with respect to each parameter in the model (2.1) is given in Table 4.

Table 4

Sensitivity indices for

Parameters	Sensitivity indices
Π	+ 1
μ	− 0.936782
τ	− 30.1718
β	+ 1
δ	− 0.999866
α	0.00135476
ξ	− 0.0647072

Numerical simulations and discussion

Various numerical simulations are carried out in this section to observe impacts of different biological parameters on the proposed measles model (2.1). In Fig. 4, the burden of the measles epidemic is observed to have declining behavior if the vaccination coverage rate is improved, whereas this burden is substantially reduced for reasonably small value of the measles contact rate as shown in Fig. 5. However, these measures are not easy to be taken in a developing country like Pakistan. Therefore, an effective step could be the improvement in the efficacy of the vaccine which reduces the burden of the epidemic as shown in Fig. 6, but the government of Pakistan will have to strive for making the vaccine cost-effective in order to be affordable by people living near the poverty line. On the other hand, improvement in the recovery rate , as shown in Fig. 7, does reduce the epidemic’s burden, but this measure is not as effective as taking efforts to reduce contacts among infectious and the susceptible ones. Moreover, the higher the rate of developing clinical symptoms more will be the number of infectious individuals as theoretically predicted. This theoretical observation is further confirmed by the infectious population in Fig. 8, but the infection seems to vanish as time goes by.

Fig. 4

Dynamical behavior of infectious population for increasing values of vaccination coverage rate while using parameters from Table 2

Fig. 5

Dynamical behavior of infectious population for decreasing values of measles contact rate while using parameters from Table 2

Fig. 6

Dynamical behavior of infectious population for increasing values of vaccine efficacy while using parameters from Table 2

Fig. 7

Dynamical behavior of infectious population for increasing values of recovery rate while using parameters from Table 2

Fig. 8

Dynamical behavior of infectious population for increasing values of rate of developing clinical symptoms while using parameters from Table 2

In order to further support the analysis and observation regarding the proposed measles model (2.1), contour plots (Figs. 9, 10, 11 and 12) for the basic reproductive number are obtained as function of some biological parameters wherein the contact rate is kept on the vertical axis in an increasing fashion and the parameters are set on the horizontal axis with their increasing values. The most significant observation from these contour plots is that the burden of the measles epidemic can effectively be reduced by reasonably increasing the values of among which the vaccine’s efficacy and vaccine’s coverage rate play important roles. Finally, the small value for the vaccine’s coverage rate is not as bad as the equally small value of the vaccine’s efficacy as depicted in Fig. 13; however, larger value of is capable enough to bring straight toward 0, whereas is not exactly 0 even for the maximum value of . Thus, the vaccine’s efficacy plays the most significant role to reduce the burden of the measles epidemic.

Fig. 9

Contour plot of the basic reproductive number as a function of vaccine efficacy and measles contact rate , whereas the remaining parameters are taken from Table 2

Fig. 10

Contour plot of the basic reproductive number as a function of vaccine coverage rate and measles contact rate , whereas the remaining parameters are taken from Table 2

Fig. 11

Contour plot of the basic reproductive number as a function of recovery rate and measles contact rate , whereas the remaining parameters are taken from Table 2

Fig. 12

Contour plot of the basic reproductive number as a function of rate of developing clinical symptoms and measles contact rate , whereas the remaining parameters are taken from Table 2

Fig. 13

Behavior of the basic reproductive number for increasing values of a vaccine coverage rate and b efficacy of vaccine

Conclusion

These research findings are about development of a new continuous time-invariant system for the measles epidemic under vaccination approach. In this regard, we divided the population into five classes of susceptible, vaccinated, exposed, infectious and recovered human population. Parameters involved in the system are obtained with assistance of parameter estimation technique under nonlinear least squares fitting strategy which later produced best fitted curve for the infectious class of the system to the curve of real experimental measles cases obtained from WHO from the month January 2019 to October 2019, in Pakistan. Thus, the proposed measles system is validated having reasonably small relative error value of 1.4685e−01.

Two unique steady-state solutions are computed for the measles system which are shown to be locally and globally asymptotically stable via Routh–Hurwitz stability theory and Lyapunov functions, respectively, under different constraints imposed upon . Thus, the epidemic is said to have died out for , whereas it persists in the case when In addition, the measles system undergoes backward bifurcation wherein a stable endemic equilibrium is found to coexist with a stable measles-free equilibrium for . The minimum fraction of population that must be vaccinated to achieve herd immunity is determined which shows that both the vaccine efficiency and the vaccinated fraction must be sufficiently high for elimination of measles. The sensitivity analysis shows that the parameter for vaccine efficacy is the one taken to be care of which is further confirmed in various numerical simulations carried out. In addition, there is need to attain higher vaccine coverage rate which is vital to preventing the disease spread. Hence, vaccine efficacy and its coverage both seem to have substantially positive roles for effective control and elimination of measles burden in the suffering community.

Future directions

In future, we aim to work on fractional-order versions of the measles model proposed in this paper and solve them using the techniques followed in [47-50]. Fractional-order operators including Weyl, Riesz, Riemann–Liouville, Caputo, Caputo–Fabrizio, Atangana–Baleanu, Atangana–Gomez, fractal–fractional and others have capability to capture complex and anomalous behavior of dynamical systems that describe a physical or natural phenomenon. The non-local nature of these operators retains memory of the underlying processes which proves to be fruitful in case of epidemiological models since such models are designed to comprehend transmission dynamics of an epidemic which, in turn, has characteristics of memory. Thus, the future works will be devoted to the use of the above operators from fractional calculus for improving the measles system introduced in the present research study.