Vaccination Thresholds and Bifurcation Analysis in a Measles SIR Model: A Comparative Numerical Study Using Variational Iteration Method and Runge-Kutta Method

Measles is a highly contagious viral disease that poses a significant public health challenge, particularly in regions with low vaccination coverage. Despite the availability of an effective vaccine, measles continues to cause outbreaks worldwide, especially in areas where immunization efforts are insufficient or inconsistent. Understanding the transmission dynamics of measles is crucial for designing effective control strategies and informing public health policies.

Agusto [1] explored the mathematical model for ebola virus and its stimulations. Althaus et al. [2] estimated the basic reproduction number for the Ebola virus through transmission rate. Chowell et al. [3] offered a comprehensive overview of mathematical models applied to Ebola’s transmission and containment. D’Silva and Eisenberg [4] used spatial invasion framework, incorporating local and cross regional transmission using gravity-based model approach. Li [5] developed a model for Ebola virus in Sierra Leone, emphasizing the influence of human movement on viral dissemination. Njankou and Nyabadza [6] constructed an optimal model control to access the awareness campaigns and medical interventions. Rachab and Torres [7] adopted simulation driven optimal control models to evaluate vaccine deployment strategies during the 2014 epidemic.

Dia et al. [8] introduced a modified Susceptible- Exposed-Infected-Recovered frame work to detailed the epidemic behaviour in Liberia. Meanwhile Sedelnikov [9] explored the fuzzy logic in disciplines like space science and engineering. Coltart et al. [10] traced the trajectory of 2013-2016 west African Ebola, outbreak highlighting the importance of data centric decision making and control measures. Barros et al. [11] studied the use of fuzzy dynamic system in epidemic modelling, while Farahi and Barati [12] applied fuzzy methodologies in time delayed system in epidemiology. Das and Pal [13] focused on Susceptible–Infected–Recovered (SIR) models with treatment strategies under fuzzy condition. Verma et al. [14] introduced fuzzy based model to analyze influence transmission and control. Sadhukhan et al. [15] explored optimal harvesting in a trophic model within a fuzzy logic framework. Bhavithra et al. [16, 17] examined parameters like basic reproduction number, backward, forward, transcritical bifurcation analysis as well as sensitivity analysis for SEIR model using sinusoidal function.

Sedelnikov et al. [18] also addressed interpretation challenges in telemetry data from the Advanced Industrial Science and Technology small satellite. Ahmad et al. [19] designed a fractional order COVID 19 model using Caputo derivatives. Yadav et al. [20] reviewed common epidemic model such as SI, SIR, SEIR analyzing key parameters particularly basic reproduction number forecasting in COVID 19 trends. Gurjar et al. [21] investigated the SIR model for COVID 19 in Indonesia, emphasising how face to face learning and vaccination affected transmission. Dwivedi et al. [22] constructed a nonlinear deterministic model to analyse the co-circulation of Japanese encephalitis and dengue. Bifurcation analysis for novel fuzzy SIR model has been studied by Bhavithra et al. [23].

Mathematical modeling plays a vital role in the study of infectious diseases by helping to predict disease spread, assess intervention strategies, and estimate key epidemiological thresholds. Among various modeling approaches, the SIR model is widely used for analyzing diseases like measles due to its simplicity and effectiveness in capturing core transmission dynamics. The basic reproduction number, often denoted as R-naught, serves as a critical threshold parameter that determines whether an infectious disease can invade and persist in a population. If the value is greater than one, the disease can spread; if it is less than one, the infection will eventually die out.

Several studies have extended the basic SIR model by including factors such as vaccination, birth and death rates, time delays, and age structure. Vaccination, in particular, has been incorporated into many models to determine the critical vaccination coverage needed to achieve herd immunity and eradicate the disease. Previous works have also employed various numerical methods to solve these models, including Runge-Kutta methods, finite difference schemes, and spectral collocation techniques. However, there is a growing need for efficient and accurate numerical approaches that can handle the nonlinear nature of these models. In this study, we develop a measles transmission model based on the SIR framework with the inclusion of a vaccination term. We analyze the system to determine the equilibrium points and their stability, and we perform a bifurcation analysis to understand the transition between disease-free and endemic states. We also apply the Variational Iteration Method and Runge Kutta method to obtain accurate numerical approximations of the system's behavior. The results provide insights into the role of vaccination in controlling measles and demonstrate the effectiveness of the proposed numerical approach.

The critical vaccination threshold, or the critical vaccination proportion, is a concept in epidemiology that indicates the fraction of a population that needs to be vaccinated to prevent an infectious disease from spreading within that population. This threshold is crucial for achieving herd immunity, where indirect protection from infectious disease is provided to unvaccinated individuals when a significant portion of the population is vaccinated.

