Stability and Threshold Dynamics in a Seasonal Mathematical Model for Measles Outbreaks with Double-Dose Vaccination

Abstract

Measles is a highly contagious viral disease that can lead to serious complications, including death, particularly in young children. In this study, we developed a mathematical model that incorporates a seasonal transmission parameter to examine the measles transmission dynamics. We define the basic reproduction number (${\mathcal{R}}_0$) and show its utility as a threshold parameter for global dynamics and the existence of periodic solutions. The model was applied to the measles outbreak that occurred in Pakistan from 2019 to 2021 and provided a good fit to the observed data. Our estimate of the basic reproduction number was found to be greater than one, indicating that the disease will persist in the population. The findings highlight the need to increase vaccination coverage and efficacy to mitigate the impact of the epidemic. The model also shows the long-term behavior of the disease, which becomes endemic and recurs annually. Our simulations demonstrate that a shorter incubation period accelerates the spread of the disease, while a higher vaccination coverage rate reduces its impact. The importance of the second dose of the measles vaccine is emphasized, and a higher vaccine efficacy rate can also help bring ${\mathcal{R}}_0$ below one. Our study provides valuable information for the development and implementation of effective measles control strategies. To prevent future outbreaks, increasing vaccination coverage among the population is the most effective way to reduce the transmission of measles.

1. Introduction

Measles, a highly contagious respiratory illness, continues to threaten global public health despite the availability of a safe and effective vaccine. Caused by the rubeola virus, the disease is easily spread through the air when an infected person sneezes, coughs, or talks, and can remain in the air for up to two hours [1,2]. Close contact with infected bodily fluids, such as saliva, mucus, or blood, can also cause virus transmission. Measles is spread easily from person to person, and an individual can contract the virus by breathing in air contaminated with the virus or by touching a surface contaminated with the virus and then touching their mouth, nose, or eyes [2]. Measles is a viral illness that causes symptoms such as fever, cough, runny nose, and red and watery eyes that usually appear 7–14 days after exposure to the virus. It also causes a red rash on the face that spreads throughout the body and lasts 3–5 days. Other symptoms include muscle pain, fatigue, and loss of appetite. In severe cases, it can lead to complications such as pneumonia, encephalitis, and deafness [1,2,3,4]. Measles can be diagnosed through laboratory tests. Measles is dangerous, especially in vulnerable populations such as young children, pregnant women, and those with weak immune systems. Despite the availability of a safe and cost-effective vaccine, measles remains one of the leading causes of death in young children globally [1,2]; an estimated 139,000 people died from measles in 2019. However, measles continues to be one of the leading causes of death and illness in many parts of the world, particularly in low- and middle-income countries with limited access to healthcare and vaccination. The top ten countries with global measles outbreaks as of January 2023 are India, Yemen, Somalia, Zimbabwe, Pakistan, Ethiopia, Indonesia, Nigeria, Angola, and Afghanistan [5]. These outbreaks were caused by low vaccination coverage, weak health systems, and outbreaks in densely populated urban areas. Vaccination is an effective method to prevent the spread of measles. Two doses of the measles vaccine provide 97% protection [6,7].

Mathematical modeling has played an important role in understanding and predicting the spread of measles. In [8], a mathematical model for the transmission dynamics of the measles epidemic was presented and validated using real incidence data from Pakistan. The study found stability conditions (depending on ${\mathcal{R}}_0$), proposed a strategy to control spread, calculated the sensitivity of ${\mathcal{R}}_0$, and suggested ways to improve vaccine efficacy and coverage. The authors of [9] studied a modified SVEIR measles model with double-dose vaccination to control measles outbreaks in Bangladesh. The model has two equilibria, disease-free and endemic, and the latter persists if ${\mathcal{R}}_0$ remains above one. The local stability analysis supports this result, and a sensitivity analysis shows the transmission rate as the most influential factor. The study concludes with simulations that provide insight to decision makers on how to effectively tackle measles outbreaks. Reference [10] describes a mathematical model to examine the impact of preventive measures on the control of measles. The study shows that the disease-free state is globally stable if the reproduction number ${\mathcal{R}}_0\le 1$, meaning the measles will die out eventually. However, when ${\mathcal{R}}_0>1$, the disease will spread. The authors recommend lowering effective contact with infected individuals and increasing vaccination rates with high-efficiency vaccines to reduce the prevalence of measles. Numerical simulations support these findings. Several autonomous mathematical models have been used to study the transmission dynamics of measles in various countries, with the aim of providing information on the control and prevention of outbreaks [11,12,13,14,15,16].

In recent years, several seasonal mathematical models have been developed to study the measles transmission dynamics and the impact of interventions, such as vaccination. For example, in 2016, a study by Dalziel et al. [17] studied the incidence of measles in 80 major cities in the US and the UK and found that small differences in seasonal patterns caused chaotic patterns in US cities, but not in UK cities. The study highlights the impact of seasonal patterns on disease transmission and the importance of understanding them to prevent persistent chaotic epidemics. Bai and Liu [18] studied a discrete-time measles model with a periodic transmission rate by using the basic reproduction number as the threshold for disease persistence. The parameters were estimated from demographic and epidemiological data, and numerical simulations showed the seasonal fluctuation of measles in China. Xue et al. [19] used the model to investigate the periodic outbreak in mainland China; they recommended enhancing vaccination and implementing an optimal control system to minimize infections. The authors of [20] analyzed the seasonal spread of measles in China using a mathematical model that accounted for periodic transmission fluctuations. The model was used to analyze the dynamics of the disease depending on ${\mathcal{R}}_0$, and to simulate the monthly data of reported cases of measles in China. The results show that increasing immunization, improving vaccine management, and raising public awareness can reduce the incidence of measles.

These studies demonstrate the usefulness of mathematical models in understanding the spread of measles and informing the development of public health interventions to prevent and control the disease. In this paper, we present a mathematical model to study the transmission dynamics of measles with double-dose vaccinations. The model includes a seasonal transmission parameter, which is incorporated to reflect the fluctuations in the number of cases observed during the outbreak in Pakistan from January 2019 to December 2021. To study the dynamics of our time-periodic model, we will employ the techniques used in previous studies on periodic epidemic models [21,22,23,24,25], later applied in several periodic epidemic models (see, e.g., [26,27,28,29,30]).

In the following section, we introduce a periodic model that includes vaccinations and shows the existence of a disease-free periodic solution. In Section 3, we derive the basic reproduction number for our model and demonstrate that depending on its value, either the disease-free periodic solution will be globally stable or the disease will persist in the population. In Section 4, we provide numerical simulations for both cases that support our theoretical results. By fitting monthly measles data from Pakistan from January 2019 to December 2021, we estimate some unknown parameters and determine the basic reproduction number. We will then conduct numerical evaluations to examine the impact of the seasonal contact rate, increasing vaccination coverage, vaccine efficacy, and other significant parameters on the transmission dynamics of measles. Finally, we present our discussion in Section 5.

