Vaccination strategies in a stochastic $\mathscr{SIVR}$ epidemic model

1.   Introduction

Infectious diseases have been a persistent and serious threat to human health, often resulting in pandemics and outbreaks that profoundly impact society, economics, and politics. A comprehensive understanding of disease dynamics is essential for public health authorities and researchers worldwide to develop effective control strategies. Mathematical models, such as the susceptible-infected-vaccinated-recovered $\mathscr{SIVR}$ model, play a vital role in studying disease transmission dynamics ^[1]. These models provide valuable insights into how diseases spread through populations and how interventions, such as vaccination campaigns or public health measures, can mitigate their impact.

In recent years, various modifications to the $\mathscr{SIVR}$ model have been proposed to better reflect real-world epidemiological complexities, including time-dependent vaccination rates, waning immunity, and heterogeneous population structures ^[2,3]. In addition, tochastic epidemic models have gained attention due to their ability to capture randomness in disease transmission ^[4]. Unlike deterministic models, stochastic models provide a more realistic framework by accounting for unpredictable fluctuations in infection rates, intervention effectiveness, and environmental factors ^[5]. However, many existing models still assume deterministic vaccination strategies, overlooking the impact of randomness on vaccine availability, public compliance, and external constraints.

To address this limitation, we incorporate stochastic perturbations into the vaccination term of the $\mathscr{SIVR}$ model. Specifically, we incorporate white noise to model random fluctuations in vaccination rates. This approach aligns with recent advanceshastic epidemiological modeling, where noise-driven fluctuations significantly influence disease outcomes ^[6,7].

The stochastic differential equation (SDE) model introduces random fluctuations in epidemic dynamics, enhancing realism compared to the deterministic ordinary differential equation (ODE) model. Unlike the ODE model, the SDE model incorporates stochastic noise, representing real-world uncertainties in disease transmission and recovery rates. Even if the deterministic model suggests the persistence of the disease ($R_0 > 1$), stochastic effects can drive the infection to extinction by causing fluctuations that bring the number of infected individuals to zero. The SDE model also exhibits mean-reverting behavior, in contrast to the ODE model, which predicts a stable equilibrium.

Our model extends existing stochastic epidemic frameworks by explicitly incorporating random fluctuations in vaccination and transmission rates. Unlike prior deterministic models ^[8,9], which assume constant vaccination rates, our stochastic approach accounts for variations due to logistical constraints and behavioral responses. This allows for a more accurate prediction of disease dynamics under uncertainty ^[10,11].

The proposed stochastic $\mathscr{SIVR}$ model introduces several key innovations and contributions. Incorporating white noise into vaccination dynamics to enhance real-world applicability. In addition, it derives a refined reproduction number $R^*_0$ that accounts for stochastic effects. A Lyapunov-based proof is provided to establish global stability and persistence conditions. Furthermore, numerical simulations are performed to illustrate the role of randomness in disease extinction ^[12].

The stochastic $\mathscr{SIVR}$ model was developed to extend the classical $\mathscr{SIR}$ model in several ways. Unlike the $\mathscr{SIR}$ model, the $\mathscr{SIVR}$ model explicitly accounts for vaccinated individuals. Stochastic perturbations can reduce the reproduction number $R^*_0$ below 1, facilitating earlier disease eradication. While the deterministic $\mathscr{SIR}$ model predicts smooth transitions, the stochastic $\mathscr{SIVR}$ model introduces fluctuations, making it more applicable to real outbreaks. Additionally, the stochastic $\mathscr{SIVR}$ model captures irregularities in vaccination rates caused by logistical challenges, which are not considered in the $\mathscr{SIR}$ model ^[12].

By addressing critical gaps in existing literature, our work advances stochastic epidemiology by providing a robust framework for evaluating public health interventions under uncertainty. In this study, we begin with the deterministic $\mathscr{SIVR}$ model ^[13] and systematically extend it into a stochastic framework by introducing perturbations in key parameters using stochastic noise. This approach allows us to realistically capture the inherent randomness in vaccination rates and disease transmission dynamics.

2.   Stochastic model formulation

The human host population, denoted by $\mathscr{N}$, is structured into four dynamic compartments: susceptible ($\mathscr{S}$), infected ($\mathscr{I}$), vaccinated ($\mathscr{V}$), and recovered ($\mathscr{R}$), with the total population $\mathscr{N}$ evolving over time $t$ as:

The below deterministic model (2.1) is described by the following parameters ^[13]:

● $A$: Recruitment rate of susceptible individuals.

● $\beta$: Transmission rate of infections.

● $\rho$: Fraction of susceptible individuals vaccinated.