The formula to find the critical vaccination proportion $\left(P_c\right)$ is derived from the reproduction number $\left(R_0\right)$, which is the average number of cases one infected individual will generate over the course of their infectious period in a completely susceptible population. The formula is: $P_c=1-\frac{1}{R_0}=1- \frac{\gamma+\pi}{\beta}$. If $P_c>p$, measles free equilibrium is locally stable with the coordinates $\mathrm{M}_0=(1-\mathrm{p}, 0,0)$.

Theorem 1:

Let S denote the susceptible population, I the infected population, and R the recovered population. Suppose the susceptible population S is four times the infection rate I, and the recovered population R is equal to the infection rate I. Let $\pi$ represent the fertility rate, γ denote the recovery rate and p denote the vaccination coverage for newborns. The following relationships hold for the rate of recovery γ as a function of the fertility rate π and the vaccination coverage p:

1. When $p=\mathbf{1}$ (newborn vaccination is 100%): then $\gamma \propto -5 \pi$.
2. When $p=0.5$ (newborn vaccination is 50%): then $\gamma \propto =2 \pi$.
3. When $p=1 / 3$ (newborn vaccination is approximately 33.33%): then $\gamma \propto-\pi$.

Here, $\propto$ denotes proportionality. These relationships describe how the rate of recovery $\boldsymbol{\gamma}$ varies with changes in the vaccination coverage $\mathbf{p}$ relative to the fertility rate $\boldsymbol{\pi}$.

Proof:

From our model we have,

$\frac{d S}{d t}=(1-\mathrm{p}) N(t) \pi-\beta \mathrm{SI}-\pi \mathrm{S}$           (1)

$\frac{d I}{d t}=\beta \mathrm{SI}-(\gamma+\pi) I$              (2)

$\frac{d R}{d t}=N(t) p \pi+\gamma I-\pi R$                 (3)

Dividing Eq. (1) by Eq. (2), we get

$\frac{d S(t)}{d I(t)}=\frac{(1-\mathrm{p}) \mathrm{N}(\mathrm{t}) \pi-\beta \mathrm{SI}-\pi \mathrm{S}}{\beta \mathrm{SI}-(\gamma+\pi) I}$

$\frac{4 d I(t)}{d I(t)}=\frac{(1-\mathrm{p}) \pi 6 I(t)-4 \beta I(t) I(t)-4 \pi I(t)}{4 \beta I(t) I(t)-(\gamma+\pi) I(t)}$

$4=\frac{(1-p) \pi 6-4 \beta I(t)-4 \pi}{4 \beta I(t)-(\gamma+\pi)}$

$\begin{gathered}4[4 \beta I(t)-(\gamma+\pi)]=(1-\mathrm{p}) \pi 6-4 \beta I(t)-4 \pi \\ 16 \beta I(t)-4(\gamma+\pi)=(1-\mathrm{p}) \pi 6-4 \beta I(t)-4 \pi \\ 20 \beta I(t)=6(1-\mathrm{p}) \pi-4 \pi+4(\gamma+\pi)\end{gathered}$

$I(t)=\frac{6(1-p) \pi-4 \pi+4(\gamma+\pi)}{20 \beta}$                    (4)

Dividing Eq. (1) by Eq. (3), we get

$\begin{gathered}\frac{d S(t)}{d R(t)}=\frac{(1-\mathrm{p}) \mathrm{N}(\mathrm{t}) \pi-\beta \mathrm{SI}-\pi \mathrm{S}}{N(t) p \pi+\gamma I-\pi R} \\ \frac{4 d I(t)}{d I(t)}=\frac{(1-\mathrm{p}) \pi 6 I(t)-4 \beta I(t) I(t)-4 \pi I(t)}{6 I(t) p \pi+\gamma I(t)-\pi I(t)} \\ 4=\frac{(1-\mathrm{p}) \pi 6-4 \beta I(t)-4 \pi}{6 p \pi+\gamma-\pi} \\ 4[6 p \pi+\gamma-\pi]=(1-\mathrm{p}) 6 \pi-4 \beta \mathrm{I}(\mathrm{t})-4 \pi \\ 4 \beta \mathrm{I}(\mathrm{t})=(1-\mathrm{p}) 6 \pi-4 \pi-4[6 p \pi+\gamma-\pi] \\ I(t)=\frac{6 \pi-4 \gamma-30 p \pi}{4 \beta}\end{gathered}$         (5)

Now Eqs. (4) and (5) we get,

