Stochastic analysis of Mpox epidemiology with vaccination strategies and environmental persistence

Categories: Disease & Virus

December 23, 2025

Stochastic analysis involves modeling systems that incorporate inherent randomness, where outcomes are influenced by probabilities and uncertainty. This approach is particularly useful for studying phenomena such as disease spread, where variability plays a key role. On the other hand, deterministic analysis focuses on systems with predictable behavior, where future states are determined by initial conditions and governing equations, without any randomness involved. This method provides precise predictions in scenarios where the system’s behavior can be accurately modeled without considering uncertainty. Both approaches offer valuable insights depending on the nature of the system being studied.

Incorporating stochasticity is especially important for diseases like Mpox, where outbreaks may be sporadic, case numbers are often low, and transmission routes can vary. The random nature of contact patterns, environmental exposures, and behavioral factors introduces uncertainty that can dramatically influence outbreak dynamics. By modeling these random effects, the stochastic framework provides a more nuanced understanding of disease spread, helping to identify the conditions under which the disease might persist, fade out, or re-emerge unexpectedly. This added layer of analysis enhances the predictive power and reliability of the model in public health planning and risk assessment.

Analysis of continuous-time Markov chain (CTMC)

Formulation of CTMC model

Mpox transmission dynamics are inherently stochastic, influenced by factors like human contact patterns, viral shedding, immune responses, and control measures such as vaccination. The Continuous Time Markov Chain (CTMC) model captures this variability by simulating state transitions over continuous time. It accounts for individual level heterogeneity and time dependent events such as disease progression, recovery etc. The CTMC model’s Markov property assumes future states depend solely on the current state, providing a realistic representation of Mpox transmission dynamics. In this section we extend the deterministic ODE system for Mpox model (3) into a CTMC (continuous time Markov chain) model. For convenience same notation for state variables and model parameters are used⁵². Let, $\tilde{{\mathcal {L}}}(t)=(S_h(t), V(t), I_h(t), R_h(t), S_r(t),I_r(t), C_s(t)), {\rm {the\,\, random\,\, variables,}}$ $S_h, V, I_h, R_h, S_r,I_r, C_s\in {0, 1, 2,…, N},$ and $t\in [0,\infty ).$ Assume that there is a sufficiently small time interval $\Delta t > 0$ such that at most one event can occur within this interval. Let, $\Delta \tilde{{\mathcal {L}}}(t)=\left( \Delta S_h(t), \Delta V(t), \Delta I_h(t), \Delta R_h(t), \Delta S_r(t), \Delta I_r(t), \Delta C_s(t)\right)$, where, $\Delta S_h=S_h(t+\Delta t)-S_h(t)$ etc. Various events along with their corresponding transition probabilities are given in Table 5.

Table 5 Possible changes of the state variables during the time interval $\Delta t$ with corresponding transition probabilities for the CTMC model.

Analysis of stochastic solutions

The CTMC (Continuous-Time Markov Chain) model is introduced to capture the inherent stochasticity and random variations in the disease transmission and recovery processes, which deterministic models cannot account for. This approach provides a more realistic representation of the dynamics.

In the beginning, starting with around 695 individuals, started decreasing rapidly as individuals move to the V class, and $I_h$ class. The rate of decline slows down within 12-50 days and stabilizes at a low level, with minor fluctuations after 125 days (Fig. 22a). Initially the vaccinated population rises quickly, reaching a peak of about 520-580 individuals. Over time (after 50 days), the population begins to decline across all paths and stabilizes at a low level, with minor fluctuations after 200 days (Fig. 22b). Figure 22c illustrates that starting at a very low level, the infected population increases rapidly, peaking around 200-270 individuals between 70-100 days across all paths and after the peak, the population steadily declines, with all paths. As time progresses, the recovery count increases steadily due to the recovery of infected individuals after day 50. The rate of recovery continues to rise near day 180 as the population moves to the recovered state (Fig. 22d). The deterministic model (black dashed curve) closely follows the average trend of the stochastic paths, and the stochastic paths displaying variations due to randomness in states.