2. Seasonal Model for Measles Transmission

We split the human population into six distinct compartments to study the dynamics of measles transmission: susceptibles $S\left(t\right)$, individuals who received the first dose of the MMR vaccine $V_1\left(t\right)$, individuals who received the second dose of the MMR vaccine $V_2\left(t\right)$, exposed or asymptomatic $E\left(t\right)$, symptomatic or infectious $I\left(t\right)$, and recovered individuals $R\left(t\right)$. The total size of the human population at any given time t is represented by $N\left(t\right)=S\left(t\right)+V_1\left(t\right)+V_2\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right).$

Figure 1 shows the movement of individuals between compartments; our model equations are given by

$$\begin{array}{cc}\hfill \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}t=& \Lambda -\beta \left(t\right)I\left(t\right)S\left(t\right)-(\theta +\mu )S\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}{V}_1\left(t\right)}{\mathrm{d}t}=& \theta S\left(t\right)-(1-\alpha )\beta \left(t\right)I\left(t\right)V_1\left(t\right)-({\sigma }_1+\mu )V_1\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}{V}_2\left(t\right)}{\mathrm{d}t}=& {\sigma }_1{V}_1\left(t\right)-({\sigma }_2+\mu )V_2\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}=& \beta \left(t\right)I\left(t\right)S\left(t\right)+(1-\alpha )\beta \left(t\right)I\left(t\right)V_1\left(t\right)-(\nu +\mu )E\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \nu E\left(t\right)-\gamma I\left(t\right)-(\delta +\mu )I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}R\left(t\right)}{\mathrm{d}t}=& \gamma I\left(t\right)+{\sigma }_2{V}_2\left(t\right)-\mu R\left(t\right),\hfill \end{array}$$

(1)

where $\beta \left(t\right)$ denotes the time-dependent effective contact rate. In (1), we assume that $\beta \left(t\right)$ is a continuous and positive $\omega$-periodic function. The human recruitment rate and the natural death rate are indicated by $\Lambda$ and $\mu$, respectively. Susceptible individuals who received the first dose of the vaccine move at a rate of $\theta$ to the vaccinated compartment $V_1$. Since the first dose of the vaccine is ineffective in protecting against measles, individuals vaccinated with only one dose may become infected by contact with symptomatic individuals at a rate of $(1-\alpha )\beta \left(t\right)$, while the rest of the population moves to the class of those vaccinated with two doses at a rate of ${\sigma }_1$. It should be noted that $\alpha =0$ indicates a vaccine that does not provide protection at all, while $\alpha =1$ represents a perfect vaccine ($0<\alpha <1$). At a rate of ${\sigma }_2$, individuals vaccinated with two doses are transferred to the recovered class. A person who contracts the disease will progress from susceptible (S) to exposed (E) at the rate of $\beta \left(t\right)SI$. Following the incubation period, an exposed individual progresses at a rate of $\nu$ to the symptomatically infected class (I). Following the period of infection, recovered individuals enter the category of R. For those in compartment I, there is an additional disease-induced death rate, denoted by $\delta$. The description of the model parameters can be found in

Table 1. Description of parameters for model (1).

Parameters	Description
$\Lambda$	Recruitment rate of susceptible humans
$\mu$	Natural death rate
$\beta \left(t\right)$	Periodic transmission rate of symptomatic infection
$\alpha$	Efficacy of vaccine
${\sigma }_1$	Rate of vaccination with the second dose
${\sigma }_2$	Progression rate from $V_2$ to R
$\theta$	Coverage rate of vaccines for the entire population
$\nu$	Incubation rate
$\gamma$	Recovery rate
$\delta$	Disease-induced death rates for humans

.

Clearly, any solution to system (1) starting from non-negative initial values is non-negative, so we discuss the bounds of the solution of this system first. Define

$$\varphi =\left(S\left(0\right),V_1\left(0\right),V_2\left(0\right),E\left(0\right),I\left(0\right),R\left(0\right)\right)\in {\mathbb{R}}_{+}^6.$$

(2)

In case the disease is not present in the population, for the total human population $N\left(t\right)$ with a positive initial condition $\varphi \in {\mathbb{R}}_{+}^6$, we have the equation

$$N\left(t\right)t=\Lambda -(\theta +\mu )N\left(t\right),$$

(3)

from which we obtain the following:

$$N\left(t\right)=N\left(0\right)e^{-(\theta +\mu )t}+\frac{\Lambda }{\theta +\mu }(1-e^{-(\theta +\mu )t}),$$

(4)

with an arbitrary initial value $N\left(0\right)$. Equation (4) has a unique equilibrium $N^{*}=\frac{\Lambda }{\theta +\mu }$ in ${\mathbb{R}}_{+}$. Consequently, $|N\left(t\right)-N^{*}|\to 0$ as $t\to \infty$ and $N^{*}$ is globally attractive in ${\mathbb{R}}_{+}$.

Letting $\Omega \left\{\left(S,V_1,V_2,E,I,R\right)\in {\mathbb{R}}_{+}^6:S+V_1+V_2+E+I+R⩽N^{*}\right\}$, we have the following result.

Proposition 1. 

The region Ω is positively invariant with respect to the system (1). In particular, $\left(S\left(t\right),V_1\left(t\right),V_2\left(t\right),E\left(t\right),I\left(t\right),R\left(t\right)\right)$ is positive for all $t>0$ if the initial values $\varphi >0$ at $t=0$.

Proof. 

To prove that solutions of system (1) are non-negative, it is sufficient to demonstrate that for non-negative initial values of $\varphi \in {\mathbb{R}}_{+}^6$, the solution $(S\left(t\right),V_1\left(t\right),V_2\left(t\right),E\left(t\right),I\left(t\right),R\left(t\right))$ remains positive for all $t>0$. Let

$$n\left(t\right)=min\{S\left(t\right),V_1\left(t\right),V_2\left(t\right),E\left(t\right),I\left(t\right),R\left(t\right)\},\forall t>0.$$

Clearly, $n\left(0\right)>0$, and let us suppose there is some $t_1>0$, such that $n\left(t_1\right)=0$ and $n\left(t\right)>0,\forall t\in [0,t_1)$. If we also assume that $n\left(t_1\right)=S\left(t_1\right)$, then using the first equation of system (1), we can conclude that

$$\frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}\ge -(\beta \left(t\right)I\left(t\right)+\theta +\mu )S\left(t\right),\forall t\in [0,t_1].$$

Then,

$$0=S\left(t_1\right)\ge S\left(0\right)e^{-{\int }_0^{t_1}(\beta \left(s\right)I\left(s\right)+\theta +\mu )ds}>0.$$