● $(1 - \rho)$: Fraction of susceptible individuals not vaccinated.

● $\delta_1$: Rate of waning vaccination in vaccinated individuals.

● $\mu$: Natural death rate.

● $\alpha$: Disease-induced death rate in infected individuals.

● $\gamma$: Recovery rate of infected individuals.

In the deterministic model (2.1), the basic reproduction number, denoted as $R_0$, represents the average number of secondary infections generated by a single infected individual in a fully susceptible population. This metric plays a crucial role in evaluating the potential for an outbreak and devising effective control measures.

The expression for $R_0$ in the context of the $\mathscr{SIVR}$ model is given by ^[13]:

To account for real-world uncertainties, we extend the deterministic model by incorporating stochastic perturbations in the transmission and vaccination processes. These perturbations are modeled using white noise, which captures fluctuations in vaccine availability, public compliance, and environmental factors.

The stochastic $\mathscr{SIVR}$ model is formulated as follows:

Here, $\sigma_1$ and $\sigma_2$ represent the intensity of stochastic perturbations affecting the transmission rate $\beta\rightarrow \beta+\sigma_1\mathrm{d}\mathscr{B}_1(t)$ and recovery rate $\gamma\rightarrow \gamma+\sigma_2\mathrm{d}\mathscr{B}_2(t)$, respectively. The terms $\mathrm{d}\mathscr{B}_1(t)$ and $\mathrm{d}\mathscr{B}_2(t)$ denote independent standard Brownian motions.

By incorporating these stochastic perturbations, the model provides a more realistic representation of disease dynamics, capturing unpredictable variations that occur in real-world epidemiological settings. This stochastic approach enables better assessment of disease control strategies under uncertainty.

3.   Existence and uniqueness

Theorem 3.1. Let the initial condition be $\zeta_{(0)} = (\mathscr{S}_{(0)}, \mathscr{I}_{(0)}, \mathscr{V}_{(0)}, \mathscr{R}_{(0)})$, where $\mathscr{S}_{(0)}$, $\mathscr{I}_{(0)}$, $\mathscr{V}_{(0)}$, and $\mathscr{R}_{(0)}$ are non-negative. Then, the stochastic system given by (2.2) has a unique, global, non-negative solution $(\mathscr{S}_{(t)}, \mathscr{I}_{(t)}, \mathscr{V}_{(t)}, \mathscr{R}_{(t)})$ for all $t \geq 0$. Moreover, this solution remains in the non-negative quadrant $\mathbb{R}_+^4$ with probability 1 almost surely ^[14,15,16].

Proof. Since the coefficients of the stochastic differential equations in (2.2) satisfy the Lipschitz condition, for any given initial value $\zeta_0 = \big(\mathscr{S}_0, \mathscr{I}_0, \mathscr{V}_0, \mathscr{R}_0 \big) \in \mathbb{R}_+^3$, there exists a unique local solution $(\mathscr{S}_{(t)}, \mathscr{I}_{(t)}, \mathscr{V}_{(t)}, \mathscr{R}_{(t)})$ on $t \in [0, \tau_e)$, where $\tau_e$ is the explosion time.

To ensure global stability, we must show that $\tau_e = \infty$. Let $a_0 \geq 0$ be sufficiently large so that $\zeta_{(0)}$ lies within the interval $\left[\frac{1}{a_0}, a_0\right]$.

For each integer $a \geq a_0$, we define the following stopping time:

Now, the total number of individuals in the system satisfies the inequality:

Solving the above equation, we obtain:

Thus, the total population remains constrained over time.

Now, suppose that

where $\varphi$ denotes the empty set, as usual. By the definition of stopping time, $\tau_a$ is increasing as $a \to \infty$. Suppose that

Now, we must show that $\tau_{\infty} = \infty$. Suppose that there exist a pair of constants $T > 0$ and $\epsilon \in (0, 1)$ such that:

Then, there exists an integer $a_1 \geq a_0$ such that:

Finally, we define a $\mathcal{C}^2$-function $\mathcal{V}: \mathbb{R}_+^4 \to \mathbb{R}_+$, by

Since ${\mathrm{log}}\; \hbar\leq\hbar-1$, for all $\hbar > 0$, $\mathcal{V}$ is clearly non-negative. Consider $a\; \geq\; a_0, \; T > 0$, and use Ito's formula on (3.8), we obtain:

where $\mathscr{L}\; d\mathcal{V}(\mathscr{S}+\mathscr{I}+\mathscr{V}+\mathscr{R}):\mathbb{R}_+^4\rightarrow \mathbb{R}_+$ is

Hence consider the following inequalities:

For all $a\geq a_1$, let $\tau_{a}\leq\; T = \Omega_{a}$ then $P(\Omega_{a})\geq\epsilon$, by corresponding to each $\omega\in\Omega_{a}$, we can find $\mathscr{S}(\tau_{a}, \omega), \mathscr{I}(\tau_{a}, \omega), \mathscr{V}(\tau_{a}, \omega), \mathscr{R}(\tau_{a}, \omega)$ whose value is either $a$ or $\frac{1}{a}$, consequently, $\mathcal{V}(\mathscr{S}(\tau_{a}), \mathscr{I}(\tau_{a}), \mathscr{V}(\tau_{a}), \mathscr{R}(\tau_{a}))$ is not less than ($\mathrm{log}\; a+\frac{1}{a}-1$) or (${a}-1-\mathrm{log}{a}$), so

Taking into account Eqs (3.6) and (3.11), we have

where $\mathbf{1}_{\Omega_{\omega}}$ is the indicator function of ${\Omega}_a$. Suppose $a \to \infty$ leads to a contradiction, then:

Therefore, we must have:

Hence, the assertion follows. □

3.1.   Remark

From Theorem (3.1), it is evident that for any initial condition $\zeta_{(0)} \in \mathbb{R}+^4$, there exists a unique global solution $\zeta_{(t)} \in \mathbb{R}+^4$ almost surely for model (2.2). This solution satisfies the inequality

Solving this inequality yields

If $\frac{\mu}{A-(1-\rho)} \geq \mathscr{N}_{(0)}$, then $\frac{\mu}{A-(1-\rho)} \geq \mathscr{N}_{(t)}$ almost surely. Thus, we define

assuming that $\zeta_{(0)} \in \Lambda^*$ always holds.

4.   Extinction and persistence

To determine the conditions for extinction and persistence in the context of model (2.2), we introduce key concepts that will underpin our main findings. Let $\langle z(t)\rangle = \frac{1}{t}\left[\int_0^t z_{(s)} ds\right]$ denote the time average of $z(t)$, and define the stochastic basic reproduction number as

This modified reproduction number $R^*_0$ incorporates the stochastic effects on disease transmission, adjusting for the variability introduced by the noise term $\sigma_2$ affecting the recovery rate $\gamma$.

4.1.   Definition

The condition for the model (2.2) to be persistent in mean is^[17]

We now prove our main result, showing the condition for which extinction of model (2.2) holds.

Theorem 4.1. Suppose the solution of Eq (2.2) is $\zeta_{(t)} = (\mathscr{S}_{(t)}, \mathscr{I}_{(t)}, \mathscr{V}_{(t)}, \mathscr{R}_{(t)})$ with $\zeta_{(0)} \in \Omega^*$. If $\alpha(\mu+\gamma) > \sigma_2^2\beta$ and $R^*_0 < 1$, then $\lim_{t \to \infty} \left(\frac{\mathrm{log} \mathscr{I}_{(t)}}{t} \right) < 0$. Additionally, $\lim_{t\rightarrow \infty}\; \mathscr{S}_{(t)} = \frac{A-(1-\rho)}{\mu}$, $\lim_{t \to \infty} \mathscr{I}_{(t)} = 0$, $\lim_{t \to \infty} \mathscr{V}_{(t)} = 0$, and $\lim_{t \to \infty} \mathscr{R}_{(t)} = 0$.

Proof. Integrating Eq (2.2), we obtain

Integrating $\mathscr{S}(s)$ over $[0, t]$, we obtain

where

As $\Theta(t)\rightarrow 0$ as $t\rightarrow \infty$, we can consider the second equation of (2.2) and apply Ito's formula to it, yielding

Integrating from 0 to $t$, we have

where

By the law of large numbers ^[18,19,20] $\Phi(t) = 0$ as $t\rightarrow \infty$; therefore;

This implies

From the third equation of model (2.2), it follows that

By applying L'Hopital's Rule to the above expression ^[21], we have

From the fourth equation of model (2.2), it follows that

By applying L'Hopital's Rule to the above expression, we have

From (3.2) it follows that

This concludes the proof. □

5.   Persistence of the disease

In this section, we will investigate conditions for the persistence of the disease.

Theorem 5.1. Assume $\mu > \left(\frac{1}{2}\sigma_1^2\vee\frac{1}{2}\sigma_2^2\right)$. Let $(\mathscr{S}_{(t)}, \mathscr{I}_{(t)}, \mathscr{V}_{(t)}, \mathscr{R}_{(t)})$ be the solution of the stochastic model (2.2), with the initial values given by $(\mathscr{S}_{(0)}, \mathscr{I}_{(0)}, \mathscr{V}_{(0)}, \mathscr{R}_{(0)})\in\mathbb{R}^4_+$. If $R_0^* > 1$, then