Figure 23a depicts that the susceptible rodent population decreases rapidly at first and stabilizes after day 60 at low levels due to transitions to infected class $I_r.$ The number of infected rodents increases initially as susceptible rodents become infected, leading to a rapid increase in cases. After reaching a peak of about 80-120 individuals between day 25–50, the infected rodents gradually declines, leading to a steady downward trend in infection levels over time (Fig. 23b). In Fig. 23c, contamination steadily increases as the infection spreads. This rise in contamination is primarily driven by the virus shedding rate from infected humans and rodents, which continuously release the virus into the environment. As the number of infected humans and rodents increases, the environmental viral load accumulates, leading to higher contamination levels. However, from day 150 onward, contamination begins to decline. As the number of infected humans and rodents declines over time, the overall virus shedding into the environment decreases, leading to a reduction in environmental contamination. This decrease is further driven by an increased virus decay rate, influenced by environmental factors such as temperature, humidity, and UV exposure, as well as intervention measures like disinfection and sanitation efforts. These factors accelerate viral degradation and reducing environmental contamination and overall transmission risk. Overall, the deterministic model provides an average trend, while the CTMC paths reveal the natural randomness in rodent infection and environmental contamination dynamics.

Stochastic differential equation (SDE) and disease dynamics

Formulation of SDE model

To incorporate the uncertainty and randomness in disease dynamics, we construct a stochastic differential equation (SDE) model derived from the CTMC framework. Unlike traditional deterministic models that describe disease spread through fixed patterns, SDE models capture real-world variability by introducing random fluctuations within the system^53,54. In general, the SDE system is formulated as:

$$\begin{aligned} d\tilde{{\mathcal {L}}}={\mathcal {U}}(t)dt+{\mathcal {G}}(t)dW(t), \end{aligned}$$

(22)

where ${\mathcal {U}}(t)$ refers to the deterministic drift vector, and ${\mathcal {G}}(t)$ is the diffusion matrix satisfies ${\mathcal {G}}(t){\mathcal {G}}^T(t)={\mathcal {C}}(t),$ with ${\mathcal {C}}(t)$ being the covariance matrix reflecting the magnitude of uncertainty. Specifically, ${\mathcal {G}}$ is structured as ${\mathcal {G}}= (\mathcal {\breve{A}}|\mathcal {\breve{B}}),$ where $\mathcal {\breve{A}}$ and $\mathcal {\breve{B}}$ are both $7\times 8$ matrices, respectively as follows:

$$\begin{aligned} & \mathcal {\breve{A}}= \begin{bmatrix} \sqrt{\mu _h N_h} & -\sqrt{\xi _1 S_h}& \sqrt{\xi _2 V}& -\sqrt{\mu _h S_h} & 0& -\sqrt{\lambda _h S_h}& 0& 0\\ 0& \sqrt{\xi _1 S_h}& -\sqrt{\xi _2 V}& 0& -\sqrt{\mu _h V} & 0& 0& 0 \\ 0& 0& 0& 0& 0& \sqrt{\lambda _h S_h}& -\sqrt{\Gamma I_h}& -\sqrt{(\mu _h+\mu _{hd})I_h}\\ 0& 0& 0& 0& 0& 0& \sqrt{\Gamma I_h}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ \end{bmatrix}, \end{aligned}$$

(23)

$$\begin{aligned} & \mathcal {\breve{B}}= \begin{bmatrix} 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ -\sqrt{\mu _hR_h}& 0& 0& 0& 0& 0& 0& 0\\ 0& \sqrt{\mu _r N_r}& -\sqrt{\lambda _r S_r}& -\sqrt{\mu _r S_r}& 0& 0& 0& 0\\ 0& 0& \sqrt{\lambda _r S_r}& 0& -\sqrt{\mu _r I_r}& 0& 0& 0\\ 0& 0& 0& 0& 0& \sqrt{\chi _hI_h}& \sqrt{\chi _rI_r}& -\sqrt{\mu _sC_s}\\ \end{bmatrix}. \end{aligned}$$