The statement leads to a contradiction, indicating that $S\left(t\right)$ cannot become negative for any $t>0$. We obtain a similar result for $V_1\left(t\right)$ and $V_2\left(t\right)$ in an analogous way.

As for the last three compartments of (1), suppose that there exists a minimal $t>0$, such that one of $E\left(t\right)$, $I\left(t\right)$, and $R\left(t\right)$ is equal to zero. Suppose first that this compartment is $E\left(t\right)$. If $n\left(t_1\right)=E\left(t_1\right)$, since $I\left(t\right)\ge 0$, $V_1\left(t\right)\ge 0$ and $S\left(t\right)\ge 0$ for all $t\in [0,t_1]$, from the fourth equation of (1), it follows that

$$\frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}\ge -(\nu +\mu )E\left(t\right),\forall t\in [0,t_1].$$

Thus,

$$0=E\left(t_1\right)\ge E\left(0\right)e^{-(\nu +\mu )t_1}>0,$$

which also leads to a contradiction. We can reach the same contradictory conclusion if we assume that either $n\left(t_1\right)=I\left(t_1\right)$ or $n\left(t_1\right)=R\left(t_1\right)$. Thus, it follows that $S\left(t\right)\ge 0$, $V_1\left(t\right)\ge 0$, $V_2\left(t\right)\ge 0$, $E\left(t\right)\ge 0$, $I\left(t\right)\ge 0$, and $R\left(t\right)\ge 0$ for all $t>0$ where equality holds if the initial values equal to zero.

In relation to Equation (4), we can express $N\left(t\right)$ as $N\left(0\right)e^{-(\theta +\mu )t}+\frac{\Lambda }{\theta +\mu }(1-e^{-(\theta +\mu )t})$, where $N\left(0\right)=S\left(0\right)+V_1\left(0\right)+V_2\left(0\right)+E\left(0\right)+I\left(0\right)+R\left(0\right)$. As a result, $N\left(t\right)$ remains bounded for all $t\ge 0$ and

$$\underset{t\to \infty }{lim\; sup}N\left(t\right)=\frac{A}{\theta +\mu }.$$

This indicates that $S\left(t\right)$, $V_1\left(t\right)$, $V_2\left(t\right)$, $E\left(t\right)$, $I\left(t\right)$, and $R\left(t\right)$ are also bounded for all t greater than zero. Thus, this concludes the proof. □

It is evident that system (1) possesses a singular disease-free equilibrium, denoted as

$$\begin{array}{cc}\hfill E_0=& \left(S^{*},V_1^{*},V_2^{*},0,0,R^{*}\right)\hfill \\ \hfill =& \left(\frac{\Lambda }{\theta +\mu },\frac{\theta \Lambda }{(\theta +\mu )({\sigma }_1+\mu )},\frac{{\sigma }_1\theta \Lambda }{(\theta +\mu )({\sigma }_1+\mu )({\sigma }_2+\mu )},0,0,\frac{{\sigma }_1\theta \Lambda }{\mu (\theta +\mu )({\sigma }_1+\mu )({\sigma }_2+\mu )}\right).\hfill \end{array}$$

3. Threshold Dynamics

We describe the dynamics of our model depending on the basic reproduction number. First, we derive the basic reproduction number ${\mathcal{R}}_0$ by using the general procedure presented by Wang and Zhao [24]. We show that the unique disease-free equilibrium $E_0$ is globally asymptotically stable and that the disease goes extinct if the basic reproduction number ${\mathcal{R}}_0<1$. We also investigate the persistence of the disease and prove the existence of a positive periodic solution of model (1) if the basic reproduction number ${\mathcal{R}}_0>1$.

3.1. Local Stability

As per the method and techniques laid out by Wang and Zhao [24], we demonstrate the local stability of the disease-free periodic equilibrium, $E_0$, in model (1), based on the basic reproduction number, ${\mathcal{R}}_0$.

Linearizing the system (1) at the disease-free equilibrium $E_0$, we obtain the equations for the exposed and infectious human population, respectively:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}=& \beta \left(t\right)I\left(t\right)S^{*}+(1-\alpha )\beta \left(t\right)I\left(t\right)V_1^{*}-(\nu +\mu )E\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \nu E\left(t\right)-\gamma I\left(t\right)-(\delta +\mu )I\left(t\right).\hfill \end{array}$$

To begin with, we can check that model (1) meets the conditions (A1)–(A7) stated in [24]. We introduce the $2\times 2$ matrix functions $F\left(t\right)$ and $V\left(t\right)$ given as

$$\begin{array}{cc}\hfill F\left(t\right)=& \left[\begin{array}{cc}0& \beta \left(t\right)S^{*}+(1-\alpha )\beta \left(t\right)V_1^{*}\\ 0& 0\end{array}\right],V\left(t\right)=\left[\begin{array}{cc}\nu +\mu & 0\\ -\nu & \gamma +\delta +\mu \end{array}\right].\hfill \end{array}$$

Suppose $\mathcal{Z}(t,s),t⩾s$, is the $2\times 2$ matrix solution of the following initial value problem

$$\left\{\begin{array}{c}\frac{\mathrm{d}\mathcal{Z}(t,s)}{\mathrm{d}t}=-V\left(t\right)\mathcal{Z}(t,s),\forall t⩾s,\hfill \\ \mathcal{Z}(s,s)=I,\hfill \end{array}\right.$$

where I stands for the $2\times 2$ identity matrix.

Let us assume that ${\mathcal{C}}_{\omega }$ is the ordered Banach space of $\omega$-periodic functions from $\mathbb{R}$ to ${\mathbb{R}}^2$, with the usual maximum norm ${\parallel ·\parallel }_{\infty }$ and introduce the positive cone

$${\mathcal{C}}_{\omega }^{+}\{\varphi \in {\mathcal{C}}_{\omega }:\varphi \left(t\right)⩾0,\forall t\in \mathbb{R}\}.$$

Let $\varphi \left(s\right)\in {\mathcal{C}}_{\omega }^{+}$ be the initial distribution of infected individuals, periodic with period $\omega$ in s. The quantity $F\left(s\right)\varphi \left(s\right)$ represents the rate of new infections caused by those infected at time s. For $t⩾s$, the expression $\mathcal{Z}(t,s)F\left(s\right)\varphi \left(s\right)$ gives us the distribution of individuals who were infected at time s and still infected at time t. Therefore,

$$\psi \left(t\right) :={\int }_{-\infty }^t\mathcal{Z}(t,s)F\left(s\right)\varphi \left(s\right)ds={\int }_0^{\infty }\mathcal{Z}(t,t-a)F(t-a)\varphi (t-a)da,$$

provides the distribution of the total of new infections at time t originating from all infected individuals $\varphi \left(s\right)$ introduced at any previous time $s⩽t$. Following that, we define the linear next infection operator $\mathcal{L}:{\mathcal{C}}_{\omega }\to {\mathcal{C}}_{\omega }$ as follows

$$\left(\mathcal{L}\varphi \right)\left(t\right)={\int }_0^{\infty }\mathcal{Z}(t,t-a)F(t-a)\varphi (t-a)da,\forall t\in \mathbb{R},\varphi \in {\mathcal{C}}_{\omega }.$$

(5)

Following [24], we define the basic reproduction number of (1) as the spectral radius $\rho \left(\mathcal{L}\right)$ of the operator $\mathcal{L}$; therefore, ${\mathcal{R}}_0\rho \left(\mathcal{L}\right)$.

With reference to Theorem 2.1 and Theorem 2.2 presented in [24], we derive the following results regarding ${\mathcal{R}}_0$ and the locally asymptotically stable disease-free equilibrium $E_0$ of (1).

Theorem 1. 

The following statements are valid:

1. 

${\mathcal{R}}_0=1$ if and only if $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)=1$;

