Modeling transmission dynamics of measles in Nepal and its control with monitored vaccination program

1.   Introduction

Measles is one of the acute and highly contagious viral diseases found as early as the 7th century ^[1]. It is transmitted either by direct contact with infectious droplets or by airborne spread ^[2,3], mainly among children under five years. Before developing vaccines, measles epidemics used to occur every 2 to 5 years, resulting in an annual 20 to 30 million infections and at least 1 million deaths worldwide ^[4,5]. Measles incidence and related deaths decreased globally during 2000–2016 ^[6,7]. However, the measles cases have begun to increase again after 2017. The number of cases in 2019 reached the highest in the past 23 years. The death due to measles also increased by 22% in 2017 and climbed up to 50% in 2019 ^[6,7].

Measles can be successfully prevented with two doses of vaccines received at the recommended ages of 9 months (first dose) and 15 months (second dose) ^[8]. However, due to various reasons, such as poor health systems, lack of access to vaccination, and fear or skepticism about vaccines, the proper implementation of vaccination programs has been a huge issue. Improper and insufficient vaccines have presumably caused frequent measles outbreaks in developed and developing countries, including Nepal. In 2017, the global coverage of the vaccine's first dose was about 85%, significantly lower than the 95% recommended by the WHO to achieve herd immunity ^[5]. Even among developed countries, such as the United States, France, United Kingdom, Argentina, Italy, Japan, Canada, Germany, Australia, and Chile, the vaccine coverage has not met the WHO-recommended threshold for herd immunity ^[5]. Efforts by the anti-vaccine activists in the world ^[9,10] might have partially contributed to the low vaccine coverage and eventually to occasional outbreaks, such as the one in New York in 2018–2019 ^[11,12].

In Nepal, the monovalent vaccine, known as MCV1, against measles was first introduced in three districts in 1979. It was later expanded to the whole country in 1989 ^[13,14,15]. Measles was one of the major causes of childhood death before 2007, presumably because of the low coverage of MCV1. After the MCV1 coverage was increased from 81% to 88% during 2007-2014 along with Supplementary Immunization Activities (SIAs), the suspected measles incidence declined by 13% ^{[14,15,16,17,18,19]}. However, the measles cases began to increase in Nepal in 2017 (99,247, and 430 measles cases in 2017, 2018, and 2019, respectively) ^[20,21]. In 2019 measles outbreaks occurred in Morang, Dang, Kapilvastu, Kathmandu, and Lalitpur districts ^{[20,21,22,23]}. Notably, a frequent measles outbreak has been reported in districts with low vaccine coverage like Kapilvastu, where 95% of children were not vaccinated in 2016 ^[24,25,26]. Even in vaccinated people, the protection level may be reduced due to improper vaccination timing and incomplete doses. Monitored vaccination programs may help achieve success in avoiding measles epidemics. Such monitored vaccination programs promote the timely completion of the vaccination, thereby increasing the chance of complete immunity gain. Mathematical modeling is a valuable tool for identifying the ideal monitored vaccination programs in the context of Nepal.

Many SIR (susceptible, infected, recovered) based deterministic models, including age-structure and immigration-impact, have already been developed for the transmission dynamics of measles ^{[27,28,29,30,31,32]}. Moreover, some SEIR (susceptible, exposed, infected, recovered) based models ^{[33,34,35,36]} have been developed, including continuous-time linear vaccination-based control strategy, meta-populations, and immunization in pregnant women. These basic models have also been extended to include immunity, vaccination, age-dependent vaccination, time-dependent vaccine efficacy, therapy, quarantine, and treatment ^{[2,37,38,39,40,41,42,43,44]}. However, modeling has not been extensively explored to study monitored vaccination programs, especially in the context of Nepal.

Regarding measles in Nepal, some descriptive, analytic, and retrospective studies ^{[14,45,46,47]} have provided insights into the progress in measles control, the Case Fatality Rate (CFR) of measles, and the genetic type of the Asian measles virus. Many of these studies have incorporated the vaccination but lack a monitored vaccination program. In this study, we develop a novel deterministic model, incorporating monitored vaccination programs with two classes of vaccinated individuals, one under the monitored program and another without the monitored program. The model is validated using two-decade-long measles data from Nepal. Our model analysis establishes the local and global stability of disease-free equilibrium, the existence of endemic equilibrium, and the uniform persistence of the disease. Furthermore, we carry out model simulations to properly evaluate the impact of monitored vaccination programs on the short-term and long-term trend of measles transmission in Nepal.

2.   Mathematical model

We develop a transmission dynamics model of measles for the population that includes all the newborns and children under 15 years. As mentioned earlier, implementing a proper vaccination program is often difficult, particularly for children whose parents have poor health knowledge and have fear or skepticism about vaccines. To improve the effectiveness of vaccination programs, we introduce a monitored vaccination program into our model. The program mainly focuses on asserting the completion of vaccines timely and accurately by children under this program. Because of extra care and regular follow-up, the children under this program are expected to have less susceptibility to infection and a higher rate of achieving immunity than those under regular (un-monitored) vaccination programs. To formulate the model, we divide the total population considered $(N)$ into six mutually exclusive compartments: susceptible $(S)$, un-monitored vaccinated ($U_V$), monitored vaccinated ($M_V$), exposed $(E)$, infectious $(I)$, and immune $(I_M)$.

The schematic diagram showing the flow of individuals from and to the compartments during the measles transmission dynamics is presented in Figure 1. The dynamical system equations representing the model are as follows:

Here, $\Lambda$ represents the recruitment rate of susceptible, i.e., newly born children. The measles infection occurs with the per capita transmission rate of $\beta I/{N}$, transferring susceptible individuals to the exposed class. The parameters $\sigma$ represents the rate of progression of individuals from the exposed class to the infectious class, and $\eta$ is the recovery rate of infectious individuals from the disease. Since only the children up to the age of 15 years are considered in the study, we assumed $\mu$ is the rate at which the children become older than 15 years leaving from the dynamics. We take $\delta$ to represent the disease-induced death rate.

The parameters $\alpha_1$ and $\alpha_2$ represent the rate of un-monitored vaccination and the monitored vaccination, respectively. As discussed in Pantha et al. ^[49], for practical purposes, the values of $\alpha_1$ and $\alpha_2$ for un-monitored and monitored vaccination programs aiming to cover $\zeta_1\%$ and $\zeta_2\%$ of children in $t_1$ and $t_2$ years, respectively, can be estimated using $\alpha_1 = -\ln{(1-\zeta_1/100)}/t_1$ and $\alpha_2 = -\ln{(1-\zeta_2/100)}/t_2$, respectively. The individual in the monitored vaccinated group is expected to be properly monitored to ensure the timely completion of the vaccination. In contrast, individuals in the un-monitored vaccinated group have more likelihood of not completing the vaccination in time, possibly delaying the gain of complete immunity. Therefore, we assumed that the un-monitored and monitored vaccinated children become immune at different rates, $\gamma_1$ and $\gamma_2$, respectively. Our model assumes individuals recovered from natural infection and those completing vaccination have similar immunity. Therefore, we include both of them in the same class, namely the immune class ($I_M$). The vaccinated children may also be infected, but at lesser infectivity rates $(1-\varepsilon_1)\beta$ and $(1-\varepsilon_2)\beta$ for un-monitored and monitored vaccinated, respectively. Since the monitoring service providers are expected to counsel individuals in the monitored vaccinated program for timely completion of vaccine and prevention practices, we also expect $0\leq \varepsilon_1 < \varepsilon_2 \leq 1$.

3.   Parameter estimation, data fitting, and model validation

3.1.   Data source

The publicly available data used in this work is obtained from the official site of the World Health Organization (WHO) ^[48]. The data includes the reported measles cases in Nepal from 2000 to 2019. Since 5% of the reported cases belong to the aged 15 and above ^[48], we deducted 5% of the cases from the data. The Crude Birth Rate (CBR) and Infant Mortality Rate (IMR) of Nepal are used from the "Nepal population growth rate 1950–2020" ^[50].

3.2.   Parameter estimation

Nepal's population under fifteen years was 9,807,000 in 2000 (taken as the base year) and 8,460,000 in 2019 ^[48]. The actual population size in individual $S, \; U_V, \; M_V, \; E, \; I, \; I_M,$ classes is not available. It is recorded that $77\%$ of the population was vaccinated in 2000 (the base year) ^[48]. For our base case simulation, among the unvaccinated 23% of the total population, we assumed that $22\%$ were in the susceptible class ($S(0) = 2,157,540$). We took $30\%$ of the vaccinated population were in the un-monitored vaccinated class ($U_V(0) = 2,206,575$). Since the monitored vaccination program was not present in 2000, we took $M_V(0) = 0$. From the data, the recorded cases were 8927 (after reducing 5% over 15 years from the total recorded 9397), among which we assume $E(0) = 300$ in the exposed class and $I(0) = 340$ in the infectious class. The remaining population $(N(0)-S(0)-U_V(0)-E(0)-I(0))$ is included in the immune class $(I_M(0) = 5,442,245)$.

The recruitment rate $(\Lambda = 612,328)$ is the annual average birth rate, which is calculated by using the Crude Birth Rate (CBR) and Infant Mortality Rate (IMR) from the 2000-2019 data ^[50]. Since only the population below fifteen years is considered in the study, we used $\mu = {1}/{15} = 0.0667$ per year. It is given that the incubation period of measles is 10–14 days on average ^[4,51], and thus the disease progression rate from the exposed class to the infectious class is taken as $\sigma = 1/12 \times 365 \approx 30$ per year. Also, since it takes about 18 days (range between 7 to 23 days) to recover from the disease ^[51], we used $\eta = 1/18 \times 365\approx 20$ per year. As per WHO guidelines, children are vaccinated with the first dose at the age of 9 months and the second dose at 15 months ^[8], giving a six-month interval between the two doses. Since the monitored vaccinated individuals ($M_V$) are expected to complete them in time, we took $\gamma_2 = 0.5$ per year. From the data ^[52,53], we estimated the average disease-induced death rate to be $\delta = 0.01$ per year. The remaining parameters of the model, $\beta, \; \alpha_1, \; \alpha_2, \; \gamma_1, \; \varepsilon_1$, and $\varepsilon_2$, are estimated by fitting the model to the measles case data from Nepal.

3.3.   Data fitting and model validation

From the model, the yearly new infections at time $t$ can be calculated using $h(t) = \sigma E(t)$, which we obtained using the numerical solutions of the system (2.1–2.6). Then we estimated the parameters with the help of the nonlinear regression method ^[54], which minimizes the following sum of the square residuals:

where $h_{k}$ denotes the model predicted yearly new infection, $\bar{h_k}$ denotes yearly new infection data, and $n$ is the number of data points used for the model fitting. For each estimated parameter, we also computed the confidence limits using the standard errors from the sensitivity matrix $(\mathscr{S})$ based on the complex-step derivative technique described in the previous study ^{[55,56,57,58]}. Our model is consistent with the yearly incidence cases observed in Nepal (Figure 2). In addition, we also show that the model prediction of the cumulative cases agrees well with the cumulative data, thereby validating our model to describe the measles epidemic in Nepal.

We note that while making all six parameters free in the data fitting process, we obtained negative lower limits of some confidence intervals. To tackle this issue, we needed to fix some of the parameters as done previously ^[55,56,57]. Since the parameters $\alpha_2$ and $\varepsilon_2$ are the two least sensitive parameters identified from the sensitivity matrix, $\mathscr{S}$, we fixed $\alpha_2 = 0.02$ per year and $\varepsilon_2 = 0.90$ at their best estimate values. We also note that taking the different values of $\alpha_2$ and $\varepsilon_2$ did not significantly affect estimates of other parameters, as expected, because of the least sensitivity. Then, we estimated the remaining only four parameters, $\beta, \; \alpha_1, \; \gamma_1, \; \mbox{and }\varepsilon_1$, from the further data fitting. Here the ratio of the data to the free parameters is 5:1, which lies within the recommended range of 5:1 to 10:1 for a reasonable parameter estimate ^[59]. Furthermore, to analyze the identifiability of the estimated four parameters, we use the sensitivity-based method ^[60], in which we obtained the rank of the matrix $\mathscr{S}^T\mathscr{S}$. For our case with four parameters estimated, the obtained full rank (Rank = 4) of the matrix $\mathscr{S}^T\mathscr{S}$ confirms that these four parameters are practically identifiable for the model and the data used in this study. All the computations were carried out in MATLAB (The Math Works. Inc.) using its various routines, including "ode45" (ODE solver) and "fmincon" (minimizer).

Our estimates show that the transmission rate $(\beta)$, the rate of un-monitored vaccination $(\alpha_1)$, the rate of recovery from un-monitored vaccination $(\gamma_1)$, and the effectiveness of un-monitored vaccination $(\varepsilon_1)$ are 63.0238 (95% CI: 63.0188–63.0288), 0.31 (95% CI: 0.0902–0.5298), 0.1 (95% CI: 0.0858–0.1142), and 0.5082 (95% CI: 0.4926–0.5238), respectively (Table 1). We note that the confidence interval of the parameter $\alpha_1$ appears to be large compared to other parameters for this particular data set. However, we consider this parameter a control parameter and vary widely for the analysis of the vaccination program; thus, one set of confidence intervals does not significantly impact the main results of our study.

4.   Model analysis

4.1.   Basic properties of model

We first establish the biological or epidemiological validation of the model by proving the non-negativity and the boundedness of the solution of the system of equations.

4.1.1.   Positivity of the solutions

Theorem 4.1. If $S(0) > 0, \; U_V(0)\geq 0, \; M_ V(0)\geq 0, \; E(0)\geq 0, \; I(0) \geq 0$, and $I_M(0) \geq 0$, then the set of solution $\{S(t), \; U_V(t), \; M_V(t), \; E(t), \; I(t), \; I_M(t)\}$ of the system (2.1–2.6) is positive for all $t \geq 0$.

Proof. From (2.1), $\dfrac {dS}{dt} > -(\beta I/N +\mu+\alpha_{1}+\alpha_{2})S$, which implies $S(t) > S(0)exp { {\left(-\int_{0}^{t}\left({\beta I(s)/N(s) +\mu+\alpha_{1}+\alpha_{2}}\right)ds\right)}}$. Since $S(0) > 0, \; S(t) > 0, \; \forall t > 0$ confirming the positivity of $S(t)$. Similarly from (2.2–2.6), we have $\dfrac{dU_V}{dt} \geq -\left(\gamma_{1}+\mu+ {(1-\varepsilon_1) \beta I}/{N}\right)U_V$, $\dfrac{dM_V}{dt} \geq -(\mu+\gamma_2+(1-\varepsilon_2)\beta I/N)M_V$, $\dfrac{dE}{dt} \geq -(\mu+\sigma)E$, $\dfrac{dI}{dt} \geq -(\mu+\eta+\delta)I$, and $\dfrac{dI_M}{dt}\geq -\mu I_M$. Then we get $U_V(t)\geq U_V(0) exp { {\left(-\int_{0}^{t}\left((1-\varepsilon_1) \beta I(s)/N(s) + \mu+\gamma_{1}\right)ds\right)}} \geq 0$, $M_V (t)\geq M_V(0) exp { {\left(-\int_{0}^{t}\left((1-\varepsilon_2)\; \beta I(s)/N(s) + \mu+\gamma_{2}\right)ds\right)}}\geq 0$, $E(t)\geq E(0)exp \left({ {-(\mu+\sigma)t}}\right) \geq 0$, $I(t)\geq I(0)exp \left({ {-(\mu+\eta+\delta)t}}\right)\geq 0$, and $I_M(t)\geq I_M(0)exp \left({ {-\mu t}}\right) \geq 0, \forall t > 0$ for $U_V(0)\geq 0, \; M_ V(0)\geq 0, \; E(0)\geq 0, \; I(0) \geq 0, \mbox{ and }\; I_M(0) \geq 0$. Therefore, $U_V(t), \; M_V(t), \; E(t), \; I(t), \; I_M(t)\geq 0, \; \forall t > 0$, showing the positivity of the solution set of the system (2.1–2.6).