Proof. Given that $R_0^* > 1$, according to [22, Lemmas 5.1 and 5.2], we establish the following relation:

Utilizing Eq (4.8), we derive:

The third equation of (2.2) suggests that

Using Eq (5.3) in the above expression yields

Proceeding from the fourth equation of model (2.2), we deduce:

Employing (4.8) implies that

□

6.   Numerical scheme and results

In this section, we present a normal approximation for the stochastic $\mathscr{S}\mathscr{I}\mathscr{V}\mathscr{R}$ model (2.2). We consider the parameter values $\mathscr{S}_{(0)} = 50$, $\mathscr{I}_{(0)} = 10$, $\mathscr{V}_{(0)} = 20$, $\mathscr{R}_{(0)} = 30$, $\beta = 0.05$, $A = 0.0123$, $\rho = 0.24$, $\mu = 0.112$, $\alpha = 0.3$, $\delta_1 = 0.011$, and $\gamma = 0.312$ from the literature ^[13].

For the mathematical analysis, we first describe the corresponding deterministic numerical model (2.1). The discretization of this model leads to the following numerical scheme ^[23]:

Moving forward, we delve into the graphical depiction of the $\mathscr{SIVR}$ model dynamics captured by (2.2). The following figures portray how the system's behavior evolves under different stochastic perturbation strengths, offering insights into the model's robustness and potential real-world implications.

The generated Figure 1 illustrates the evolution of the susceptible $\mathscr{S}$, infected $\mathscr{I}$, vaccinated $\mathscr{V}$, and recovered $\mathscr{R}$ populations under different noise levels $\sigma_1$. The results demonstrate how stochastic perturbations influence disease dynamics, leading to either disease persistence or extinction, depending on the intensity of fluctuations.

For low noise levels $\sigma_2 = 0.0, 0.01$, the infected population $\mathscr{I}$ remains persistent over time, indicating that the disease does not die out but fluctuates around an endemic state. However, for higher noise levels $\sigma_2 = 0.02, 0.03$, fluctuations become more pronounced, and in several instances, the infected population decreases toward zero, suggesting that stochastic effects can drive the disease to extinction.

The graphical results further illustrate and confirm these findings. The susceptible population $\mathscr{S}$ increases slightly over time as vaccination and recovery reduce the number of infected individuals. The infected population $\mathscr{I}$ decreases more rapidly when stochastic effects are stronger, illustrating how noise can drive infection levels toward extinction. Meanwhile, the vaccinated $\mathscr{V}$ and recovered $\mathscr{R}$ populations fluctuate but remain relatively stable, highlighting the role of vaccination in mitigating infection spread. These results highlight the importance of incorporating stochasticity into epidemic models, as real-world uncertainties in transmission and intervention strategies can lead to substantial variations in disease outcomes.

The Figure 2 illustrates the trajectories of the system in the susceptible-infected ($\mathscr{S}$-$\mathscr{I}$) plane, demonstrating how the disease evolves over time under different noise levels $\sigma_2$. These trajectories provide valuable insights into the long-term behavior of the epidemic and the role of stochastic perturbations in shaping disease dynamics.

For low noise levels $\sigma_2 = 0.0, 0.01$, the infected population $\mathscr{I}$ exhibits cyclic behavior, meaning the disease persists within the system without fading out completely. The trajectory stabilizes over time, implying a long-term endemic state, where fluctuations do not lead to disease elimination.

In contrast, for higher noise levels $\sigma_2 = 0.02, 0.03$, the infected population $\mathscr{I}$ gradually declines towards zero, indicating disease extinction due to the increasing influence of stochastic effects. The randomness disrupts sustained transmission, causing occasional fluctuations that push infection levels below the threshold required for persistence, ultimately leading to eradication.

From a biological and public health perspective, these findings highlight the critical role of stochastic effects in disease dynamics. When stochasticity is weak, natural fluctuations are insufficient to eliminate the infection, meaning stronger intervention strategies such as vaccination campaigns and public health measures are required to prevent long-term persistence. However, when stochastic effects are strong, random fluctuations alone can drive the disease toward extinction, emphasizing the importance of considering variability in epidemic modeling. These results reinforce the need for adaptive disease control policies that account for uncertainties in disease transmission and intervention efficacy.