$\frac{6(1-\mathrm{p}) \pi-4 \pi+4(\gamma+\pi)}{20 \beta}=\frac{6 \pi-4 \gamma-30 p \pi}{4 \beta}$

$6(1-p) \pi-4 \pi+4(\gamma+\pi)=5[6 \pi-4 \gamma-30 p \pi]$

$24 \gamma=24 \pi-144 p \pi, 24 \gamma=(24-144 p) \pi$

Now when $p=1$ (newborn vaccination is 100%): $\gamma= \frac{(24-144) \pi}{24} \propto-5 \pi$.

When $p=0.5$ (newborn vaccination is 50%): then $\gamma= \frac{(24-72) \pi}{24} \propto-2 \pi$.

When $p=1 / 3$ (newborn vaccination is approximately 33.33%): then $\gamma=\frac{(24-48) \pi}{24} \propto-\pi$.

Hence proved.

Theorem 2:

Let S denote the susceptible population, I the infected population, and R the recovered population. Suppose the susceptible population S is twice the infection rate I, and the recovered population R is equal to the infection rate I. Let π represent the fertility rate, γ denote the recovery rate and p denotes the vaccination coverage for newborns. The following relationships hold for the rate of recovery γ as a function of the fertility rate π and the vaccination coverage p:

1. When $p=0$ (newborn vaccination is 0%): then $\gamma \propto \pi$.
2. When $p=1$ (newborn vaccination is 100%): then $\gamma \propto -3 \pi$.
3. When $p=1 / 6$ (newborn vaccination is approximately 16.66%): then $\gamma \propto \pi / 3$.
4. When $p=1 / 4 \mathrm{~m}, \mathrm{~m}>1$: then $\gamma \propto \frac{\mathrm{m}-1}{\mathrm{~m}} \pi$.

Here, $\propto$ denotes proportionality. These relationships describe how the rate of recovery $\gamma$ varies with changes in the vaccination coverage $p$ relative to the fertility rate $\pi$.

Proof:

From our model we have,

$\frac{d S}{d t}=(1-\mathrm{p}) N(t) \pi-\beta \mathrm{SI}-\pi \mathrm{S}$                (6)

$\frac{d I}{d t}=\beta S I-(\gamma+\pi) I$                    (7)

$\frac{d R}{d t}=N(t) p \pi+\gamma I-\pi R$                 (8)

Dividing Eq. (6) by Eq. (7), we get

$\begin{gathered}\frac{d S(t)}{d I(t)}=\frac{(1-\mathrm{p}) \mathrm{N}(\mathrm{t}) \pi-\beta \mathrm{SI}-\pi \mathrm{S}}{\beta \mathrm{SI}-(\gamma+\pi) I} \\ \frac{2 d I(t)}{d I(t)}=\frac{(1-\mathrm{p}) \pi 4 I(t)-2 \beta I(t) I(t)-2 \pi I(t)}{2 \beta I(t) I(t)-(\gamma+\pi) I} \\ 2=\frac{4(1-\mathrm{p}) \pi-2 \beta I(t)-2 \pi}{2 \beta I(t)-(\gamma+\pi)} \\ 2[2 \beta I(t)-(\gamma+\pi)]=4(1-\mathrm{p}) \pi-2 \beta I(t)-2 \pi \\ 6 \beta I(t)=4 \pi-4 \pi p+2 \gamma \\ I(t)=\frac{4 \pi-4 \pi p+2 \gamma}{6 \beta}\end{gathered}$               (9)

Dividing Eq. (6) by Eq. (8), we get

$\begin{gathered}\frac{d S(t)}{d R(t)}=\frac{(1-\mathrm{p}) \mathrm{N}(\mathrm{t}) \pi-\beta \mathrm{SI}-\pi \mathrm{S}}{N(t) p \pi+\gamma I-\pi R} \\ \frac{2 d I(t)}{d I(t)}=\frac{(1-\mathrm{p}) \pi 4 I(t)-2 \beta I(t) I(t)-2 \pi I(t)}{4 I(t) p \pi+\gamma I(t)-\pi I(t)} \\ 2=\frac{(1-\mathrm{p}) \pi 4-2 \beta I(t)-2 \pi}{4 p \pi+\gamma-\pi} \\ 2[4 p \pi+\gamma-\pi]=(1-p) 4 \pi-2 \beta I(t)-2 \pi \\ 8 p \pi+2 \gamma-2 \pi=(1-p) 4 \pi-2 \beta I(t)-2 \pi \\ 2 \beta I(t)=4 \pi-2 \gamma-12 p \pi \\ I(t)=\frac{4 \pi-2 \gamma-12 p \pi}{2 \beta}\end{gathered}$                 (10)