(24)

The system of stochastic differential equations for model (3) is,

$$\begin{aligned} {\left\{ \begin{array}{ll} dS_h(dt) ={\mathcal {U}}_1dt+\sqrt{\mu _h N_h}W_1 +\sqrt{\xi _2 V}W_3-\sqrt{\lambda _hS_h}W_6-\sqrt{\xi _1 S_h}W_2-\sqrt{\mu _h S_h}W_4,\\ dV(dt) = {\mathcal {U}}_2dt+\sqrt{\xi _1 S_h}W_2- \sqrt{\xi _2 V}W_3- \sqrt{\mu _h V}W_5 ,\\ dI_h(dt) = {\mathcal {U}}_3dt+\sqrt{\lambda _hS_h}W_6-\sqrt{\Gamma I_h}W_7-\sqrt{(\mu _h+\mu _{hd}) I_h}W_8,\\ dR_h(dt) ={\mathcal {U}}_4dt+\sqrt{\Gamma I_h}W_7-\sqrt{\mu _h R_h}W_9, \\ dS_r(dt) = {\mathcal {U}}_5dt+\sqrt{\mu _r N_r}W_{10}- \sqrt{\lambda _rS_r}W_{11}-\sqrt{\mu _r S_r}W_{12},\\ dI_r(dt) = {\mathcal {U}}_6dt+\sqrt{\lambda _rS_r}W_{11}-\sqrt{\mu _r I_r}W_{13},\\ dC_s(dt) = {\mathcal {U}}_7dt+\sqrt{\chi _hI_h}W_{14}+\sqrt{\chi _rI_r}W_{15}- \sqrt{\mu _s C_s}W_{16}. \end{array}\right. } \end{aligned}$$

(25)

where the column vector ${\mathcal {U}}=({\mathcal {U}}_1, {\mathcal {U}}_2, {\mathcal {U}}_3, {\mathcal {U}}_4, {\mathcal {U}}_5, {\mathcal {U}}_6, {\mathcal {U}}_7)^T$ is

$$\begin{aligned} {\mathcal {U}}= \begin{bmatrix} \Lambda _h+\xi _2 V-\lambda _hS_h-(\xi _1+\mu _h) S_h\\ \xi _1 S_h- (\xi _2+\mu _h) V\\ \lambda _hS_h-\Gamma I_h-(\mu _h+\mu _{hd}) I_h\\ \Gamma I_h-\mu _h R_h \\ \Lambda _r- \lambda _rS_r-\mu _r S_r\\ \lambda _rS_r-\mu _r I_r\\ \chi _hI_h+\chi _rI_r- \mu _s C_s \end{bmatrix}, \end{aligned}$$

(26)

and $W(t)=(W_1(t),…,W_{16}(t))^T$ is a vector of sixteen independent Wiener processes, corresponding to the sixteen events represented in Table 5. Here $W_1$, $W_2$,…,$W_8$ corresponds to each of the column vector of $\mathcal {\breve{A}}$ respectively and $W_9$,…,$W_{16}$ corresponds to each of the column vector of $\mathcal {\breve{B}}$ respectively. The SDE system (22) is numerically solved using the Euler-Maruyama method²⁴, a standard and effective method for integrating stochastic differential equations. A constant time step size of 0.1 is used to ensure a balance between numerical accuracy and computational efficiency. To verify convergence and robustness of the simulation results, multiple runs are performed under identical parameter conditions. This approach enabled the model to reflect realistic variability in Mpox transmission and control dynamics, capturing fluctuations in peak timing, magnitude, and disease fade-out scenarios.

Analysis of stochastic solutions

The susceptible human class $(S_h)$ represents individuals at risk of contracting Mpox. Figure 24a depicts that, initially, the number of susceptible humans starts around 635 individuals, then drops sharply in the first 10 days as the disease begins to spread. From day 10 to 100, the decline continues more gradually, showing noticeable fluctuations due to stochastic effects. By day 100, the susceptible population decreases below 100 in all five paths. After day 150, the numbers stabilize at a low level, fluctuating between 40 and 70. The variability in the paths illustrates the impact of randomness on the timing and rate of infection, even though all trajectories follow a similar trend.

Vaccinated individuals are less likely to contract the infection, though the level of protection depends on vaccine efficacy; imperfect vaccination may still contribute to disease transmission and impact overall outbreak dynamics. Figure 24b depicts that, initially, the number of vaccinated humans rapidly increases and peaks between day 20 and day 50, reaching about 480 to 540 individuals. This early rise reflects intense vaccination efforts or immune response. However, after peaking, the population steadily declines as vaccinated individuals either move to other compartments or due to waning immunity. By day 150, the vaccinated population falls between approximately 100 to 150 individuals and then fluctuates near these values. Different stochastic paths show varying rates of decline, indicating diverse individual responses and intervention effectiveness.

Figure 25a presents the stochastic evolution of the infected human population. All five paths begin near zero and increase steadily, peaking between days 90 and 100. Peak infections range from 170 to 275, demonstrating variation among stochastic outcomes. After the peak, infection levels gradually decline. The deterministic curve lies near the center of the paths, indicating that while average dynamics are captured, SDE reveal possible deviations due to random fluctuations.

The recovered class comprises individuals who have successfully overcome the infection and gained immunity. Figure 25b shows that, initially, this group is empty, bur recovery accelerates after day 50 as more individuals transition out of the infected class. From day 50 to day 150, the recovered population increases significantly, reaching values between 110 and 130 by day 170. The upward trend eventually levels off as the infection subsides and fewer individuals remain to recover. All five stochastic paths reflect consistent long term behavior, with some showing faster recoveries than others. Figures 24 and 25 show that the deterministic curve remains close to the mean of the stochastic paths, validating the model’s reliability under uncertainty.

Figure 26a shows the dynamics of the susceptible rodent population under five SDE solution paths and a deterministic curve. All paths show a rapid decline from an initial value of around 95 to below 10 within the first 40 days, suggesting a swift disease spread. After this point, the population stabilizes at a low level, fluctuating slightly due to stochastic effects. The deterministic solution closely follows the average trend of the stochastic paths.

Figure 26b illustrates the progression of infected rodents over time using five SDE paths. Infections increase rapidly, peaking between days 35 and 60, with maxima between 70 and 109. After peaking, infections decline slowly and display more variability than in the human population. Some paths remain around 40 even after 200 days, indicating prolonged infection persistence in rodents. The deterministic solution again provides a mean trajectory, surrounded by the stochastic paths.

Figure 26c compares five sample paths of the SDE model with the deterministic solution for environmental contamination in Mpox transmission. All paths peak between days 110-160 with virus loads ranging from 480 to 610. The deterministic path smoothly peaks around day 120 at 510. Stochastic fluctuations reflect real-world variability, emphasizing environmental viral load’s critical role in sustaining transmission dynamics.

The stochastic differential equation (SDE) framework offers notable advantages over the continuous-time Markov chain (CTMC) model by incorporating continuous-time stochasticity to reflect natural variability in population dynamics. This results in smoother, biologically plausible transitions that more accurately represent real-world epidemic behavior, particularly in large populations. In contrast, the CTMC model’s discrete transitions may oversimplify complex processes. Thus, the SDE approach provides a more robust and realistic foundation for epidemiological modeling and intervention analysis.

Remark 4

The analysis of the existence and uniqueness of positive solutions, together with the stochastic dynamics of both the disease-free and endemic equilibria, is presented in the Supplementary Material.

Source link