4.1.2.   Boundedness and invariant region

Adding all differential equations (2.1–2.6), we get ${dN}/{dt} = \Lambda- \mu N-\delta I \leq \Lambda- \mu N,$ which provides $N(t) \leq N(0)e^{ {-\mu t}} + {\Lambda}/{\mu}(1-e^{ {-\mu t}})$ and $\limsup\limits_{t\rightarrow \infty} N(t)\leq {\Lambda}/{\mu}$, showing that the total population is ultimately bounded by ${\Lambda}/{\mu}$. Hence the solution set bounded by $\Lambda/\mu$ is positively invariant in the feasible region

4.2.   Existence of Equilibria

In this section we discuss the disease-free equilibrium point (DFE) and the endemic equilibrium point of the system (2.1–2.6).

4.2.1.   Disease free equilibrium point and formulation of the vaccinated reproduction number

Setting $E = 0$ and $I = 0$ in the system, we obtain DFE: $E^0 = (S^0, \; U_V^0, \; M_V^0, \; 0, \; 0, \; I_M^0)$, where $S^0 = \dfrac{\Lambda }{\alpha _1+\alpha _2+\mu }$, $U_V^0 = \dfrac{\alpha _1 \Lambda }{\left(\alpha _1+\alpha _2+\mu \right) \left(\gamma _1 +\mu \right)}$, $M_V^0 = \dfrac{\alpha _2 \Lambda }{\left(\alpha _1+\alpha _2+\mu \right) \left(\gamma _2+\mu \right)}$, and $I_M^0 = \dfrac{\Lambda \left(\alpha _1 \gamma _1(\mu +\gamma _2) +\alpha _2 \gamma _2(\mu + \gamma _1)\right)}{\mu \left(\alpha _1+\alpha _2+\mu \right) \left(\gamma _1+\mu \right) \left(\gamma _2+\mu \right)}.$ Note that our model implies non-zero immune class $I_M^0$ at the DFE due to the immunity gained through vaccination.

The vaccinated reproduction number, denoted by $R_v$, is defined as the average number of secondary cases generated by a single infectious case introduced into the mixed population with susceptible and vaccinated status. Here we formulate the vaccinated reproduction number using the Next Generation method ^{[61,62,63,64,65]}.

Following the Next Generation Matrix method ^[62], we divide the system into two groups, infected $\vec{x} = (x_i, \; i = 1, 2) = (E, \; I)$ and uninfected $\vec{y} = (y_j, \; j = 1, 2, 3, 4) = (S, \; U_V, \; M_V, \; I_M)$, as follows:

where $\mathscr{F}_i(\vec{x}, \vec{y})$ is the rate of appearance of new infections in the compartment $i$, and $\mathscr{V}_i(\vec{x}, \vec{y})$ is the difference between the transfer of individuals out of and into the compartment $i$ ($i = 1, 2$). Here we have

It is easy to verify that $\mathscr{F}_i$ and $\mathscr{V}_i$ satisfy the conditions A(1)–A(5) of ^[62]. As the process provided in ^[62], we obtained the Jacobians of $\mathscr{F} \mbox{ and } \mathscr{V}$ at the disease free equilibrium point $(D \mathscr{F}(E^0) = F, \; D \mathscr{V}(E^0) = V)$:

Clearly, $F$ is non-negative and $V$ is non singular M-matrix. The next generation matrix for our model is given by:

Then the vaccinated reproduction number $R_v$ is given by:

Using the expression of DFE obtained above, we get

Furthermore, as the system (4.1) satisfies all the conditions of [62,Theorem 2], it follows the local stability of the disease-free equilibrium for $R_V < 1$. Moreover, for the verification purpose, we also prove the local stability of DFE by showing the negative real part of eigenvalues in Subsection 4.3.1 below.

4.2.2.   Endemic equilibrium point

Let

Then solving the full system (2.1–2.6) equated to zero, we get

Substituting (4.3) into (4.2) and after some simplification, we obtain the following equation in terms of $\lambda^*$:

where,

Note that $\lambda^* = 0$ corresponds to the disease-free equilibrium point, and the endemic equilibrium points are given by the solutions of $A_3 \lambda^{*3}+A_2 \lambda^{*2}+A_1 \lambda^{*} + A_0 = 0$. Here, $A_3 > 0$ and $A_0 < 0$ for $R_v > 1$, so the equation (4.4) has at least one positive root for $R_v > 1$. The positiveness of $\lambda^*$ implies the positiveness of $I^*$ and $E^*$. Thus, we conclude that if $R_v > 1$, the system (2.1–2.6) has at least one endemic equilibrium point given by (4.3).

4.3.   The global dynamics analysis

In this section, we first establish the local and global stability of the disease-free equilibrium when $R_v < 1$ and then prove the persistence of the disease when $R_v > 1$.

4.3.1.   Local stability analysis of the disease free equilibrium point

Theorem 4.2. The disease-free equilibrium point of the system (2.1–2.6) is locally asymptotically stable if $R_v < 1$ and unstable if $R_v > 1$.

Proof. Jacobian of the system (2.1–2.6) at the disease free equilibrium point is $J = \left(\begin{array}{cc} A_{3\times 3} & B_{3\times 3} \\ C_{3\times 3} & D_{3\times 3} \\ \end{array} \right)$ where,