2. 

${\mathcal{R}}_0>1$ if and only if $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)>1$;

3. 

${\mathcal{R}}_0<1$ if and only if $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)<1$,

where ${\Phi }_{F-V}\left(t\right)$ is the monodromy matrix of the linear ω-periodic system $\frac{\mathrm{d}u}{\mathrm{d}t}=\left[F\right(t)-V(t\left)\right]u.$

The following theorem, based on the previous discussion, addresses the local stability of the disease-free solution $E_0$ of (1).

Theorem 2. 

If ${\mathcal{R}}_0<1$, the disease-free solution $E_0$ of (1) is locally asymptotically stable and unstable if ${\mathcal{R}}_0>1$.

Proof. 

The following is the Jacobian matrix of (1) calculated at $E_0$

$$J\left(t\right)=\left[\begin{array}{cc}F\left(t\right)-V\left(t\right)& 0\\ A\left(t\right)& B\end{array}\right],$$

where

$$\begin{array}{cc}\hfill A\left(t\right)=& \left[\begin{array}{cc}0& \beta \left(t\right)S^{*}\\ 0& (1-\alpha )\beta \left(t\right)V_1^{*}\\ 0& 0\\ 0& 0\end{array}\right]\mathrm{and}B=\left[\begin{array}{cccc}-(\theta +\mu )& 0& 0& 0\\ \theta & -({\sigma }_1+\mu )& 0& 0\\ 0& {\sigma }_1& -({\sigma }_2+\mu )& 0\\ 0& 0& {\sigma }_2& -\mu \end{array}\right].\hfill \end{array}$$

According to [31], $E_0$ is LAS if $\rho \left({\Phi }_B\left(\omega \right)\right)<1$ and $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)<1$. B is a constant matrix, and its eigenvalues are ${\lambda }_1=-(\theta +\mu )<0$, ${\lambda }_2=-({\sigma }_1+\mu )<0$, ${\lambda }_3=-({\sigma }_2+\mu )<0$ and ${\lambda }_4=-\mu <0$. Since ${\lambda }_i$, $i=1,2,3,4$ are negative, we have $\rho \left({\Phi }_B\right)<1$. Therefore, the stability of $E_0$ depends on $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)$. Thus, $E_0$ is locally asymptotically stable if $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)<1$, and unstable if $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)>1$. Thus, we finish the proof using Theorem 1. □

3.2. Global Stability of the Disease-Free Equilibrium

Theorem 3. 

If ${\mathcal{R}}_0<1$, then the disease-free periodic solution $E_0$ of (1) is globally asymptotically stable.

Proof. 

From Theorem 2, it is clear that if ${\mathcal{R}}_0>1$, then $E_0$ is unstable and $E_0$ is locally asymptotically stable when ${\mathcal{R}}_0<1$. Therefore, it is simply sufficient to demonstrate that $E_0$ is globally attractive when ${\mathcal{R}}_0<1$. By Theorem 1, it flows that ${\mathcal{R}}_0<1$ if and only if $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)<1$. Since $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)$ is continuous, we can choose $\epsilon >0$ small enough such that $\rho \left({\Phi }_{M_{\epsilon }}\left(\omega \right)\right)<1$, the $2\times 2$ matrix function $M_{\epsilon }\left(t\right)$ given by

$$M_{\epsilon }\left(t\right)=\left[\begin{array}{cc}-\nu -\mu & \beta \left(t\right)(S^{*}+\epsilon )+(1-\alpha )\beta \left(t\right)(V_1^{*}+\epsilon )\\ \nu & -\gamma -\delta -\mu \end{array}\right].$$

From the first equation of (1), we have

$$\begin{array}{cc}\hfill \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}=& \Lambda -\beta \left(t\right)I\left(t\right)S\left(t\right)-(\theta +\mu )S\left(t\right).\hfill \end{array}$$

Because $I\left(t\right)\ge 0$ from Proposition 1, it can be shown that

$$\frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}\le \Lambda -(\theta +\mu )S(t),$$

which implies that

$$\underset{t\to \infty }{lim\; sup}S\left(t\right)\le \frac{\Lambda }{\theta +\mu }=S^{*}.$$

Thus, for an arbitrary positive $\epsilon$, there is $t_1>0$, such that $S\left(t\right)\le S^{*}+\epsilon$ for all $t>t_1$. A similar analysis can be done for $V_1\left(t\right)$, from the second equation of (1) and Proposition 1, we have

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{V}_1\left(t\right)}{\mathrm{d}t}=& \theta S\left(t\right)-(1-\alpha )\beta \left(t\right)I\left(t\right)V_1\left(t\right)-({\sigma }_1+\mu )V_1\left(t\right)\hfill \\ \hfill \le & \theta S\left(t\right)-({\sigma }_1+\mu )V_1\left(t\right),\hfill \end{array}$$

which implies that

$$\underset{t\to \infty }{lim\; sup}{V}_1\left(t\right)\le \frac{\theta S^{*}}{{\sigma }_1+\mu }=\frac{\theta \Lambda }{(\theta +\mu )({\sigma }_1+\mu )}=V_1^{*}.$$

Hence, there exists $t_1>0$, such that $V_1\left(t\right)\le V_1^{*}+\epsilon$ for all $t>t_1$, for an arbitrary positive $\epsilon$. From system (1), for $t>t_1$, we obtain

$$\begin{array}{cc}\hfill \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}\le & \beta \left(t\right)I\left(t\right)(S^{*}+\epsilon )+(1-\alpha )\beta \left(t\right)I\left(t\right)(V_1^{*}+\epsilon )-(\nu +\mu )E\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}\le & \nu E\left(t\right)-\gamma I\left(t\right)-(\nu +\mu )I\left(t\right),\hfill \end{array}$$

utilizing the comparison system described below

$$\begin{array}{cc}\hfill \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}=& \beta \left(t\right)I\left(t\right)(S^{*}+\epsilon )+(1-\alpha )\beta \left(t\right)I\left(t\right)(V_1^{*}+\epsilon )-(\nu +\mu )E\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \nu E\left(t\right)-\gamma I\left(t\right)-(\nu +\mu )I\left(t\right),\hfill \end{array}$$