Now equating Eq. (6) and Eq. (10) we get,

$\begin{gathered}\frac{4 \pi-4 \pi p+2 \gamma}{6 \beta}=\frac{4 \pi-2 \gamma-12 p \pi}{2 \beta} \\ 4 \pi-4 \pi p+2 \gamma=3[4 \pi-2 \gamma-12 p \pi] \\ 8 \gamma=(1-4 p) 8 \pi ; \gamma=(1-4 p) \pi\end{gathered}$

Now, when $p=0$ (newborn vaccination is 0%): then $\gamma= \pi$.

When $p=1$ (newborn vaccination is 100%): then $\gamma= (1-4) \pi=-3 \pi$.

When $p=1 / 6$ (newborn vaccination is approximately 16.66%): then $\gamma=\left(1-\frac{2}{3}\right) \pi=\pi / 3$.

When $p=1 / 4 m, m>1$: then $\gamma=\left(1-\frac{4}{4 m}\right) \pi=\frac{m-1}{m} \pi$

Hence proved.

In SIR mathematical model, the bifurcation parameter is often related to basic reproduction number, which is a key determinant of disease transmission dynamics.it is a dimensionless number that integrates various biological environment and behavioural factors affecting the spread of infectious disease. The bifurcation parameter of our model is found to be $\beta^*=\frac{\gamma+\pi}{1-p}$.

5.1 Transcritical bifurcation

In the context of the basic SIR model used in epidemiology to describe the spread of infectious diseases, a transcritical bifurcation refers to a type of critical point where the stability of the disease-free equilibrium (DFE) and the endemic equilibrium swap as a parameter crosses a critical value.

Let us consider the set of equations Let us consider the set of equations $\frac{d S}{d t}=(1-\mathrm{p}) \pi-\beta \mathrm{SI}- \pi S$:

$\frac{d I}{d t}=\beta S I-(\gamma+\pi) I ; \frac{d R}{d t}=p \pi+\gamma I-\pi R$

The Jacobian of the above equation is found as

$J=\left[\begin{array}{ccc}-\beta \mathrm{I}-\pi & -\beta \mathrm{S} & 0 \\ \beta \mathrm{I} & \beta \mathrm{S}-(\gamma+\pi) & 0 \\ 0 & \gamma & -\pi\end{array}\right]$

Now we find the Jacobian at $\left(M_0, \beta^*\right)$

$J_{\left(M_0, \beta^*\right)}=\left[\begin{array}{ccc}-\pi & -\beta^*(1-\mathrm{p}) & 0 \\ 0 & 0 & 0 \\ 0 & \gamma & -\pi\end{array}\right]$

The characteristic equation is given by $|J-\lambda I|=0$

$\left|\begin{array}{ccc}-\pi-\lambda & -\beta^*(1-\mathrm{p}) & 0 \\ 0 & -\lambda & 0 \\ 0 & \gamma & -\pi-\lambda\end{array}\right|=0$

$(-\pi-\lambda)(-\lambda)(-\pi-\lambda)=0$

The eigen values are $0,-0.2,-0.2$.

Since it has one zero eigen value and 2 negative eigen value, hence our model does not have saddle point.

We can proceed further, let $\lambda_1=0, \lambda_2=\lambda_3=-0.2$, now we find the eigen vector $X=\left(x_1, x_2, x_3\right)^T$ corresponding to $\lambda_1$ which satisfy $J X=\lambda_1 X$.

$\left[\begin{array}{ccc}-\pi & -\beta^*(1-p) & 0 \\ 0 & 0 & 0 \\ 0 & \gamma & -\pi\end{array}\right]\left[\begin{array}{l}x_1 \\ x_2 \\ x_3\end{array}\right]=0$

$\begin{gathered}-\pi x_1-\beta^*(1-\mathrm{p}) x_2=0 \\ \gamma x_2-\pi x_3=0\end{gathered}$

Solving the above equations, we get $\pi x_1=-\beta^*(1-\mathrm{p}) x_2$.

Thus $x_1=\frac{-\beta^*(1-\mathrm{p})}{\pi} x_2$.

Also, $x_3=\frac{\gamma}{\pi} x_2$.

Let $A_1=\frac{-\beta^*(1-\mathrm{p})}{\pi}$ and $A_2=\frac{\gamma}{\pi}$.

Hence $x_1=A_1 x_2$ and $x_3=A_2 x_2$.

Thus $X=\left(\begin{array}{c}A_1 x_2 \\ x_2 \\ A_2 x_2\end{array}\right)$.