The eigenvalues of the matrix $J$ are

Since all the eigenvalues are negative if $R_v < 1$, $E^0$ is locally asymptotically stable if $R_v < 1$.

4.4.   Global stability of disease free equilibrium point

In the following theorem, we prove that $R_v < 1$ asserts the global stability of the DFE.

Theorem 4.3. If $R_v < 1$, the disease-free equilibrium point $E^0$ of the system (2.1–2.5) is globally asymptotically stable in $\mathbb R^6$.

Proof. The Jacobian corresponding to (2.4) and (2.5) at the disease-free equilibrium point $E^0$ is

The spectral bound of $J_0$ is defined as: $s(J_0) = max\{Re(\lambda): \lambda\; \mbox{is an eigenvalue of } J_0\}$. Using the theorem 4.2, we conclude:

1. $R_v = 1 \mbox{ if and only if } s(J_0) = 0$.

2. $R_v < 1 \mbox{ if and only if } s(J_0) < 0$.

3. $R_v > 1 \mbox{ if and only if } s(J_0) > 0$.

If $R_v < 1$ then $s(J_0) < 0$. Also, we can obtain a sufficiently small positive number $\gamma$ such that $s(J_{\gamma}) < 0$, where

From the subsections 4.1 and 4.2, $N(t)\rightarrow \Lambda/\mu$ and $S(t)\rightarrow S^0$ as $t\rightarrow \infty$ for all $\gamma > 0$. Therefore, there exists $t_1 > 0$ such that $\forall t \geq t_1$, we have $S\leq (S^0 +\gamma), \; U_v \leq (U_v^0 +\gamma), \; M_v \leq (M_v^0+\gamma), \; N\geq (\Lambda/\mu -\gamma)$. From (2.4) and (2.5), it follows that

Clearly, the system (4.6) has the Jacobian $J_{\gamma}$, which is irreducible with non-negative off-diagonal elements. Then $s(J_{\gamma})$ is simple and associated with strongly positive eigenvector $\tilde {v}, \; \forall \; t\geq t_1$ ^[66]. For any solution $\psi(t)$ of the system (4.6) with non-negative initial value $\psi(0)$, there is a sufficiently large positive number $\zeta > 0$ such that $(E(t_1), I(t_1))\leq \zeta \tilde {v}$. It is easy to see that ${V}(t) = \zeta e^{ {s(J_\gamma)(t-t_1)}} \tilde{v}, \forall\; t\geq t_1$ is a solution of ($4.6$) with $V(t_1) = \zeta \tilde v$. Then by the comparison principle [66,Theorem B.1], it follows that $(E(t), I(t))\leq \zeta \tilde {v}, \; \; \forall t\geq t_1$. Since $s(J_{\gamma}) < 0$, we get

which implies that (2.1) is asymptotic to the following equation:

Hence, we get $\lim\limits_{t\rightarrow \infty} S(t) = \dfrac{\Lambda}{\left(\mu +\alpha_{1}+\alpha_{2}\right)}$. Solving the system of differential equations (2.2 & 2.3) together with the help of (2.1), we get

where,

Thus, we get

In the limiting case, we can further verify that $\lim\limits_{t\rightarrow \infty} I_M(t) = I_M^0$. Hence if $R_v < 1$, the disease-free equilibrium point $E^0$ of the system (2.1–2.5) is globally asymptotically stable.

4.4.1.   Uniform persistence

The disease is endemic if the system is uniformly persistent. The system (2.1–2.6) is said to be uniformly persistent if there exists $\zeta > 0$ such that $(S(t), \; U_V(t), \; M_V(t), \; E(t), \; I(t), \; I_M(t))\in \Gamma$ satisfying

Here it is enough to consider the reduced system (2.1-2.5). We first define the sets

We have

which implies that $\partial \Gamma^o$ is relatively closed in $\Gamma = (S, U_V, M_V, E, I)\in \mathbb R^5.$

To establish the disease persistence for $R_v > 1$, we now prove the following theorem.

Theorem 4.4. If $R_v > 1$, then the system (2.1–2.5) is uniformly persistent with respect to $(\Gamma^0, \partial \Gamma^0)$ in the sense that there is a positive constant $\zeta > 0$ such that every solution $(S(t), \; U_V(t), \; M_V(t), \; E(t), \; I(t))$ of (2.1–2.5) with $(S(0), \; U_V(0), \; M_V(0), \; E(0), \; I(0)) \in \Gamma^0$ satisfies$\lim\limits_{t\rightarrow \infty}\mathit{\mbox{Inf}} E\geq \zeta, \; \; \lim\limits_{t\rightarrow \infty}\mathit{\mbox{Inf}} I \geq \zeta$.

Proof. Let $\psi(t)p$ be the solution function of system (2.1–2.5) with initial value $p$. We can show the solution $(\psi(t)p), \; t > 0$ is uniformly persistent with respect to $(\Gamma^o, \partial\Gamma^o)$ ^[67]. For any $(S(0), \; U_V(0), \; M_V(0), \; E(0), \; I(0))\in \Gamma^0$, (2.1–2.3) provide