After perturbation, the path of $\mathscr{S}_{(t)}$ in model (2.2) demonstrates noticeable changes as $\sigma_1$ and $\sigma_2$ vary from 0.00 to 0.03. Figure 3 illustrates these changes, showing the trajectories of $\mathscr{S}(t)$ for different values of $\sigma_1$ and $\sigma_2$. At $\sigma_1 = \sigma_2 = 0.00$, the path remains close to the unperturbed trajectory, indicating minimal stochastic influence. As $\sigma_1$ and $\sigma_2$ increase, the trajectory exhibits greater fluctuation, reflecting the increased impact of stochasticity on the system. The perturbed paths diverge further from the unperturbed trajectory, suggesting that higher levels of noise lead to greater variability in the susceptible population. This observation underscores the importance of accounting for stochastic effects when modeling infectious disease dynamics.

The evolution of $\mathscr{I}_{(t)}$ in model (2.2) exhibits distinct patterns as the levels of stochasticity, characterized by $\sigma_1$ and $\sigma_2$, increase. When both $\sigma_1$ and $\sigma_2$ are 0.00, $\mathscr{I}_{(t)}$ follows a smooth trajectory, mirroring the deterministic dynamics of the system. However, as stochastic perturbations are introduced with $\sigma_1 = \sigma_2 = 0.01, 0.02, 0.03$, the behavior of $\mathscr{I}_{(t)}$ becomes more erratic, displaying larger fluctuations around the deterministic solution. These fluctuations indicate the impact of randomness on the infected population's dynamics, emphasizing the need to consider stochastic factors in modeling infectious diseases for more accurate predictions. This trend is visually depicted in Figure 4.

The behavior of the vaccinated population, $\mathscr{V}_{(t)}$, in model (2.2) under varying levels of stochasticity, represented by $\sigma_1$ and $\sigma_2$, reveals notable trends. When $\sigma_1$ and $\sigma_2$ are both 0.00, $\mathscr{V}_{(t)}$ follows a deterministic trajectory, unaffected by stochastic fluctuations. However, as stochastic perturbations are introduced with $\sigma_1 = \sigma_2 = 0.01, 0.02, 0.03$, the behavior of $\mathscr{V}_{(t)}$ becomes more erratic, showing increased variability around the deterministic solution. This increased variability reflects the impact of randomness on the vaccinated population's dynamics, highlighting the importance of incorporating stochastic elements in modeling vaccination strategies for infectious diseases. This trend is visually illustrated in Figure 5.

The dynamics of the recovered population, $\mathscr{R}_{(t)}$, in model (2.2) are also influenced by stochastic perturbations, as indicated by varying levels of $\sigma_1$ and $\sigma_2$. Initially, when both $\sigma_1$ and $\sigma_2$ are 0.00, $\mathscr{R}_{(t)}$ follows a smooth, deterministic trajectory. However, as stochasticity is introduced with $\sigma_1 = \sigma_2 = 0.01, 0.02, 0.03$, the behavior of $\mathscr{R}_{(t)}$ becomes more erratic, exhibiting increased fluctuations around the deterministic solution. These fluctuations underscore the impact of randomness on the dynamics of the recovered population, suggesting that stochastic factors play a significant role in shaping the long-term behavior of the system. Figure 6 visually depicts this trend.

7.   Conclusions

In conclusion, the stochastic $\mathscr{SIVR}$ model (2.2) provides valuable insights into the dynamics of infectious diseases, incorporating both deterministic and stochastic elements. The analysis demonstrates that stochastic perturbations, represented by $\sigma_1$ and $\sigma_2$, have a significant impact on the long-term trends of the susceptible, infected, vaccinated, and recovered populations. These perturbations introduce increased variability in the population dynamics, highlighting the importance of stochasticity in disease transmission and control. The findings underscore the need to consider stochastic factors in epidemiological models to enhance our understanding and management of infectious diseases. Future research could focus on refining the model's parameters and validating its predictions using empirical data.

Author's contributions

Shah Hussain: Writing original draft, software, methodology, formal analysis; Naveed Iqbal: Writing review and editing, visualization, validation, supervision; Elissa Nadia Madi: Investigation, formal analysis, conceptualization; Thoraya N. Alharthi: Visualization, validation, supervision; Ilyas Khan: Methodology, analysis, writing review and editing, supervision. All the authors have equal contributions in this article.

Use of Generative-AI tools declaration

All authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

All authors declare no conflict of interest in this paper.

Acknowledgments

The authors extend their appreciation to the Deanship of Postgraduate Studies and Scientific Research at Majmaah University for funding this research work through the project number (R-2025-1535).

The authors are also thankful to the Deanship of Graduate Studies and Scientific Research at the University of Bisha for supporting this work through the Fast-Track Research Support Program.