(6)

system (6) can be rewritten as

$$\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}=M_{\epsilon }\left(t\right)U\left(t\right),$$

(7)

where $U\left(t\right)={\left(E\left(t\right),I\left(t\right)\right)}^T$.

From ([25], Lemma 2.1), there is an $\omega$-periodic positive function $p\left(t\right)$, such that $U\left(t\right)=p\left(t\right)e^{\xi t}$ is a solution of (7) and $\xi =\frac{1}{\omega }ln\rho \left({\Phi }_{M_{\epsilon }}\left(\omega \right)\right)<0$ since $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)<1$ when ${\mathcal{R}}_0<1$. Therefore, we have $U\left(t\right)\to 0$ as $t\to \infty$, which implies that the zero solution of system (6) is globally asymptotically stable. For any ${(E\left(0\right),I\left(0\right))}^T\in {\mathbb{R}}_{+}^2$, we can choose $n^{*}>0$ s.t. ${(E\left(0\right),I\left(0\right))}^T\le n^{*}p\left(0\right)$. By utilizing the comparison principle ([32], Theorem B.1), we have ${(E\left(t\right),I\left(t\right))}^T\le n^{*}U\left(t\right)$ for all $t>0$, where $n^{*}U\left(t\right)$ is also the solution of system (7). Hence, we obtain

$$\underset{t\to \infty }{lim}{\left(E\left(t\right),I\left(t\right)\right)}^T={\left(0,0\right)}^T.$$

From the equations of $S^{\prime }\left(t\right),V_1^{\prime }\left(t\right),V_2^{\prime }\left(t\right),R^{\prime }\left(t\right)$ in system (1) and considering the limit system, we have

$$\begin{array}{cc}\hfill S^{\prime }\left(t\right)=& \Lambda -(\theta +\mu )S\left(t\right),\hfill \\ \hfill V_1^{\prime }\left(t\right)=& \theta S\left(t\right)-({\sigma }_1+\mu )V_1\left(t\right),\hfill \\ \hfill V_2^{\prime }\left(t\right)=& {\sigma }_1{V}_1\left(t\right)-({\sigma }_2+\mu )V_2\left(t\right),\hfill \\ \hfill R^{\prime }\left(t\right)=& {\sigma }_2{V}_2\left(t\right)-\mu R\left(t\right).\hfill \end{array}$$

By the theory of asymptotically autonomous systems [33], and system (1), we have ${lim}_{t\to \infty }S\left(t\right)=S^{*},{lim}_{t\to \infty }{V}_1\left(t\right)=V_1^{*},{lim}_{t\to \infty }{V}_2\left(t\right)=V_2^{*}$ and ${lim}_{t\to \infty }R\left(t\right)=R^{*},$ and the proof is complete. □

3.3. Existence of Positive Periodic Solutions

Define

$$\begin{array}{cc}\hfill \mathcal{X} :=& \left\{(S,V_1,V_2,E,I,R)\in {\mathbb{R}}_{+}^6\right\},\hfill \\ \hfill {\mathcal{X}}_0 :=& \left\{(S,V_1,V_2,E,I,R)\in \mathcal{X}:E>0,I>0\right\},\hfill \end{array}$$

and

$$\partial {\mathcal{X}}_0 :=\mathcal{X}\backslash {\mathcal{X}}_0=\left\{(S,V_1,V_2,E,I,R)\in \mathcal{X}:E=0\mathrm{or}I=0\right\}.$$

Let $P:{\mathbb{R}}_{+}^6\to {\mathbb{R}}_{+}^6$ denote the Poincaré map corresponding to (1), then P is given by

$$P\left(x^0\right)=u(\omega ,x^0),\mathrm{for}{x}^0\in {\mathbb{R}}_{+}^6,$$

where $u(t,x^0)$ is the unique solution of (1) with the initial condition $x^0\in \mathcal{X}$.

Lemma 1. 

If ${\mathcal{R}}_0>1$, then there exists a $\sigma >0$, such that for any $\varphi \in {\mathcal{X}}_0$ defined in (2) with $\parallel \varphi -E_0\parallel ⩽\sigma$, we obtain

$$\underset{m\to \infty }{lim\; sup}d\left(P^m\left(\varphi \right),E_0\right)⩾\sigma .$$

Proof. 

If ${\mathcal{R}}_0>1$, we know from Theorem 1 that $\rho \left({\Phi }_{F-V}\left(\omega \right)\right)>1$. Therefore, we can select $\kappa >0$ small enough such that we have $\rho \left({\Phi }_{F-V-M_{\kappa }}\left(\omega \right)\right)>1$, where $M_{\kappa }\left(t\right)$ is the $2\times 2$ matrix function given by

$$M_{\kappa }\left(t\right)=\left[\begin{array}{cc}0& \beta \left(t\right)\kappa +(1-\alpha )\beta \left(t\right)\kappa \\ 0& 0\end{array}\right].$$

In the following, we claim

$$\underset{m\to \infty }{lim\; sup}d\left(P^m\left(\varphi \right),E_0\right)⩾\sigma .$$

(8)

By contradiction, suppose that (8) does not hold. Then,

$$\underset{m\to \infty }{lim\; sup}d\left(P^m\left(\varphi \right),E_0\right)<\sigma ,$$

(9)

for some $\varphi \in {\mathcal{X}}_0$. Without loss of generality, we may assume

$$d\left(P^m\left(\varphi \right),E_0\right)<\sigma ,\forall m\ge 0.$$

By the continuous dependence of the solutions w.r.t. the initial values, we have

$$\parallel u(t_1,P^m\left(\varphi \right)-u(t_1,E_0)\parallel <\kappa ,\forall m\ge 0,t_1\in \left[0,\omega \right].$$

For any $t⩾0,$ let $t=m\omega +t_1,$ where $t_1\in \left[0,\omega \right)$ and $m=\left[\frac{t}{\omega }\right]$, which is the largest integer less than or equal to $\frac{t}{\omega }$. Then, we obtain

$$\parallel u\left(t,\varphi \right)-u\left(t,E_0\right)\parallel =\parallel u\left(t_1,P^m\left(\varphi \right)\right)-u(t_1,E_0)\parallel <\kappa ,$$

for all $t⩾0$, which implies that $S^{*}-\kappa \le S\left(t\right)\le S^{*}+\kappa$, $V_1^{*}-\kappa \le V_1\left(t\right)\le V_1^{*}+\kappa$, $V_2^{*}-\kappa \le V_2\left(t\right)\le V_2^{*}+\kappa$, $0\le E\left(t\right)\le \kappa$ and $0\le I\left(t\right)\le \kappa$. Then for $\parallel \varphi -E_0\parallel ⩽\sigma$, we obtain

$$\begin{array}{cc}\hfill \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}\ge & \beta \left(t\right)I\left(t\right)(S^{*}-\kappa )+(1-\alpha )\beta \left(t\right)I\left(t\right)(V_1^{*}-\kappa )-(\nu +\mu )E\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}\ge & \nu E\left(t\right)-\gamma I\left(t\right)-(\nu +\mu )I\left(t\right),\hfill \end{array}$$