$\mbox{ where }\; a = \alpha_1 +\alpha_2 +\mu +{\beta I}/{N}, \; b = \Lambda, \; a^* = \gamma_1 +\mu+ {(1-\varepsilon_1)\beta I}/{N}$, $b^* = \alpha_1 S, \; a^o = \gamma_2 +\mu+{(1-\varepsilon_2)\beta I}/{N}$, $\mbox{ and }\; b^o = \alpha_2 S$. Here $\Lambda > 0$ implies $S(t) > 0$. This follows $U_V(t) > 0, \; M_V(t) > 0, \; \forall \; t > 0$. Thus the non-diseased variables $S(t), \; U_V(t), \mbox{and}\; M_V(t)$ are positive.

Defining the two sets,

we claim that $\omega(p) = \{E^0\}, \; \forall p\in M_{\partial}$. If $p(t)\in M_{\partial}$ then $\psi(t)p\in {\partial}\Gamma_0$ which implies $I = 0$ and $E = 0.$ From (4.7, 4.8) and (2.1–2.5), as $t \rightarrow \infty$ we obtain $S(t), U_V(t)$, and $M_V(t)$ approaching to

which is the equilibrium point $E^0$. Hence $\omega(p) = \{E^0 \}, \; \; \forall p\in M_{\partial}$.

Using the theorem 4.3 (Condition 3, $R_v > 1$ for $s(J_0) > 0$), we can get a sufficiently small $\rho > 0$ such that the perturbation $J_\rho$ on $J_0$ satisfies $s(J_{\rho}) > 0$, (see ^[67], [68,Section II.5.8]) for $R_v > 1$, where

Now we claim that disease-free equilibrium point $E^0$ is uniform weak repeller with any solution $\psi(t)$. For this we need to show

On contrary, suppose there exists a $p_0 \in \Gamma^0$ such that $\lim\limits_{t\rightarrow \infty}\mbox{ Sup }\left\|\psi(t)p_0 -E^0\right \| < \rho$. From the subsections 4.1 and 4.2, $N(t)\rightarrow \Lambda/\mu$ and $S(t)\rightarrow S^0$ as $t\rightarrow \infty$ implying that for all $\rho > 0$ there exists $t_2 > 0$ such that $\forall t \geq t_2,$ $S\geq (S^0 -\rho), \; U_v \geq (U_v^0 -\rho), \; M_v \geq (M_v^0-\rho), \; N\leq (\Lambda/\mu +\rho)$. Then (2.4) and (2.5) follows that,

We consider the following auxiliary equations:

Here, $J_{\rho}$ is the Jacobian of the system (4.12), is irreducible with non-negative off-diagonal elements, then $s(J_{\rho})$ is simple and associated with strongly positive eigenvector $\tilde {v}, \forall\; t\geq t_2$, $(E(t_2), I(t_2)) > 0$. Thus there is a positive number $\xi > 0$ such that $(E(t_2), I(t_2))\geq \xi \tilde {v}$ holds. It is easy to see that ${V} = \xi e^{ {s(J_{\rho})(t-t_2)}} \tilde{v}, \; \; \forall\; t\geq t_2$ is a solution of (4.12) with ${V}(t_2) = \xi\tilde{v}$. Hence by comparison principle [66,Theorem B.1], we get

Also, for $R_v > 1$, we have $s(J_{\rho}) > 0$, implying $\lim\limits_{t\rightarrow \infty} E(t) = \infty, \; \lim\limits_{t\rightarrow \infty} I(t) = \infty,$ which is a contradiction. This proves that the solution repels from $E^0$. It follows that forward orbit of any solutions in $M_{\partial}$ converges to $E^0$ is isolated in $\mathbb R^5$.

Now, we define a stable set of $E^0$:

Clearly, $W^s(E^0)\cap \Gamma^0 = \emptyset$. It follows that there is no cycle in $M_{\partial}$ from $E^0$ to $E^0$. Applying [69,Theorem 1.3.1], we conclude that the system (2.1–2.5) is uniformly persistent, i.e., there exists $\zeta > 0$ satisfying (4.9).

5.   Numerical results

5.1.   Basic dynamics of measles in Nepal

We first present the properties of basic dynamics of measles in Nepal, particularly vaccinated reproduction numbers and long-term dynamics.

5.1.1.   Computation of the reproduction numbers

Using the estimated parameters (Table 1), we obtain the reproduction number of measles in the presence of vaccination in Nepal to be $R_v = 1.0098$. While $R_v > 1$ is consistent with the currently ongoing endemic of measles, the actual magnitude we observed is significantly lower than the previously estimated reproduction number between $5$ and $18$ in other places ^[70]. The low value of $R_v$ is expected as its estimate is based on the parameters influenced by the vaccination program. We now compute the time-dependent effective reproduction number, $R_t$, which describes the time-varying average number of secondary cases. The value of $R_t$ allows us to track whether the epidemic at time $t$ is in an increasing $(R_t > 1)$ or decreasing $(R_t < 1)$ trend. For our model, the effective reproduction number is given by

Using the estimated parameters (Table 1), we obtained the pattern of $R_t$ as shown in Figure 3. As predicted by our model, the value of the effective reproduction number remains less than unity, indicating the epidemic is in decreasing trend until 2030, after which $R_t$ increases and remains greater than unity until 2095, showing the increasing trend of the disease in the period 2030–2095 (Figure 3).

5.1.2.   Longterm dynamics