Now we find the eigen vector $Y=\left(y_1, y_2, y_3\right)^T$. corresponding to $\lambda_1$ which satisfy $J^T Y=\lambda_1 Y$.

$\left[\begin{array}{ccc}-\pi & 0 & 0 \\ -\beta^*(1-\mathrm{p}) & 0 & \gamma \\ 0 & 0 & -\pi\end{array}\right]\left[\begin{array}{l}y_1 \\ y_2 \\ y_3\end{array}\right]=0$

$-\pi y_1=0,-\beta^*(1-\mathrm{p}) y_1+\gamma y_3=0,-\pi y_3=0$

Solving the above equations, we get

$y_1=y_3=0 ; Y=\left(\begin{array}{c}0 \\ y_2 \\ 0\end{array}\right)$

Now we write system of equations in equations in the vector form as $\frac{d Y}{d t}=g(x)$, where

$Y=(S, I, R)^T ; g=\left(g_1, g_2, g_3\right)^T$

Now to calculate $\frac{d g}{d \beta}=g_\beta$.

$\begin{aligned} & \text { Let } g_1=(1-p) \pi-\beta S I-\pi S ; g_2=\beta S I-(\gamma+\pi) I ; g_3 = p \pi+\gamma I-\pi R ; \frac{d g_1}{d \beta}=-S I ; \frac{d g_2}{d \beta}=-S I ; \frac{d g_3}{d \beta}=0 ;g_\beta\left(M_0, \beta^*\right)=\left(\begin{array}{l}0 \\ 0 \\ 0\end{array}\right) .\end{aligned}$

Second condition to proceed Sotomayor theorem is also verified.

Now to check $Y^T\left[D g_\beta\left(M_0, \beta^*\right)\right] X \neq 0$,

$D g_\beta\left(M_0, \beta^*\right)=\left[\begin{array}{ccc}0 & -(1-\mathrm{p}) & 0 \\ 0 & 1-\mathrm{p} & 0 \\ 0 & 0 & 0\end{array}\right]$

Now,

$\begin{aligned} Y^T\left[D g_\beta\left(M_0, \beta^*\right)\right] X & =\left(y_1, y_2, y_3\right)\left[\begin{array}{ccc}0 & -(1-\mathrm{p}) & 0 \\ 0 & 1-\mathrm{p} & 0 \\ 0 & 0 & 0\end{array}\right]\left[\begin{array}{l}x_1 \\ x_2 \\ x_3\end{array}\right] = (1-p) x_2 y_2 \neq 0\end{aligned}$

According to Sotomayor theorem, when the parameter $\beta$ passing through $\beta^*$ then transcritical bifurcation occurs at measles free equilibrium point.

The comparison between the RK and VIM for the SIR model reveals consistent differences in the estimated dynamics of the susceptible and recovered populations over a 30-day period are shown in Figures 2 and 3. With identical initial conditions (S(0) = 1, I(0) = 0, R(0) = 0), the RK method yields a faster decline in the susceptible population and a slower increase in the recovered population compared to VIM. Throughout the simulation, VIM consistently predicts higher values for the susceptible class and more rapid accumulation in the recovered class, indicating a relatively accelerated disease resolution.

22.png

222.png

303.png

Figure 3. Comparison of basic reproduction number by VIM and RK

Table 2. Comparison of susceptible case by VIM and RK method

Table 3. Comparison of recovered case by VIM and RK method

Table 4. Comparison of basic reproduction number by VIM and RK

Additionally, the basic reproduction number (R₀) obtained through both methods demonstrates a decreasing trend over time; however, values obtained via VIM are marginally higher than those from RK at each step, suggesting a slightly slower perceived reduction in transmission under VIM. These discrepancies highlight the influence of numerical technique on the simulation outcomes and suggest that VIM may offer a more conservative estimate of disease spread dynamics. Both methods effectively replicate the qualitative behavior of the epidemic, yet the variation in their outputs underscores the importance of method selection in epidemiological modeling. Simulations further show that increasing vaccination coverage significantly alters the trajectory of measles transmission. At low coverage levels, the model predicts large outbreaks, while increasing vaccination rates toward the critical threshold leads to a sharp decline in infection levels. At the critical vaccination threshold, the system exhibits a transcritical bifurcation, where the disease-free equilibrium becomes stable and the endemic equilibrium vanishes. For a basic reproduction number of 10, the estimated critical vaccination coverage is 90%, and simulations near this threshold show that even slight deviations (89.9%) can result in minor outbreaks. This demonstrates the importance of maintaining high vaccination coverage to eliminate disease persistence and prevent resurgence, emphasizing the sensitivity of measles control to precise vaccination rates (Tables 2-4).