Consider the auxiliary linear system

$$\begin{array}{cc}\hfill \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}=& \beta \left(t\right)I\left(t\right)(S^{*}-\kappa )+(1-\alpha )\beta \left(t\right)I\left(t\right)(V_1^{*}-\kappa )-(\nu +\mu )E\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \nu E\left(t\right)-\gamma I\left(t\right)-(\nu +\mu )I\left(t\right).\hfill \end{array}$$

(10)

By ([25], Lemma 2.1) there exists a positive, $\omega$-periodic function $p_2\left(t\right)$, such that $h\left(t\right)=(E\left(t\right),I\left(t\right))=e^{{\xi }_{2}t}{p}_2\left(t\right)$ is a solution of (10) and ${\xi }_2=\frac{1}{\omega }ln\rho \left({\Phi }_{F-V+M_{\kappa }}\left(\omega \right)\right)>0$. As ${\mathcal{R}}_0>1$ and $\rho \left({\Phi }_{F-V-M_{\kappa }}\left(\omega \right)\right)>1$, if $h\left(0\right)>0$, $h\left(t\right)\to \infty$ as $t\to \infty$. Applying the comparison principle ([32], Theorem B.1), we have $E\left(0\right)>0$, $I\left(0\right)>0$, ${lim}_{t\to \infty }E\left(t\right)=\infty$ and ${lim}_{t\to \infty }I\left(t\right)=\infty$. This is a contradiction to $E\left(t\right)<\kappa$, $I\left(t\right)<\kappa$ and Equation (9). Therefore, Equation (8) is true and the proof is complete. □

Theorem 4. 

Assume that ${\mathcal{R}}_0>1$. Then system (1) has at least one positive periodic solution.

Proof. 

We need to prove that P is uniformly persistent w.r.t. $({\mathcal{X}}_0,\partial {\mathcal{X}}_0)$, as from this, applying ([34], Theorem 3.1.1), it follows that the solution of (1) is uniformly persistent with respect to $({\mathcal{X}}_0,\partial {\mathcal{X}}_0)$. First, we show that ${\mathcal{X}}_0$ and $\partial {\mathcal{X}}_0$ are positively invariant in regard to the system (1). For $\varphi \in {\mathcal{X}}_0$, solving (1) for all $t>0$, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill S\left(t\right)=& e^{{\int }_0^t-(\beta \left(s\right)I\left(s\right)+\theta +\mu )ds}\left[S\left(0\right)+\Lambda {\int }_0^t{e}^{{\int }_0^s(\beta \left(r\right)I\left(r\right)+\theta +\mu )dr}ds\right]>0,\hfill \end{array}\end{array}$$

(11)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E\left(t\right)=& e^{-(\nu +\mu )t}\left[E\left(0\right)+{\int }_0^t\left(S\left(s\right)+(1-\alpha )V_1\left(s\right)\right)\beta \left(s\right)I\left(s\right)e^{(\nu +\mu )s}ds\right]>0,\hfill \end{array}\end{array}$$

(12)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill I\left(t\right)=& e^{-(\gamma +\nu +\mu )t}\left[I\left(0\right)+\nu {\int }_0^{t}E\left(s\right)e^{(\gamma +\nu +\mu )s}ds\right]>0,\hfill \end{array}\end{array}$$

(13)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill R\left(t\right)=& e^{-\mu t}\left[R\left(0\right)+{\int }_0^t\left(\gamma I\left(s\right)+{\sigma }_2{V}_2\left(s\right)\right)e^{\mu s}ds\right]>0.\hfill \end{array}\end{array}$$

(14)

Thus, ${\mathcal{X}}_0$ is a positively invariant set. Since $\mathcal{X}$ is also positively invariant and $\partial {\mathcal{X}}_0$ is relatively closed in $\mathcal{X}$, it gives $\partial {\mathcal{X}}_0$ as positively invariant. Moreover, from Proposition 1, it follows that system (1) is point dissipative. Let us introduce

$$M_{\partial }=\left\{x^0\in \partial {\mathcal{X}}_0:P^m\left(x^0\right)\in \partial {\mathcal{X}}_0,\forall m⩾0\right\},$$

where $x^0=\varphi$. We first show that

$$M_{\partial }=\left\{\left(S\left(0\right),V_1\left(0\right),V_2\left(0\right),E\left(0\right),I\left(0\right),R\left(0\right)\right):S>0,V_1>0,V_2>0,R>0\right\}.$$

(15)

Let us note that $M_{\partial }\supseteq \left\{\left(S\left(0\right),V_1\left(0\right),V_2\left(0\right),E\left(0\right),I\left(0\right),R\left(0\right)\right):S>0,V_1>0,V_2>0,R>0\right\}$. It suffices to prove that

$$M_{\partial }\subset \left\{\left(S,V_1,V_2,0,0,R\right):S>0,V_1>0,V_2>0,R>0\right\},$$

for arbitrary initial condition $\varphi \in \partial {\mathcal{X}}_0$, such that $E\left(n\omega \right)=0$ or $I\left(n\omega \right)=0$, for all $n⩾0$.

Let us assume, by way of contradiction, that there exists an integer $n_1⩾0$, such that $E\left(n_1\omega \right)>0$ and $I\left(n_1\omega \right)>0$. Then, by putting $t=n_1\omega$ into the place of the initial time $t=0$ in (11)–(14), we have $S\left(t\right)>0$, $V_1\left(t\right)>0$, $V_2\left(t\right)>0$, $E\left(t\right)>0$, $I\left(t\right)>0$, $R\left(t\right)>0$. This is in contradiction with the positive invariance of $\partial {\mathcal{X}}_0$.

From Lemma 1, P is weakly uniformly persistent w.r.t. $({\mathcal{X}}_0,\partial {\mathcal{X}}_0)$. Proposition 1 guarantees the existence of a global attractor of P. Then $E_0$ is an isolated invariant set in $\mathcal{X}$ and $W^s\left(E_0\right)\cap {\mathcal{X}}_0=\varnothing$. Each solution in $M_{\partial }$ tends to $E_0$; $E_0$ is clearly acyclic in $M_{\partial }$. Based on Theorem 1.3.1 and Remark 1.3.1 in [34], we can deduce that P is uniformly (strongly) persistent w.r.t. $({\mathcal{X}}_0,\partial {\mathcal{X}}_0)$. Hence, there exists an $\epsilon >0$, such that