In this section, we present our model prediction for the long-term dynamics of measles transmission in Nepal. If the current trend continues, there will still be 178 cases in 2023 (Figure 4), indicating an obstacle to the measles elimination program set by the government of Nepal. Our model predicts that the cases remain at a significantly low level until 2033 but persist at a low level without being eradicated. After 2033, the resurgence of the outbreak will begin and reach the peak value (7670) in 2097. The dynamic is relatively slow, reaching a steady-state only after about 2100. The resurgence of measles in Nepal predicted by our model supports the worldwide trend of the epidemic, which shows an eventual resurgence in many places such as the UK, the US (31 states including New York), and the Philippines. We also predict the contribution of susceptible (unvaccinated), un-monitored vaccinated, and monitored vaccinated to the resurgence of the measles cases. As per our model prediction, the major contribution to the resurgence is from the susceptible and the un-monitored vaccinated groups. Note that the un-monitored vaccinated group can be significantly high in Nepal, as shown by the data that while 92% are vaccinated with MCV1, only about 76% are vaccinated with MCV2 in 2019 ^[48]. Therefore, proper implementation of the monitored vaccination program may be needed to avoid the resurgence of the disease.

5.2.   Sensitivity analysis

5.2.1.   Sensitivity of parameters to ${R_v}$

We first observe the local sensitivity of $R_v$ to each of the parameters. For this, we obtain the sensitivity index $S_x$, given by

where $x$ is the parameter of which the sensitivity is to be obtained. Based on $S_x$, we found that $R_v$ is highly sensitive to $\beta$. The parameter $\gamma_1$ also affects $R_v$ more than the other parameters, and the $\varepsilon_2$ affects less while the effect of $\sigma$ and $\alpha_2$ are negligible (Figure 5).

We also extend the analysis to the global sensitivity by using Latin Hypercube Sampling (LHS) ^[71], taking 1000 sample points from the global parameter space. We compute the partial rank correlation coefficients to identify the most influential parameters. We observed that in the global parameter space, the parameters $\varepsilon_1, \; \alpha_1$, and $\beta$ are the most strongly effective to $R_v$, followed by $\alpha_2, \varepsilon_2, \gamma_1$, and $\sigma$, while $\gamma_2$ is less effective (Figure 5).

5.2.2.   Sensitivity of parameters to the dynamics

In this section, we use Latin Hypercube Sampling ^[71] from the global parameter space to identify the sensitivity level of the peak value of the infected class and that of the time at the peak of the epidemic. The computed partial rank correlation coefficient corresponding to each parameter is presented in Figure 6. Our analysis shows that the peak value of the infected class is highly correlated to $\beta$ (positive correlation). The peak value is moderately affected by $\alpha_1, \alpha_2$, $\gamma_1$, and $\sigma$ and is weakly correlated to the parameters $\gamma_2$, $\varepsilon_2$, and $\varepsilon_1$ (Figure 6). Similarly, the peak time of the epidemic is mostly affected by the parameters $\beta$ and $\varepsilon_1$, while it is less influenced by the parameters $\alpha_1$ and $\gamma_1$. The parameters $\sigma, \gamma_2, \alpha_2$, and $\varepsilon_2$ have the least effect on the peak time.

5.3.   Effects of monitored and un-monitored vaccination

We use five parameters, the un-monitored vaccination rate $(\alpha_1)$, the monitored vaccination rate $(\alpha_2)$, the immunization rate of un-monitored vaccination $(\gamma_1)$, the effectiveness of un-monitored vaccination $(\varepsilon_1)$, and the effectiveness of monitored vaccination $(\varepsilon_2)$ to evaluate the effects of the vaccination programs.

5.3.1.   Effects on the measles eradication

Note that $R_v < 1$ can be associated with a condition required to eradicate the disease (Section 4.3). We now use our $R_v$ formulation to identify vaccine-related parameters that can reduce $R_v$ below one (Figure 7). We found that at the current level of $\alpha_1$, either the monitored vaccination rate $\alpha_2$ needs to rise by 250% $(\alpha_2 \geq.07$) or the effectiveness of un-monitored vaccination $\varepsilon_1$ needs to rise by 4% $(\varepsilon_1 \geq 0.52)$ to make $R_v < 1$. On the other hand, with the current effectiveness of monitored vaccination ($\varepsilon_2$), the effectiveness of un-monitored vaccination ($\varepsilon_1$) needs to be greater than 0.52 to make $R_v < 1$ (Figure 7).

If the un-monitored vaccination rate ($\alpha_1$) is increased by 20% ($\alpha_1 \approx 0.42 \; \mbox{per year}$), the activity of the monitored vaccination program can be somewhat relaxed at the current level of $\alpha_2 \approx 0.02 \; \mbox{per year}$ for $R_v < 1$. On the other hand, if $\alpha_1$ is decreased by 26% ($\alpha_1 \approx 0.23 \; \mbox{ per year}$), the activity of the monitored vaccination program needs to be raised, making the level of $\alpha_2$ exceed more than $0.1 \mbox{ per year}$ to achieve $R_v < 1$.

5.3.2.   Effects on the dynamics

In the absence of monitored vaccination, the model predicts that the cases will rise and reach the peak value of 85,750 in 2038 for $\alpha_2 = 0$. With the reduction of un-monitored vaccination by 50% or $\alpha_1 = 0.15$, our model predicts that the cases will rise and reach the peak value of 647,300 in 2028 (Figure 8). If the immunity rate of un-monitored vaccination is decreased by 50%, i.e., $\gamma_1 = 0.05$, the cases rise, leading to the peak value of about 452,800 in 2030. Similarly, if the effectiveness of un-monitored vaccination is decreased by 20% $(\varepsilon_1 = 0.4)$, the cases may rise and reach the peak value of 268,200 in 2027.

