key: cord-0011309-oxjeyw5i authors: Memon, Zaibunnisa; Qureshi, Sania; Memon, Bisharat Rasool title: Mathematical analysis for a new nonlinear measles epidemiological system using real incidence data from Pakistan date: 2020-04-28 journal: Eur Phys J Plus DOI: 10.1140/epjp/s13360-020-00392-x sha: a58e8f6574084d5b03fc98bada5516acf68c19ef doc_id: 11309 cord_uid: oxjeyw5i 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 [Formula: see text] and [Formula: see text] , respectively. While using Lyapunov functions, the equilibria are found to be globally asymptotically stable under the former conditions on [Formula: see text] . However, backward bifurcation shows coexistence of stable endemic equilibrium with a stable measles-free equilibrium for [Formula: see text] . A strategy for measles control based on herd immunity is presented. The forward sensitivity indices for [Formula: see text] 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. 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. a e-mail: zaib.memon@faculty.muet.edu.pk (corresponding author) 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. 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 . The equations describing the model are: Fig. 1 Flow diagram for measles model specified in (2.1) with force of infection λ = β I . 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, 0 < τ < 1 (τ = 1 means perfect vaccine, while τ = 0 represents a vaccine that offers no protection at all). The initial conditions of the model (2.1) are of the form It can be easily shown that the solution of model (2.1) subject to the initial conditions (2.2) exists and is nonnegative for all t ≥ 0. Further, is the positively invariant region for the model (2.1). 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 = 207862518 × μ ≈ 260479. 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 1 10 10 k=1 4685e − 01 is used to measure goodness of the fit which is further confirmed by reasonably small relative error's value (1.4685e−01). Fitting the proposed measles model to the real statistical data using parameters from Table 2 3 Equilibria and stability Given that the total population N is a constant, it is possible to obtain disease-free equilibrium X 0 by equating to zero the right side of equations in system (2.1) as , 0, 0, 0 . Next, we compute the basic reproductive ratio R 0 using only the two equations corresponding to compartments E and I from system (2.1) using the next-generation matrix method [37] . Table 2 The next-generation matrix is the product of matrices F and V −1 where the matrix F represents the rate of infection transmission in these compartments, and the matrix V describes all other transfers across the compartments. The matrices F, V and V −1 are given as The basic reproductive ratio, defined as the spectral radius of the matrix FV −1 , is obtained as . It is easy to prove that if R 0 > 1, in addition to the disease-free equilibrium point X 0 , system (2.1) also has an endemic equilibrium point , where In this section, we analyze the model (2.1) for backward bifurcation [38, 39] . The phenomenon occurs in models having multiple endemic equilibria for R 0 < 1. Evaluating force of infection λ at the endemic equilibrium yields following quadratic equation: Hence, the following theorem is established: Clearly, a > 0 and c > or < 0 according to R 0 < or > 1, respectively. Note that case (i) of Theorem 3.1 suggests the existence of a unique endemic equilibrium for R 0 > 1. 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 R 0 , denoted as R c 0 , given by . Therefore, backward bifurcation occurs for Proof The Jacobian matrix of the system (2.1) at It is easy to verify that three of the eigenvalues of J (X 0 ) are λ 1 = −(ξ + μ) < 0, λ 2 = λ 3 = −μ < 0 (recall that 0 ≤ τ ≤ 1 and remaining parameters are nonnegative). The remaining two eigenvalues λ 4 and λ 5 can be obtained from the equation which gives a quadratic equation with λ 4 and λ 5 as its roots satisfying the following: The two conditions stated in (3.7) imply that both λ 4 and λ 5 have negative real parts provided R 0 < 1. This proves stability of X 0 . If R 0 > 1, one of λ 4 and λ 5 has positive real part. In this case, X 0 is unstable. Proof The Jacobian matrix of system (2.1) at The characteristic equation for J (X * ) is given by where , (3.9) From (3.8), an eigenvalue of J (X * ) is λ = −μ < 0. Further, it is easy to verify from the four equations in (3.9) that Hence, by Routh-Hurwitz stability criterion, the endemic equilibrium point X * is locally asymptotically stable for R 0 > 1. Theorem 3.4 The disease-free equilibrium X 0 is globally asymptotically stable in if R 0 ≤ 1. Proof Consider the Lyapunov function The derivative of this function is The largest compact invariant set in is the singleton set {X 0 }. Hence, by LaSalle's invariance principle [40] , X 0 is globally asymptotically stable in . Theorem 3.5 The endemic equilibrium X * is globally asymptotically stable in for R 0 > 1. Proof Consider the following function It may be noted that Substituting the equilibrium relations from (3.12) in the above expression for L 2 , we obtain (3.13) Using the inequality of arithmetic and geometric means, one can show that (3.14) From (3.13) and (3.14), it follows that L 2 ≤ 0 for R 0 > 1. Furthermore, the equality holds if and only if S = S * , V = V * , I * I = E E * . Substituting S = S * in the first equation of (2.1), we get V = V * , I = I * , E = E * and thus R = R * . Hence, by LaSalle's invariance principle [40] , (S, V, E, I, R) → (S * , V * , E * , I * , R * ) as t → ∞. Therefore, X * is globally asymptotically stable in for R 0 > 1. 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 X 0 (the disease-free equilibrium). Then, In the absence of vaccination, i.e., when ξ = 0, the basic reproductive number given in (3.2) reduces toR Hence, we can write It must be noted that R 0 ≤R 0 . The equality holds only when ρ = 0 (i.e., ξ = 0) or τ = 0. This implies that the vaccine, even not 100% efficient, will certainly diminish the basic reproductive number of the disease. As R 0 ≤ 1 is a necessary and sufficient condition for the elimination of measles (Theorems 3.2 and 3.4 ), it follows from Eq. (4.1) that is also a necessary and sufficient condition for measles control. Here, ρ c denotes herd immunity. Combining Theorems 3.2 and 3.4, we get the following result: Table 3 illustrates the threshold values for the herd immunity ρ c subject to different values for τ andR 0 . For instance, ifR 0 = 14 and τ = 0.97, 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. As the initial transmission and persistence of a disease are both dependent on the basic reproductive number R 0 , we perform a sensitivity analysis to identify the model parameters that have the greatest and the least impact on R 0 . 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 R 0 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 R 0 with respect to each parameter in the model (2.1) is given in Table 4 . (5.2) 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 Page 13 of 21 378 Table 2 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. 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 R 0 straight toward 0, whereas R 0 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. 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 popu- Table 2 lation 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 R 0 . Thus, the epidemic is said to have died out for R 0 < 1, whereas it persists in the case when R 0 > 1. In addition, the measles system undergoes backward bifurcation wherein a stable endemic equilibrium is found to coexist with a stable measles-free equilibrium for R 0 < 1. 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. 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] [48] [49] [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. Measles virus: cellular receptors, tropism and pathogenesis The immune response in measles: virus control, clearance and protective immunity The clinical significance of measles: a review Measles epidemic from failure to immunize World Health Organization, Eastern Mediterranean Vaccine Action Plan 2016-2020: a framework for implementation of the Global Vaccine Action Plan (No. WHO-EM/EPI/353/E). World Health Organization Progress toward measles elimination-Pakistan Mathematical model for control of measles epidemiology Mathematical analysis of effect of area on the dynamical spread of measles A mathematical investigation of vaccination strategies to prevent a measles epidemic Mathematical model for assessing the impact of vaccination and treatment on measles transmission dynamics Mathematical model for the control of measles Modelling the effects of treatment and quarantine on measles Nonstandard finite difference scheme for control of measles epidemiology A deterministic model of measles with imperfect vaccination and quarantine intervention Comparison of fractional order techniques for measles dynamics Computational modelling and optimal control of measles epidemic in human population A Susceptible-Exposed-Infected-Recovered-Vaccinated (SEIRV) Mathematical Model of Measles in Madagascar (Doctoral dissertation Age-structure and transient dynamics in epidemiological systems Seasonal transmission dynamics of measles in China Characteristics of measles epidemics in China (1951-2004) and implications for elimination: a case study of three key locations Mathematical modeling on the control of measles by vaccination: case study of KISII County Mathematical Model for the study of measles in Cape Coast Metropolis Parameterizing state-space models for infectious disease dynamics by generalized profiling: measles in Ontario Realistic population dynamics in epidemiological models: the impact of population decline on the dynamics of childhood infectious diseases: Measles in Italy as an example Modelling vaccination programmes against measles in Taiwan On the existence, uniqueness, stability of solution and numerical simulations of a mathematical model for measles disease Some mathematical models and applications used in epidemics A fractional measles model having monotonic real statistical data for constant transmission rate of the disease Effects of vaccination on measles dynamics under fractional conformable derivative with Liouville-Caputo operator Real life application of Caputo fractional derivative for measles epidemiological autonomous dynamical system Monotonically decreasing behavior of measles epidemic well captured by Atangana-Baleanu-Caputo fractional operator under real measles data of Pakistan Measles outbreak risk in Pakistan: exploring the potential of combining vaccination coverage and incidence data with novel data-streams to strengthen control World Health Organization Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission Backward bifurcation analysis of epidemiological model with partial immunity A Mathematical Model to Study the 2014-2015 Large-Scale Dengue Epidemics in Kaohsiung and Tainan Cities in Taiwan (China Mapping global sensitivity of cellular network dynamics: sensitivity heat maps and a global summation law A methodology for performing global uncertainty and sensitivity analysis in systems biology Illustration of sampling-based methods for uncertainty and sensitivity analysis Factorial sampling plans for preliminary computational experiments On sensitivity estimation for nonlinear mathematical models Parameter estimation and sensitivity analysis of dysentery diarrhea epidemic model Application of reproducing kernel algorithm for solving Dirichlet timefractional diffusion-Gordon types equations in porous media Application of residual power series method for the solution of time-fractional Schrödinger equations in one-dimensional space Numerical solutions of systems of first-order, two-point BVPs based on the reproducing kernel algorithm This manuscript has associated data in a data repository. [Authors' comment: The data used in this study is taken from Measles monthly bulletin published on World Health Organization website (http://www.emro.who.int/vpi/publications/measles-monthly-bulletin.html). See reference [36] .] Conflict of interest The authors declare that they have no conflict of interest.Funding This research did not receive any specific grant from funding agencies in the public, commercial or not-for-profit sectors.