$$\underset{t\to \infty }{lim\; inf}{\left(E\left(t\right),I\left(t\right)\right)}^T⩾{\left(\epsilon ,\epsilon \right)}^T,$$

for all $\varphi \in {\mathcal{X}}_0$. According to Theorem 1.3.6 in [34], P has a fixed point $\overline{\varphi }\in {\mathcal{X}}_0$ and, hence, system (1) has at least one periodic solution $u(t,\overline{\varphi })$ with $\overline{\varphi }=(\overline{S}\left(0\right),{\overline{V}}_1\left(0\right),{\overline{V}}_2\left(0\right),\overline{E}\left(0\right),\overline{I}\left(0\right),$$\overline{R}\left(0\right))\in {\mathcal{X}}_0.$ Now, let us prove that $\overline{S}\left(0\right)$ is positive. If $\overline{S}\left(0\right)=0$, then we have that $\overline{S}\left(0\right)>0$ for all $t>0$. However, using the periodicity of the solution, we have $\overline{S}\left(0\right)=\overline{S}\left(n\omega \right)=0$, which is a contradiction. □

4. Case Study: Measles Outbreak in Pakistan

In this part of the article, we utilize our model to look at how measles spread in Pakistan during the outbreak from January 2019 to December 2021. To show how well model (1) with periodic parameters matches seasonal fluctuation data, the simulation results are shown. Here, following [19,20,35,36], we assume that the transmission rate $\beta \left(t\right)$ is a time-periodic function with 12 months as a period, given by $\beta \left(t\right)={\beta }_0·\left(a+bsin(c+\frac{2\pi t}{12})\right),$ where $a,b$ are free adjustment parameters, c is the magnitude of forcing, and ${\beta }_0$ is the baseline value of the time-dependent contact rate.

4.1. Parameters Estimation and Curve Fitting

According to the World Bank [37], the total population in Pakistan in 2021 was $N\left(0\right)=231,402,117$, with a total life expectancy at birth of 66 years, so $\mu =\frac{1}{12\times 66}$ and the monthly human birth population can be computed as $\Lambda =\mu \times N\left(0\right)\approx 292,174$. However, according to the World Health Organization (WHO), the World Bank, and GAVI [38,39,40,41], Pakistan’s coverage of the first dose of the measles vaccine was around 81%, and the coverage of the second dose of the measles vaccine was around 55%. Hence, we set $V_1\left(0\right)=0.81\times N\left(0\right)$ and $V_2\left(0\right)=0.55\times N\left(0\right)$. The rest of the initial values were set to be $S\left(0\right)=231,402,117,E\left(0\right)=500,I\left(0\right)=286$, and $R\left(0\right)=231,532$. One dose of the MMR vaccine is approximately 93% effective at preventing measles, whereas two doses of the vaccine are about 97% effective [6,7]. Therefore, it is assumed that $\alpha$, the value representing the efficacy of the vaccine, is equal to $\alpha =0.93$.

We employed Latin hypercube sampling along with the least squares method to estimate the parameters of model (1), which provides the best appropriate fit to the data. This sampling method is used to assess the variance of many parameter values concurrently (see [42] for details). Model (1) was fitted to the observed number of measles cases in Pakistan between January 2019 and December 2021 [43], allowing us to obtain an estimate for the measles model parameters. The parameter values of model (1) were taken from the literature to help us make estimates, and the rest of the parameters were determined by fitting the model to the data. Figure 2 presents the model fit to the measles data from Pakistan, showing a reasonably good fit and a high goodness of fit, $R^2=0.749996$, indicating that it can capture the key patterns of the measles epidemic incidences between 2019 and 2021.

Table 2. Parameters for model (1) providing the best fit.

Parameters	Values	Units	Sources
$\Lambda$	292174	Month${}^{-1}$	[37,38]
$\mu$	$\frac{1}{12\times 66}$	Month${}^{-1}$	[37,38]
${\beta }_0$	$2.02\times {10}^{-8}$	Month${}^{-1}$	Fitted
a	0.92	–	Fitted
b	1	–	Fitted
c	7.9	–	Fitted
$\alpha$	0.93	Month${}^{-1}$	[6,7]
${\sigma }_1$	0.0033	Month${}^{-1}$	[45]
${\sigma }_2$	0.028	Month${}^{-1}$	Fitted
$\theta$	0.000154	Month${}^{-1}$	Fitted
$\nu$	0.57	Month${}^{-1}$	[1,2]
$\gamma$	3.78	Month${}^{-1}$	[1,2]
$\delta$	0.0465	Month${}^{-1}$	[46]

displays the parameter values that provide the best fit in Figure 2. Given the provided parameter values, we can use the method established in [44] to estimate Pakistan’s current basic reproduction number as ${\mathcal{R}}_0\approx 1.09938>1$.

It is clear from Section 3 that ${\mathcal{R}}_0$ is a parameter that serves as a threshold to determine whether the disease will remain present in the population (see Theorems 3 and 4). According to Theorem 4, if ${\mathcal{R}}_0>1$, system (1) has a positive $\omega$-periodic solution. Figure 3 demonstrates the long-term behavior and uniform persistence of measles when ${\mathcal{R}}_0\approx 1.09938>1$, the disease compartments persist, and the epidemic becomes endemic in the population and recurs annually.

Figure 4 shows that the numerical solution of our model (1) is consistent with the analytic results, which state that the unique disease-free periodic solution $E_0$ is globally asymptotically stable when ${\mathcal{R}}_0=0.544246,0.70752,0.81637<1$.

4.2. The Effects of Different Parameter Changes on the Spread of Measles

In this subsection, a series of numerical simulations are performed to examine the effects of various biological parameters on model (1). Measles is a highly contagious disease and the rate at which it spreads is determined by several factors, including the number of susceptible individuals, the length of the infectious period, and the frequency of contact between infected and susceptible individuals. Understanding the dynamics of transmission is crucial to controlling outbreaks and preventing the further spread of measles. To determine the transmission rate of measles, our mathematical model (1) is used, which takes into account these and other factors. When the transmission rate is high, the spread of the disease is rapid, but when it is low, the spread is slower, as shown in Figure 5a. Figure 6a also shows that a higher contact rate results in an increased basic reproduction number, ${\mathcal{R}}_0$.

During a measles outbreak, understanding the incubation period is crucial. This period, which is the time between being exposed to the virus and showing symptoms, can affect when an infected individual becomes contagious and how long their illness lasts. The average incubation period for measles is 14 days but can range from 7 to 21 days. During this time, the virus spreads throughout the body, but the person does not yet exhibit any symptoms. This means that a person can unknowingly spread the virus to others. To prevent the spread of measles, it is important to be aware of the incubation period and take measures such as staying at home and avoiding contact with others [1,2]. In this work, we explore the effect of the incubation period on the transmission of measles; our findings suggest that a shorter incubation period greatly increases the speed of disease spread (see Figure 5b and Figure 6b).