5.3.3.   Effects on the steady state

We explore how the different levels of vaccination-related parameters, $\alpha_1, \; \alpha_2, \; \gamma_1, \gamma_2, \varepsilon_1, \mbox{ and }\; \varepsilon_2$, (Figure 9) affect the level of steady-state of new infections (Figure 9). As expected, an increase in the vaccination rate decreases the steady-state level of new infections. Based on our model prediction, the current level of vaccination results in 7,312 new cases at a steady state. To sufficiently reduce the steady-state level of new infections, the rate of un-monitored vaccination should be increased by 6.5% $(\alpha_1 = 0.33)$, or the rate of monitored vaccination should be increased by 50% $(\alpha_2 = 0.03)$. Similarly, the immunizing rate of un-monitored vaccine needs to be increased by 100% $(\gamma_1 = 0.2)$, or its effectiveness needs to be increased by 20% $(\varepsilon_1 = 0.6)$ to bring the steady-state level of new infections to a sufficiently low level.

Furthermore, the reduction of 6.5% of the un-monitored vaccination rate $(\alpha_1)$ increases the steady-state new cases by 206% (22.35 thousand), showing the high impact of the vaccination rate. Also, lowering the monitored vaccination rate $(\alpha_2)$ by 50% increases the steady-state level of new cases by approximately 250% (25.56 thousand). Similarly, decreasing the immunity rate of un-monitored vaccination $(\gamma_1)$ by 10% increases the steady-state new cases by about 300% (29.1 thousand) while reducing the effectiveness of monitored vaccination $(\varepsilon_2)$ by 18% increases the steady-state new cases by 30% (9.528 thousand).

6.   Conclusions

In recent years, the frequent outbreaks of measles in developing and developed countries have become a significant obstacle to achieving the goal of elimination. Mainly, unvaccinated children are the victims of the disease ^[6,7,26]. Despite the Supplementary Immunization Activities (SIAs), most of the population is not vaccinated in some districts of Nepal, including Rautahat, Kapilvastu, Morang, and Bajura. Furthermore, many of those initiating the vaccination does not complete the vaccine doses correctly, making them not completely protected. Because of the resurgence of measles in Nepal since 2017 ^[20,21], WHO's deadline for the elimination goal has been extended to 2023 ^[53]. We developed a novel deterministic model validated by the data from Nepal to evaluate the monitored vaccination programs.

We thoroughly analyzed our model to formulate the vaccinated reproduction number ($R_v > 1$), the stability analysis, and the disease persistence theory. Using the model and the available data from the official websites of WHO, we estimated key parameters related to the un-monitored vaccination. Moreover, we performed a global sensitivity analysis using Latin Hypercube Sampling from the wider parameter space. Our model predicts that the measles elimination goal can be achieved if the monitored vaccination rate is increased by 50% $(\alpha_2 = 0.03)$ or the un-monitored vaccination rate is increased by 6.5% $(\alpha_1 = 0.33)$. The elimination goal can also be achieved by the high effectiveness of the un-monitored vaccine (i.e., with $\varepsilon_1 \geq 0.53$). However, if the current trend continues, our model predicts that the measles will persist causing an obstacle to the measles elimination goal of Nepal. Furthermore, our model predicts that without any additional interventions, measles transmission will continue with a rise in the epidemic after 2033 (Figure 4). Our model also predicts that the major contributor to the measles resurgence in Nepal is susceptible (unvaccinated) and unmonitored vaccinated groups, emphasizing the need for the expedition of monitored-vaccination programs. The epidemic dynamics in Nepal is quite slow, reaching the peak only in about 2097 (Figure 4). The observed slow dynamics is consistent with the previous study ^[72].

We acknowledge some limitations of our study. The parameter estimations are based on the limited data of yearly incidence cases (2000–2019). Moreover, measles is a short-term disease recovered within a month, so the daily or weekly more detailed and accurate data can help improve the prediction of our model. We could obtain the closed form of the unique endemic equilibrium but are unable to perform a detailed analysis of the endemic equilibrium. Instead, we established the disease persistence criteria. The homogeneous mixing assumed among the children of Nepal can be improved by the model with an appropriate network of contact among children.

In summary, we develop a model of measles transmission in the context of Nepal, where monitoring is critical for the successful implementation of the vaccination program. Our thorough analysis and the detailed numerical simulations of the model can provide helpful information for policymakers to design ideal monitored-vaccination programs to achieve the elimination goal of measles from Nepal.

Acknowledgments

This research was supported by the GRAID (Graduate Research Assistantships in Developing Countries) awards from the International Mathematical Union (IMU). AP acknowledges the Nepal Mathematical Society (NMS) for the NMS Ph.D. Fellowship Award 2020. KA acknowledges the Nepal Academy of Science and Technology (NAST) for Ph.D. Fellowship and the Nepal Mathematical Society (NMS) for the NMS Ph.D. Fellowship award 2020. RG acknowledges the University Grants Commission (UGC) for Ph.D. Fellowship and Nepal Mathematical Society (NMS) for the NMS Ph.D. Fellowship Award 2020. The work of NKV was supported by NSF grants DMS-1951793 and DEB-2030479 from the National Science Foundation of the USA and the UGP award from San Diego State University.

Conflict of interest

The authors declare there is no conflict of interest.