As seen in Figure 7a, when the vaccination coverage rate ($\theta$) improves, the impact of the measles epidemic decreases. This can be a challenge to achieve in countries such as Pakistan. The government could focus on improving the effectiveness of the vaccine ($\alpha$), as shown in Figure 7b, while also making sure that the vaccine is accessible to those in poverty. Finally, a low coverage rate ($\theta$) for the vaccine, as seen in Figure 8a, is not as problematic as a low efficacy rate ($\alpha$). A higher efficacy rate can help bring ${\mathcal{R}}_0$ closer to under 1, highlighting its importance in mitigating the burden of the measles epidemic (see Figure 8b).

Although increasing the recovery rate ($\gamma$), as shown in Figure 9a, can also reduce the impact of the epidemic, it is not as impactful as reducing the spread of the disease from infected to susceptible individuals. The significance of the second dose of the measles vaccine cannot be overstated. Measles is a highly infectious disease that can have severe consequences, particularly for children and individuals with weakened immune systems. While the first dose of the vaccine provides initial protection (see Figure 7a), it is the second dose that reinforces and maintains this immunity over time, as shown in Figure 9b. Our findings suggest that the overall effectiveness of the vaccine in preventing the spread of measles increases significantly with the administration of a second dose. Thus, it is crucial for individuals to receive the second dose of the measles vaccine to remain protected against this dangerous disease.

In order to determine how changes in certain parameters affect the basic reproduction number (${\mathcal{R}}_0$), we conducted a sensitivity analysis using the Latin hypercube sampling (LHS) method and calculated partial rank correlation coefficients (PRCCs) for various input parameters. Our findings, shown in Figure 10, indicate that the coverage rate $\theta$, vaccine efficiency rate $\alpha$, baseline transmission rate $\beta$, and recovery rate $\gamma$ are the most influential factors in determining ${\mathcal{R}}_0$, as opposed to parameters related to asymptomatic infection. This suggests that taking measures to reduce transmission rates, such as using personal protective equipment and practicing social distancing, promptly treating infected individuals, and increasing vaccination rates, can significantly reduce the number of new infections.

4.3. Effects of Possible Control Measures

We conducted simulations to examine the potential impacts of changing certain parameters. Specifically, we focused on the measles baseline contact rate and vaccine coverage rate, as these are the most likely to be adjusted through intervention efforts. Our objective was to determine to what extent these parameters would need to be modified to prevent future measles outbreaks. The simulations were initiated using the parameters that were fitted up until month 36, and a change was then implemented. The most effective way to reduce measles transmission is to increase the vaccination coverage rate $\theta$ among the population, as shown in Figure 11a. Developing countries such as Pakistan face difficulties in applying the measles vaccine, including limited access to health services, logistics and supply chain challenges, cultural and religious barriers, lack of awareness and education, political instability, and limited financial resources [47]. Isolation is one measure that can be applied to decrease human interaction and slow the spread of the measles virus. However, it is important to consider the possible impacts of isolation on individuals and communities, such as loss of income and social support, before implementing this measure [1,2]. From the simulations, it was observed that the decrease in the person-to-person transmission rate $\beta$ led to a decrease in the number of infected individuals, as shown in Figure 11b.

To control the spread of measles in Pakistan, a comprehensive approach is needed. This involves increasing access to the MMR vaccine, raising public awareness about vaccinations, improving healthcare access and resources, tracking outbreaks, responding quickly, implementing isolation and quarantine measures, and investing in vaccines and treatment research.

5. Discussion

In this article, we present a mathematical model to study the transmission dynamics of measles. The human population is divided into six compartments: susceptible, vaccinated (first and second dose), exposed, symptomatic infectious, and recovered. Our model (1) builds upon previous models [8,9,10] as well as other models on measles transmission dynamics by incorporating a seasonal transmission parameter. The inclusion of this parameter adds an additional layer of complexity to the model, allowing for a more detailed analysis of the disease transmission dynamics. The solution of the system is guaranteed to be non-negative based on non-negative initial values. The basic reproduction number, ${\mathcal{R}}_0$, is calculated and used to determine the stability of the disease-free equilibrium and the persistence of the disease (see Section 3), by using the general procedure presented in [24]. Our results show that if the basic reproduction number is less than one, the disease will become extinct, whereas if it is greater than one, the disease will persist with a positive periodic solution (see Figure 3 and Figure 4).

In this work, we examined the spread of measles in Pakistan during the outbreak from January 2019 to December 2021 as a case study. To match the seasonal fluctuation data, we used a mathematical model with periodic parameters. We assume that the transmission rate, which determines the rate at which the disease spreads, follows a periodic pattern of 12 months. We used Latin hypercube sampling to determine the best-fit parameters for the model, which were fit to the observed number of measles cases in Pakistan. Our model provides a reasonably good fit to the data, with a high goodness of fit ($R^2=0.749996$), as shown in Figure 2. The estimated basic reproduction number of the disease, ${\mathcal{R}}_0$, is $1.09938$, which is greater than 1, which means that the disease will remain present in the population. Our model also demonstrates the long-term behavior of the disease, which becomes endemic in the population and recurs annually.

The spread of measles is determined by various biological factors, including the number of susceptible individuals, the length of the infectious period, and the frequency of contact between infected and susceptible individuals. The transmission rate of measles can be modeled using our mathematical model, which takes into account these and other factors. The results of simulations conducted to examine the effect of various biological parameters suggest that a shorter incubation period greatly increases the speed of disease spread (see Figure 5b), while a higher vaccination coverage rate decreases the impact of the epidemic (see Figure 7a). A higher rate of vaccine efficacy can also help bring the basic reproduction number (${\mathcal{R}}_0$) below one (see Figure 8b), highlighting its importance in mitigating the spread of measles. The significance of the second dose of the measles vaccine is also emphasized (see Figure 9b).

This improvement over previous models, such as [8,9,10,18,20,48], provides a better understanding of the impact of seasonal changes on the spread of the disease. The results of our study provide valuable information to develop and implement effective measles control strategies. The most effective strategy to reduce measles transmission and prevent future outbreaks is to increase vaccination coverage among the population. By being vaccinated, individuals can protect themselves and others from becoming sick, helping to break the transmission chain. This can be achieved through various means, such as public health campaigns, school-based vaccination programs, and targeted outreach efforts to communities with low vaccination rates. The ultimate goal is to achieve herd immunity, where a high enough proportion of the population is immune to the virus so that it is unable to spread and outbreaks are prevented. Other measures, such as isolation, quarantine, and decreasing the person-to-person transmission rate, can also help control the spread of measles, but a comprehensive approach involving various factors, such as improving access to the MMR vaccine, raising public awareness, and investing in vaccine and treatment research, is needed. It is important to implement a combination of these methods for maximum effect.