A novel epidemiologically informed particle filter for assessing epidemic phenomena. Application to the monkeypox outbreak of 2022

The accurate prediction of the evolution of an epidemic event, plays an important role in reducing the transmissibility of infectious diseases. Contagious illnesses like Ebola, hepatitis, influenza, dengue etc. are constantly affecting more and more people every day, resulting in widespread health crises especially in developing nations. The 2022 outbreak of monkeypox (mpox) was first reported in the UK in May, and new cases were rapidly reported in other European countries [1]. Although Europe and the U.S. remain the epicenters of the outbreak, cases have since been detected in Asia, Oceania and elsewhere in America [2]. Structurally, this virus is related to the smallpox virus and infects the host in a similar way, while the symptoms of mpox are less severe [3].

The situation is evolving, and the world health organization claim that there will be more cases of mpox identified as surveillance expands in non-endemic countries. Current available evidence suggests that those who are most at risk are those who have had close physical contact with someone with mpox, while they are symptomatic [4]. Recent papers propose various values for the basic reproduction number ${R_0}$—the expected number of secondary cases produced by an already infected individual [5]—according to the examined region. Schneider & Eichner [6] stated that the ${R_0}$ in African villages is usually less than 1, while another estimation of a mpox's clade in the Congo Basin had shown a ${R_0}{\text{ }}$ between 1.46 and 2.67. Although, values greater than 3 should be deemed as highly pessimistic.

Kaler et al [7], mention that the reported ${R_0}$ value of mpox is between 1.10 and 2.40 in countries where exposure to Orthopoxvirus species (variola, mpox, cowpox, vaccinia) is negligible, while in Italy a ${R_0}$ of 2.43 was proposed during the early stages of the disease spread that attenuated during the following months. Du et al [8], suggest that the effective reproduction number ${R_e}$ in the United States, where we encounter the majority of reported cases, is around 1.95, while European countries displayed lower values. Riopelle et al [9] declares that data from both the United States and the Netherlands were examined, and the incubation period was estimated around 7.6 days, and around 8.5 days using only the cases in the Netherlands.

The investigation of the disease dynamics using mathematical modeling could accompany policies to suppress the diseases' outbreaks. The evolution of infectious illnesses is often described through deterministic epidemiological models, where the susceptible-infected-recovered (SIR) model represents the most simplistic approach [10, 11]. Previous studies investigated the epidemics via the SIR model [12], or its extensions, e.g. the susceptible-infective-recovered-susceptible [13], the susceptible-exposed-infective-recovered (SEIR) [14, 15] or the susceptible-exposed-infective-recovered-deceased (SEIRD) schemes [16].

Although the deterministic methods that accompany the previously mentioned epidemic models provide short-term accuracy and they do not adequately describe the fluctuating dynamics of an outbreak for long time periods. Furthermore, due to the misclassification rate of the diagnostic tests and the delays in records, it was realized that there exist uncertainties in the reported data [17, 18]. In addition, a deterministic model aiming to describe all the potential states and transitions of an epidemic will inevitably be characterized by a computationally expensive structure that could result into significant overfitting, without capturing the time varying dynamics of the pandemic. Thus, by utilizing a stochastic approach, the variability in states and parameters could also be described.

Sebbagh and Kechida [19], employed a SEIRD-extended Kalman filter (EKF), including the parameters of the model in the updating process of the EKF. Ndanguza et al [20] included in the EKF-SEIR model the estimation of the epidemic parameters, while Costa et al [21] combined also an EKF with a SEIR scheme to simulate an epidemic outbreak. Calvetti et al [22], mixed the susceptible-exposed-asymptomatic-infectious-recovered structure with the methodology of the particle filter for the estimation of the early-stages of the spread of COVID-19 in Ohio and Michigan. Storvik et al [23], proposed a sequential Monte-Carlo methodology for the estimation of the effective reproduction number for contagious diseases, applied on daily Norwegian COVID-19 data. Johndrow et al employed Poisson likelihood functions for the modeling of the cumulative number of deaths from SARS-Cov-2 [24]. Papageorgiou and Tsaklidis [25] propose an extended epidemic model in parallel with an unscented Kalman Filter for the estimation of the prevalence of SARS-Cov-2 in France. Endo et al [26], presented a tutorial of particle Monte-Carlo Markov Chain methodologies for the modeling of epidemic outbreaks.

In this article, an epidemiologically informed particle filter (EI-PF) for the estimation of mpox evolution is proposed. Our model combines the structure of the deterministic SEIRD model, with the methodology of sequential importance resampling particle filter aiming to exploit the information provided by the daily records. This implementation increases significantly our model's ability to encapsulate the stochasticity that accompanies the population dynamics of infectious illnesses, compared to the widely used deterministic schemes. The addition of a time-varying transmissibility $\left( {{\beta _t}} \right)$ and mortality $\left( {{\mu _t}} \right)$ rate, presents enhanced adaptability to the observations while allowing for a fluctuating reproductive number that is characteristic of mpox, as in the abovementioned case of Italy.

According to the SEIRD structure, we employ discrete probability distributions for the transitions of particles during the importance sampling step, while log-normal distributions of zero mean define the steps of ${\beta _t}$ and ${\mu _t}$, for $t = 1, \ldots ,{\text{ }}T$. Also, a bivariate Poisson distribution is employed to connect the observations and the hidden states; a selection that is new to the literature. During the filtering procedure and the iterative determination of particle weights, we underline the addition of extra parameters—penalty factors—to the standard likelihood function that penalize irregular behaviors, namely trajectories that are not compatible with the specificities of epidemiological cases. The inclusion of penalty factors stands out as a key novelty in this analysis. This addition not only produces estimates aligned with the specifics of epidemic spread but also improves the fitting and predictive capabilities of the methodology introduced. This phenomenon may arise due to the resampling step of the particle filter algorithm, leading to significant irregularities in the final estimations.

Although, we cannot avoid the implementation of the resampling process -due to its essential contribution to the confrontation of the degeneracy problem associated with particle filters- leading to the employment of penalty factors. The addition of the abovementioned parameters, improves the fitting/prediction of new cases and deaths, as shown in the results section. More specifically, the proposed approach takes advantage of the derivative's sign of deceased and recovered cases, since decreasing tendencies in the estimations of the corresponding cumulative numbers of these quantities should be eliminated.

We conducted tests on the mpox data, specifically on the daily infections and cumulative deaths. Our study demonstrates the appropriateness of the proposed stochastic algorithm over deterministic models, highlighting the importance of additional penalty components in achieving accurate fitting. Given that the available datasets for this epidemic only provide information on the number of infected and deceased individuals -which is a common limitation in various cases of epidemics- we believe that the SEIRD system is appropriate for studying the current outbreak. It is worth noting that using a more complex epidemiological model could result in overfitting, as it would introduce extra states that are not supported by the available data [25]. Also, the nature of the recent outbreak does not provide enough evidence for the inclusion of additional compartments. Our methodology has the potential to describe other contagious diseases like influenza, hepatitis, diphtheria, and more. This is due to its simple, yet adaptive structure.

The rest of the paper is organized as follows: in paragraph 2, we describe the mathematical scheme of the proposed particle SEIRD model reinforced with information derived from epidemic characteristics. In paragraph 3, we validate the fitting-predictive efficiency of the proposed model based on the daily cases and total deaths caused by mpox, starting from the 25th of June 2022. In paragraph 4 we discuss and in paragraph 5 we conclude the major findings of our study giving special emphasis on the advantages of the proposed stochastic epidemiological methodology.

2.1. The deterministic SEIRD compartmental model

In this section, we briefly present the structure of the deterministic SEIRD model. This step is necessary to understand the fundamental population dynamics that we later employ for the construction of the stochastic extension based on particle filtering. The SEIRD model contains five states, as it is mentioned in the introduction, and four transition rates that connect them. It consists of a system of five ordinary differential equation (1) and represents the possible transitions of the individuals inside the model's compartments. To our knowledge, there are limited case studies on the transmission of the mpox virus to newborns, and these instances have been linked to a significant risk of miscarriage and fetal loss [27]. Additionally, there is not sufficient evidence regarding the transmission of mpox to infants during the current outbreak in the United States, making it challenging to confidently determine an appropriate influx rate. Moreover, we have assumed that, for the population under study, transportation's expected inputs and outputs are roughly balanced, and thus considering a closed population. Lastly, minor impact transition rates contribute to the system's noise/uncertainty, which the particle filtering methodology effectively addresses [22]. This represents a significant advantage of stochastic approaches, allowing the use of simpler and more computationally robust compartmental schemes [25].

Figure 1 presents the transitions between the states of the epidemiological model, while table 1 provides the definitions for the notation of the deterministic SEIRD (2),

$$\begin{align} \frac{{{\text{d}}S\left( t \right)}}{{{\text{d}}t}} & = - \frac{{\beta S\left( t \right)I\left( t \right)}}{N}\nonumber\\ \frac{{{\text{d}}E\left( t \right)}}{{{\text{d}}t}} &= \frac{{\beta S\left( t \right)I\left( t \right)}}{N} - \alpha E\left( t \right)\nonumber\\ \frac{{{\text{d}}I\left( t \right)}}{{{\text{d}}t}} &= \alpha E\left( t \right) - \kappa I\left( t \right) - \mu I\left( t \right),\nonumber\\ \frac{{{\text{d}}R\left( t \right)}}{{{\text{d}}t}} &= \kappa I\left( t \right)\nonumber\\ \frac{{{\text{d}}D\left( t \right)}}{{{\text{d}}t}} &= \mu I\left( t \right)\nonumber\\ \forall t \unicode{x2A7E} 0. \end{align} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

Table 1. Parameters and state definitions of the SEIRD.

Symbol	Parameters/States
S	Susceptible
E	Exposed
I	Infective
R	Recovered
D	Deceased
$\alpha$	Incubation rate
${\text{ }}\beta$	Infection or transmission rate
${\text{ }}\kappa$	Recovery rate of infective cases
${\text{ }}\mu$	Mortality rate of infective cases

The selection of the SEIRD model for our analysis is not random, as it was based on the present spread of mpox in the population and the available data. For instance, we avoid the inclusion of extra states that are not supported by daily records, like quarantine, hospital, or intensive care unit admissions, as it is not feasible to test the effectiveness of our modeling procedure on them. In addition, there is no strong evidence in the literature, during the present epidemic outbreak, which supports the employment of extra compartments. Usually system (1) is solved with the aid of iterative algorithms like the 4th-order Runge–Kutta method that provides estimations with satisfactory efficiency [28].

Moreover, in the occasion of COVID-19, we know that the immunity caused by infection in recovered individuals attenuates after some months. This leads to the necessity of another transition rate that connects the compartments of susceptible and recovered cases. However, the recent reports about mpox do not support the possibility of re-infection, especially during the early stages of the outbreak. This fact leads to a more robust structure that is characterized by a lower complexity.

On the other hand, we believe that the employment of simplifications of the SEIRD scheme, such as other variations of the SIR model, seem inappropriate in our case. It is evident that the inclusion of the compartment of deceased cases is necessary as the daily observations are characterized by data concerning deaths from mpox. In addition, new deaths have a significant role on the selection of the likelihood function, as it is underlined in section 2.2. Finally, the incorporation of the exposed state aligns with the properties of the present epidemic outbreak due to the existence of incubation periods [9], as we have mentioned in the introduction section.

However, as we also discuss in the results section, the deterministic approach is not able to efficiently describe the stochasticity that accompanies an epidemiological phenomenon. Therefore, the transition to a stochastic point of view appears to be necessary. The employed hybrid, stochastic SEIRD particle filter aims to overcome the drawbacks of conventional approaches.

2.2. An epidemiologically informed sequential importance resampling particle filter

As we previously stated, the deterministic approaches do not efficiently handle the stochasticity that characterizes epidemics, as their transmissibility varies significantly especially during the early stages of the disease. The compartmental models like the extensions of the conventional SIR, are not perfect as they do not consider attributes such as the location where the outbreak takes place, the age of the population or additional specificities related to the immune system of infected individuals in general. An attempt to encapsule all this information to a system of differential equations would lead to unprofitable models due to high complexity. On the other hand, the daily observations are characterized by noise that is compatible with the nature of epidemics. This phenomenon is influenced by errors, delays in records of cases or deaths or even asymptomatic infections.

According to the abovementioned remarks, we aim to adapt the particle filtering methodology to the specificities of an epidemic's spread. The construction of a hybrid particle epidemiological model should effectively combine the properties of the proposed compartmental model with the available observations, while managing the drawbacks of these two sources. As a result, we aim to encapsulate to the standard SEIRD, any available information regarding the evolution of the epidemic such as the daily numbers of new cases and the total of deceased patients.

Our purpose is to estimate the posterior density of the trajectory of hidden states employing the information provided by the timeseries of measurements up to time $t$, $p({{\boldsymbol{x}}_{0:t}}|{{\boldsymbol{y}}_{0:t}})$, leading to an inverse problem. This solution is approximated by a sum of weighted samples/particles, where

$$\begin{align}p\left( {{{\boldsymbol{x}}_{0:t}}|{{\boldsymbol{y}}_{0:t}}} \right) = {\text{ }}\mathop \sum \limits_{i = 1}^M w_t^{\left( i \right)}\delta \left( {{{\boldsymbol{x}}_t} - {\boldsymbol{x}}_t^{\left( i \right)}} \right),{\text{ }}\end{align} \tag{ 2 }$$

$\delta \left( . \right)$ denotes the Dirac delta function and $M$ the number of particles. Since we cannot sample from the posterior itself, we sample from the so-called importance density function, ${f_{\boldsymbol{\theta}} }\left( {{{\boldsymbol{x}}_{0:t}}|{{\boldsymbol{y}}_{0:t}}} \right)$. The importance weights can be updated recursively as

$$\begin{align}w_t^{\left( i \right)} & = \frac{{p\left({\boldsymbol{x}}_{0:t}^{\left( i \right)}|{{\boldsymbol{y}}_{0:t}}\right)}}{{{f_{\boldsymbol{\theta }}}\left( {{\boldsymbol{x}}_{0:t}^{\left( i \right)}|{{\boldsymbol{y}}_{0:t}}} \right)}} \propto \frac{{p\left( {{{\boldsymbol{y}}_t}|{\boldsymbol{x}}_t^{\left( {\boldsymbol{i}} \right)}} \right)p\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)}} \right)p\left( {{\boldsymbol{x}}_{0:t - 1}^{\left( i \right)}|{{\boldsymbol{y}}_{0:t - 1}}} \right)}}{{{f_{\boldsymbol{\theta }}}\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{0:t - 1}^{\left( i \right)},{{\boldsymbol{y}}_{0:t}}} \right){f_{\boldsymbol{\theta }}}\left( {{\boldsymbol{x}}_{0:t - 1}^{\left( i \right)}|{{\boldsymbol{y}}_{0:t - 1}}} \right)}} \nonumber\\ & = w_{t - 1}^{\left( i \right)}\frac{{p\left( {{{\boldsymbol{y}}_t}|{\boldsymbol{x}}_t^{\left( {\boldsymbol{i}} \right)}} \right)p\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)}} \right)}}{{{f_{\boldsymbol{\theta }}}\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{0:t - 1}^{\left( i \right)},{{\boldsymbol{y}}_{0:t}}} \right)}}.\end{align} \tag{ 3 }$$

In practice, we are interested in the current filtered estimate $p({{\boldsymbol{x}}_t}|{{\boldsymbol{y}}_{0:t}})$. Hence ${f_{\boldsymbol{\theta}} }\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{0:t - 1}^{\left( i \right)},{{\boldsymbol{y}}_{0:t}}} \right),$ is assumed to be equivalent to ${f_{\boldsymbol{\theta}} }\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)},{{\boldsymbol{y}}_t}} \right).$ Then, the weights can be recursively calculated based on

$$\begin{align}w_t^{\left( i \right)} \propto w_{t - 1}^{\left( i \right)}\frac{{p\left( {{{\boldsymbol{y}}_t}|{\boldsymbol{x}}_t^{\left( {\boldsymbol{i}} \right)}} \right)p\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)}} \right)}}{{{f_{\boldsymbol{\theta }}}\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)},{{\boldsymbol{y}}_t}} \right)}}.\end{align} \tag{ 4 }$$

Considering that the states follow a first-order Markov process. The weights are then normalized while the symbol $\propto$, signifies the proportionality of the above quantities. The weights compensate for the fact that samples are drawn from the importance density ${f_{\boldsymbol{\theta}}}$ rather than the posterior pdf. The symbol $\propto$ denotes the proportionality, while index $\left( i \right)$ corresponds to the $i$th particle. The quantities $w_{t - 1}^{\left( i \right)},$ are the respective weights of the previous time step. We understand that the weights are proportional to the likelihood of observations based on each particle's values. Consequently, the procedure retains the particles that provide the highest likelihood.

Based on [29], we should specify the distribution that describes the law of transitions between the hidden states, ${f_{\boldsymbol{\theta}} }( {{\boldsymbol{x}}_t^{( i )}|{\boldsymbol{x}}_{t - 1}^{( i )},{y_t}} )$, where ${\boldsymbol{\theta}}$ denotes the parametric set of the system. Since we work with subpopulations, it seems necessary to use discrete distributions for the transitions between compartments, such as binomial and multinomial ones, depending on the number of outgoing transitions. Hence, for the compartments of susceptible and exposed individuals, we employ the following transition formulas

$$\begin{align}\left( {{d_{SE}}} \right)_t^{\left( i \right)}\sim B\left( {S_{t - 1}^{\left( i \right)},{ }1 - {e^{ - \frac{{\beta I_{t - 1}^{\left( i \right)}}}{N}{ }\Delta t}}{ }} \right),{ }\end{align} \tag{ 5 }$$

and

$$\begin{align}\left( {{d_{EI}}} \right)_t^{\left( i \right)}\sim B\left( {E_{t - 1}^{\left( i \right)},{\text{ }}1 - {e^{ - \alpha {\text{ }}\Delta t}}{\text{ }}} \right),{\text{ }}\end{align} \tag{ 6 }$$

where $1 - {e^{ - \frac{{\beta I_{t - 1}^{\left( i \right)}}}{N}{ }\Delta t}}$ and $1 - {e^{ - \alpha {\text{ }}\Delta t}}$ represent the probability of an individual transitioning from susceptible to exposed and from the exposed to infective state, respectively, during the interval $\left( {t,{\text{ }}t + \Delta t} \right)$. For the compartment of infective cases, that is characterized from more than one outgoing transitions, we utilize a multinomial distribution, given by

$$\begin{align}\left( {\begin{array}{*{20}{c}} {{d_{IR}}} \\ {{d_{ID}}} \end{array}} \right)_t^{\left( i \right)}\sim M\left( {I_{t - 1}^{\left( i \right)},\left( {\begin{array}{*{20}{c}} {\frac{\kappa }{{\kappa + \mu }}\left( {1 - {e^{ - \left( {\kappa + \mu } \right){\text{ }}\Delta t}}} \right)} \\ {\frac{\mu }{{\kappa + \mu }}\left( {1 - {e^{ - \left( {\kappa + \mu } \right){\text{ }}\Delta t}}} \right)} \end{array}{\text{ }}} \right)} \right).\end{align} \tag{ 7 }$$

In equation (5), the factors $\frac{\kappa }{{\kappa + \mu }}$, $\frac{\mu }{{\kappa + \mu }} \in \left( {0,1} \right)$ for $\kappa ,{\text{ }}\mu > 0,{\text{ }}$ can be considered as weights that ensure that the probabilities of the multinomial sum to unity. Thus, we note that

$$\begin{align}\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)}} \right)\sim {\boldsymbol{x}}_{t - 1}^{\left( i \right)} + {\text{ }}\left( {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} { - {\text{ }}\left( {{d_{SE}}} \right)_t^{\left( i \right)}} \\ {\left( {{d_{SE}}} \right)_t^{\left( i \right)} - \left( {{d_{EI}}} \right)_t^{\left( i \right)}} \\ {\left( {{d_{EI}}} \right)_t^{\left( i \right)} - \left( {{d_{IR}}} \right)_t^{\left( i \right)} - \left( {{d_{ID}}} \right)_t^{\left( i \right)}} \end{array}} \\ {\begin{array}{*{20}{c}} {\left( {{d_{IR}}} \right)_t^{\left( i \right)}} \\ {\left( {{d_{ID}}} \right)_t^{\left( i \right)}} \end{array}} \end{array}} \right),{\text{ }}\end{align} \tag{ 8 }$$

where ${\boldsymbol{x}}_t^{( i )} = {( {S_t^{( i )},{\text{ }}E_t^{( i )},{\text{ }}I_t^{( i )},{\text{ }}R_t^{( i )},{\text{ }}D_t^{( i )}} )^T},{\text{ }}$for $t \in \mathbb{N}$. The new observations, e.g. the daily new infections, can be considered as arrivals of a Poisson process in the (sub)intervals $\left( {t,{\text{ }}t + \Delta t} \right)$ [24, 30]. Thus, the respective times between new infections follow exponential distributions. For instance, if $\tilde T$ is the interarrival time between consecutive new infections and considering that the last infection took place at the time step $t$, we have that

$$\begin{align} & {\text{prob}}\{ {\text{the next infection occurs in the interval}}\left( {t,{\text{ }}t + \Delta t} \right)\}\nonumber\\ & \quad = {\text{prob}}\left\{ {{\text{the next transition to the infective state occurs in}}\left( {t,\,t + \Delta t} \right)} \right\} \nonumber\\ & \quad = {\text{prob}}\left\{ \tilde T \unicode{x2A7D} t + \Delta t|\tilde T > t\right\} = \frac{{{\text{prob}}\left\{ {t < \tilde T \unicode{x2A7D} t + \Delta t} \right\}}}{{{\text{prob}}\left\{ {\tilde T > t} \right\}}}\nonumber\\ & \quad = \,\frac{{1 - {e^{ - \alpha \left( {t + \Delta t} \right)}} - \left( {1 - {e^{ - \alpha t}}} \right)}}{{{e^{ - \alpha t}}}} = 1 - {e^{ - \alpha \,\Delta t}} .\end{align} \tag{ 9 }$$

The process of obtaining the remaining probability formulas follows a similar approach.

Regarding the selection of the likelihood function given the hidden states, denoted by $p( {{{\boldsymbol{y}}_t}|{\boldsymbol{x}}_t^{( i )}} )$, other analyses have employed univariate Poisson distributions whose rate have been based on the particles' estimated number of new cases $\alpha E_t^{\left( i \right)},$ at each time step $t \in \mathbb{N}$ [22], namely

$$\begin{align}{\text{ }}p\left( {{y_t}|{\boldsymbol{x}}_t^i} \right) = {\text{ }}\frac{{{{\left( {\alpha E_t^{\left( i \right)}} \right)}^{{y_t}}}}}{{{y_t}!}}{e^{ - \alpha E_t^{\left( i \right)}}},{\text{ }}\end{align} \tag{ 10 }$$

or the estimated daily diseased cases [31]. In (11), ${y_t}$ represents the reported count of new cases, that is later referred to as $\left( {NC_t^{{\text{obs}}}} \right)$.

Expanding this approach, in the current study, the conditional distribution connecting the time series of daily records with the estimated hidden states at time $t$, is a Bivariate Poisson distribution, due to the fact that new cases and deaths are usually correlated. Now, ${y_t} = {\left( {NC_t^{{\text{obs}}},{ }ND_t^{{\text{obs}}}} \right)^T},\forall t \in \mathbb{N}.$ Based on Karlis and Ntzoufras [32], a Bivariate Poisson distribution with time-varying parameters, ${\text{BivPois}}\left( {{\lambda _{1t}},{\lambda _{2t}},{\lambda _{3t}}} \right),$ has a joint probability function

$$\begin{align} \left( {{{\boldsymbol{y}}_t}|{\boldsymbol{x}}_t^{\left( i \right)}} \right) & \sim {\text{BivPois}}\left( {{\lambda _{1t}},{ }{\lambda _{2t}},{ }{\lambda _{3t}}} \right) \nonumber\\ & \quad = \exp \left\{ { - \left( {{\lambda _{1t}} + {\lambda _{2t}} + {\lambda _{3t}}} \right)} \right\}\frac{{\lambda _{1t}^{NC_t^{{\text{obs}}}}}}{{\left( {NC_t^{{\text{obs}}}} \right)!}}\frac{{\lambda _{2t}^{ND_t^{{\text{obs}}}}}}{{\left(ND_t^{{\text{obs}}}\right)!}}\nonumber\\ & \qquad \times \mathop {\mathop \sum \nolimits }\limits_{k = 0}^{\min \left( {NC_t^{{\text{obs}}},{ }ND_t^{{\text{obs}}}} \right)} \left( {\begin{array}{*{20}{c}} {NC_t^{{\text{obs}}}} \\ k \end{array}} \right)\left( {\begin{array}{*{20}{c}} {ND_t^{obs}} \\ k \end{array}} \right)k!{{\left( {\frac{{{\lambda _{3t}}}}{{{\lambda _{1t}}{ }{\lambda _{2t}}}}} \right)}^k}.\end{align} \tag{ 11 }$$

The parameters ${\lambda _{1t}} + {\lambda _{3t}}$ and ${\lambda _{2t}} + {\lambda _{3t}}$ describe the incidence rate of each Poisson random variable, namely ${\lambda _{1t}} + {\lambda _{3t}} = \left( {{d_{EI}}} \right)_t^{\left( i \right)}{\text{ }}$and ${\lambda _{2t}} + {\lambda _{3t}} = \left( {{d_{ID}}} \right)_t^{\left( i \right)},$ whereas the parameter of ${\lambda _{3t}}$ itself denotes the covariance between the two Poisson variables. For the estimation of the covariance—which is usually assumed to be constant ${\lambda _3}{\text{ }}$—researchers reside on recursive estimation techniques such as the Newton–Raphson algorithm or generalized least squares [33, 34]. Since the covariance of hidden states is unknown, we culminate in a conservative choice, where ${\lambda _3} = 0.05$, which provided the maximum log-likelihood estimate $\log \hat p\left( {{{\boldsymbol{y}}_{0:t}}} \right),{\text{ }}$ for the examined time interval of $t =$150 days.

As we aim to examine an epidemic outbreak, we could improve the particle selection process by providing additional information derived from the specificities of the physical phenomenon. Each particle represents a possible trajectory of the epidemic's spread in the population. As we investigate the number of deaths during the outbreak, it is important to eliminate particles that display fewer deceased or recovered cases as we move to future time steps. This phenomenon derives from the resampling step of the particle filter's algorithm.

During the importance sampling step, we sample from a density ${f_{\boldsymbol{\theta} }}({\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)})$, that is concentrated on a set where ${D_t} \unicode{x2A7E} {D_{t - 1}}$ and ${R_t} \unicode{x2A7E} {R_{t - 1}}$, due to equation (8). Although, the stochasticity of the resampling step may generate particles which represent less recovered and deceased cases on time step $t$ compared to the time step $t - 1$. This fact is evident in figure 4 of the results section, where we observe estimations of the standard particle filter that contradict with the specificities of an epidemic's evolution.

On the other hand, particles that display positive derivatives for the total number of deaths and recoveries with respect to time, should be rewarded. As a result, we add a penalty factor to (5), in order to satisfy the above properties leading to

$$\begin{align}w_t^{\left( i \right)} \propto { }w_{t - 1}^{\left( i \right)}\frac{{p\left( {{\boldsymbol{y}_t}|\boldsymbol{x}_t^{\left( i \right)}} \right)p\left( {\boldsymbol{x}_t^{\left( i \right)}|\boldsymbol{x}_{t - 1}^{\left( i \right)}} \right)}}{{{f_\theta }\left( {\boldsymbol{x}_t^{\left( i \right)}|\boldsymbol{x}_{t - 1}^{\left( i \right)},{\boldsymbol{y}_t}} \right)}} + {\omega _1}{\left( {1 - {\omega _2}} \right)^t}{ }{\frac{{{\text{d}}D}}{{{\text{d}}t}}^{\left( i \right)}} + {\omega _3}{\left( {1 - {\omega _4}} \right)^t}{ }{\frac{{{\text{d}}R}}{{{\text{d}}t}}^{\left( i \right)}}.\end{align} \tag{ 12 }$$

The ${\omega _j},j = 1, \ldots ,{\text{ }}4{\text{ }},$ are constant factors necessary to limit the influence of the derivatives of the $i$th particle ${\frac{{{\text{d}}D}}{{{\text{d}}t}}^{\left( i \right)}}{ }$ and ${\frac{{{\text{d}}R}}{{{\text{d}}t}}^{\left( i \right)}}$. The ${\omega _1},$ ${\omega _3} \in {\mathbb{R}^ + },$ and take small values aiming to balance the derivatives' influence with the expected likelihood levels. Quantities ${\omega _2}$ and ${\omega _4}$ has an attenuation role. The issue of negative derivatives is encountered more frequently during the early stages of the epidemic and smooths out as we progress to higher mortality levels. Therefore, we select $\left( {{\omega _2},{\text{ }}{\omega _4}} \right) \in {\left[ {0,1} \right]^2}$ according to the desired attenuation degree. In our case, ${\omega _1} = {\omega _3}$ and ${\omega _2} = {\omega _4}.$

Based on (8), the transition density of the hidden states for one time step is selected as the importance function ${f_{\boldsymbol{\theta}}}( {{\boldsymbol{x}}_t^{( i )}|{\boldsymbol{x}}_{t - 1}^{( i )},{{\boldsymbol{y}}_t}} ) = p({\boldsymbol{x}}_t^{( i )}|{\boldsymbol{x}}_{t - 1}^{( i )})$, which is also the most frequently utilized selection in the literature [22, 30–38]. It is important to mention that when the variance of the process model noise is very low, the samples will have minor variations, and over multiple iterations, the particles may converge to a single point in the state space. However, as demonstrated in the results section, this is not the case in our application, thanks to effective model initialization. Therefore, we can conclude that this selection of the importance function is appropriate when coupled with efficient calibration.

Considering the fluctuating dynamics of epidemic phenomena, we augment the state space by adding the time-varying parameters of the transmissibility and mortality rates. As a result, we increase the number of hidden states aiming to effectively describe the dynamics of the epidemic evolution. At each time step, the innovation of the hidden parameters is based on log-normal distributions of zero mean, ${\text{dlog}}{\beta _t}\sim N\left( {0,{ }{\sigma _1}\left( t \right)} \right)$ and ${\text{dlog}}{\mu _t}\sim N\left( {0,{ }{\sigma _2}\left( t \right)} \right)$. It is important that according to this decision, we do not force the evolution of ${\beta _t}$ and ${\mu _t}$ to display any predefined trends, while their trajectory is determined completely through the model's likelihood to the observations. Finally, for the respective variances we employ fractions of the parameters' values at the previous time step, aiming to introduce a dynamic parameter evolution.

Algorithm 1 summarizes the implementation of the epidemiologically informed sequential importance resampling particle filter, with penalty factors regarding the evolution of cumulative recovered and deceased cases (EI-PF). We denote by $M$ the number of particles, while the decision on whether the existing particles should be resampled, is based on the effective sample size measure of degeneracy [35]

$$\begin{align}{M_{{\text{eff}}}}\left( t \right) = \frac{M}{{1 + {\text{Va}}{{\text{r}}_{p\left( { \bullet |{y_{1:t}}} \right)}}\left( {w\left( {{x_t}} \right)} \right)}}.\end{align} \tag{ 13 }$$

Algorithm 1. Epidemiologically informed sequential importance resampling particle filter.

Require: $M$, ${f_\theta },{\text{ }}{p_{{y_t}|{x_t}}},{\text{ }}$ $T$, ${M_{{\text{th}}}},{ }{\omega _1},{ }{\omega _2},{ }{\omega _3},{ }{\omega _4}$

Initialize $\left\{ {{\boldsymbol{\tilde x}}_0^{\left( i \right)} = {{\left( {{\boldsymbol{x}}_0^{\left( i \right)},\beta _0^{\left( i \right)},\mu _0^{\left( i \right)}} \right)}^T},\,w_0^{\left( i \right)}} \right\}$

for $t = 1,{\text{ }}2, \ldots ,{\text{ }}T{\text{ }}$ do

/* Step 1: Importance sampling

Sample ${\boldsymbol{\tilde x}}_t^{\left( i \right)}\sim \,{f_{\boldsymbol{\theta }}}({\boldsymbol{\tilde x}}_t^{\left( i \right)}|{\boldsymbol{\tilde x}}_{t - 1}^{\left( i \right)},{{\boldsymbol{y}}_t})$

$w_t^{\left( i \right)}\, \propto \,w_{t - 1}^{\left( i \right)}\frac{{p\left( {{{\boldsymbol{y}}_t}|{\boldsymbol{\tilde x}}_t^{\left( i \right)}} \right)p\left( {{\boldsymbol{\tilde x}}_t^{\left( i \right)}|{\boldsymbol{\tilde x}}_{t - 1}^{\left( i \right)}} \right)}}{{{f_{\boldsymbol{\theta }}}\left( {{\boldsymbol{\tilde x}}_t^{\left( i \right)}|{\boldsymbol{\tilde x}}_{t - 1}^{\left( i \right)},{{\boldsymbol{y}}_t}} \right)}} + {\omega _1}{\left( {1 - {\omega _2}} \right)^t}{\frac{{dD}}{{dt}}^{\left( i \right)}} + {\omega _3}{\left( {1 - {\omega _4}} \right)^t}{\frac{{{\text{d}}R}}{{{\text{d}}t}}^{\left( i \right)}}$

for $i = 1,{\text{ }}2, \ldots ,{\text{ }}M{\text{ }}$ do

Normalization $w_t^{\left( i \right)} = { }\frac{{w_t^{\left( i \right)}}}{{{{\mathop \sum \nolimits}}_{i = 1}^Nw_t^{\left( i \right)}}}$

/* Step 2: Resampling when needed

if ${\hat M_{{\text{eff}}}}\left( t \right) < {M_{{\text{th}}}}$

Resample with replacement new $M$ particles from the existing ones based on the

weight distribution $P\left( {k\left( i \right) = l} \right) = w_t^{\left( l \right)},{\text{ }}l = 1, \ldots ,{\text{ }}M$

$\widetilde {\boldsymbol{x}}_{0:t}^{\left( i \right)} = \widetilde {\boldsymbol{x}}_{0:t}^{k\left( i \right)}$

Reset weights $w_t^{\left( i \right)} = {\text{ }}\frac{1}{M}$

end if

end for

end for

As this quantity cannot be calculated directly, we utilize the estimation [38]

$$\begin{align}{\hat M_{{\text{eff}}}}\left( t \right) = { }\frac{1}{{{{\mathop \sum \nolimits}}_{i = 1}^M{{\left( {w_t^{\left( i \right)}} \right)}^2}}}.\end{align} \tag{ 14 }$$

It is necessary to select a threshold value for the implementation of the resampling step. Usually, this threshold is a percentage of the predetermined number of particles $M.$ In our occasion, we choose ${M_{{\text{th}}}} = 0.75M,$ as a relatively frequent resampling procedure has been considered appropriate to mitigate degeneracy issues [38].

Remark. It should be noted that the attenuation parameters ${\omega _2}$, ${\omega _4}$ ∈ $\left[ {0,1} \right]$ have an extra role, except from providing a more flexible penalty determination. In equation (12), the exponent of factors $1 - {\omega _2}$ and $1 - {\omega _4}$ increases as we move to future time steps. As a result,

$$\begin{align} & w_{t - 1}^{\left( i \right)}\frac{{p\left( {{{\boldsymbol{y}}_t}|{\boldsymbol{x}}_t^{\left( i \right)}} \right)p\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)}} \right)}}{{{f_{\boldsymbol{\theta }}}\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)},{{\boldsymbol{y}}_t}} \right)}} + {\omega _1}{\left( {1 - {\omega _2}} \right)^t}\,{\frac{{dD}}{{dt}}^{\left( i \right)}} + {\omega _3}{\left( {1 - {\omega _4}} \right)^t}\,{\frac{{{\text{d}}R}}{{{\text{d}}t}}^{\left( i \right)}} \nonumber\\ & \quad \to w_{t - 1}^{\left( i \right)}\frac{{p\left( {{{\boldsymbol{y}}_t}|{\boldsymbol{x}}_t^{\left( i \right)}} \right)p\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)}} \right)}}{{{f_{\boldsymbol{\theta }}}\left( {{\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)},{{\boldsymbol{y}}_t}} \right)}},\end{align} \tag{ 15 }$$

for $t \to \infty ,$ while this phenomenon amplifies as ${\omega _2},{\text{ }}{\omega _4} \to 1$. Thus, for a sufficiently large $t$, the weight determination procedure returns exponentially fast to the standard process.

2.3. Dataset

2.4. Model comparison

In order to assess the fitting capacity of all the investigated models, we have utilized the likelihood and the root-mean-squared-error (RMSE). The RMSE serves as a standard metric for assessing the adequacy of the applied models. It has been employed in numerous epidemiological studies employing stochastic estimation methods [25, 41–46]. The RMSE contrasts the daily predictions of the model with the recorded daily observations of new cases and deaths, and it is given by

$$\begin{align}{\text{RMSE}} = \sqrt {\frac{{{{\mathop \sum \nolimits}}_{i = 1}^n{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}}}{n}} ,{ }\end{align} \tag{ 16 }$$

where ${y_i}$ and ${\hat y_i},$ denote the true observations and model's predictions, respectively. We notice that the employment of particle filtering considers observational errors, providing a significant advantage compared to deterministic compartmental methods. By using particle filters, we generate interval and not just point estimates, providing a more comprehensive understanding of the uncertainty/noise linked to epidemic phenomena.

To discover the best parameters for the deterministic SEIRD model, we employed a grid search methodology. We examined the model's ability to fit data by exploring various combinations of parameters, including 100 infection rates and 50 recovery rates. Additionally, we tested 10 different incubation rates and 200 mortality rates. To evaluate the performance of the SEIRD model, we computed the RMSE values using the reported counts for new cases and deaths, utilizing a specific parametric set $\left\{ {\beta ,\,\alpha ,\,\kappa ,\,\mu } \right\}$, for $\beta \in \left[ {0.01,{\text{ }}1} \right],\alpha \in \left[ {\frac{1}{{11}},{\text{ }}\frac{1}{3}} \right],\kappa \in \left[ {\frac{1}{{23}},\frac{1}{3}} \right]$ and $\mu \in \left[ {5 \times {{10}^{ - 6}},{\text{ }}2 \times {{10}^{ - 4}}} \right].$ The best parameter combination provides RMSE values of 392.0306 and 3.2385 for the new cases and deaths due to mpox, corresponding to $\beta = 0.58,{\text{ }}\alpha = \frac{1}{4},{\text{ }}\kappa = \frac{1}{4}$ and $\mu = 2.85 \times {10^{ - 5}}.$ The aforementioned results, accompanied by the RMSE values representing the fitting capacity of the examined stochastic models are displayed in table 3.

For the particle filtering model, we display the chosen parameter initializations in table 2. For the incubation rate, a normally distributed variable with a mean of $1/8$ was utilized, as it is assumed that after the exposure to the virus, the individual would exhibit the first symptoms after approximately 8 days [9]. Infected individuals were assumed to recover after 14–28 days [47]; thus, as an initialization for the recovery rate we selected $\kappa \sim U\left( {\frac{1}{{28}},\frac{1}{{14}}} \right)$. Throughout the analysis, these two parameters are kept constant, as there is no evidence in the existing literature to suggest their variation over time. Consequently, we emphasize that both initializations rely on information obtained from existing research, with the intention of identifying appropriate parameter choices that are supported by the literature.

Table 2. Parameter initializations for the epidemiologically informed particle filter's parameters.

Symbol	Parameters/States	Initializations
$\alpha$	Incubation rate	$N\left( {1/8,0.02} \right)$
$\beta$	Transmission rate	$U\left( {1.8,{\text{ }}2} \right)$
$\kappa$	Recovery rate of infective cases	$U\left( {1/28,1/14} \right)$
$\mu$	Mortality rate of infective cases	$U\left( {2 \times {{10}^{ - 5}},4 \times {{10}^{ - 5}}} \right)$

To our knowledge, no accurate transmission rate had been reported in the literature regarding the recent outbreak, while the mortality rate differs for various regions. Therefore, for ${\beta _0}$ and ${\mu _0}$ we opt for uniform distributions with boundaries determined by maximum likelihood estimations. This decision is made due to the lack of substantial evidence in the literature supporting a definitive parameter selection. In the context of the EI-PF model, these rates are regarded as time-dependent, offering improved fitting and predictive capabilities when compared to the particle filter methodologies described in the subsequent results section. Specifically, the maximum likelihood estimations for the EI-PF model were obtained using ${\beta _0}\sim U\left( {1.8,{\text{ }}2} \right)$ and ${\mu _0}$ ∼ $U\left( {2 \times {{10}^{ - 5}},4 \times {{10}^{ - 5}}} \right)$.

In figure 2 we present the produced scaled log-likelihoods for 225 parameter combinations. The scaling has been implemented by dividing with the absolute value of the maximum generated log-likelihood. This representation shows that the maximum scaled log-likelihood corresponds to the pair of central points $\left( {1.9,{\text{ }}3 \times {{10}^{ - 5}}} \right)$, which is associated with the distributions ${\beta _0}\sim U\left( {1.8,{\text{ }}2} \right)$ and ${\mu _0}$ ∼ $U\left( {2 \times {{10}^{ - 5}},4 \times {{10}^{ - 5}}} \right)$. More interval estimates for the ${\beta _0}$ and ${\mu _0}$ parameters have been tested, although they provided reduced log-likelihood values.

Zoom In Zoom Out Reset image size

In summary, the prior distributions for the rates $\alpha$ and $\kappa$, as shown in table 2, are chosen using values obtained from literature sources [9, 47]. Conversely, for the priors of $\beta$ and $\mu$, we have utilized a grid search method based on the maximum log-likelihood values [36], due to the absence of trustworthy literature information for selecting these values. Finally, for the evolution of the time-varying rates we employ ${\text{dlog}}\beta _t^{\left( i \right)}\sim N\left( {0,{ }0.25\beta _{t - 1}^{\left( i \right)}} \right)$ and ${\text{dlog}}\mu _t^{( i )}\sim N( {0,{ }0.25\mu _{t - 1}^{( i )}} )$, leading to a dynamic parameter estimation. Based on the abovementioned formulas, at each time step, we employ a standard deviation equal to the parameters' value on the previous step. Thus, we lead to a more adaptable parameter estimation process. As we previously mentioned, introducing parameters ${\beta _t}$ and ${\mu _t}$ in the dynamic estimation procedure leads to an augmentation of the state space. To advance the system's states by one step ahead, we draw samples for the importance density ${f_{\boldsymbol{\theta }}}({\boldsymbol{x}}_t^{\left( i \right)}|{\boldsymbol{x}}_{t - 1}^{\left( i \right)}).$ Consequently, a state's value at the time step $t,$ relies on its value on the previous time step, considering that the states' propagation follows a first-order Markov process. Ultimately, incorporating the previous values of the two parameter into the aforementioned standard deviations establishes an one step time dependence, resembling the propagation role that governs the system's states.

The choice of 25% was determined after a similar grid search approach based on the maximum log-likelihood values, similar to the case of the prior distributions of the ${\beta _0}$ and ${\mu _0}$ rates. After the application of the grid search method, we have led to the conclusion that excessively high percentages lead to notably irregular trajectories that adversely affect the state estimation procedure. Conversely, excessively low percentages result in relatively constant rates, making the dynamic estimation process ineffective.

A variety of estimations have appeared regarding the mortality rate of mpox. The World Health Organization reports a historical rate of 11%, and for the recent outbreak, the reported estimate is 3%–6% [48]. In the United States, only 1 and 12 deaths have been observed through 21 September and 21 November 2022 respectively, while the cumulative infected cases were 24 040 and 28 935, indicating a death rate of about 0.004% and 0.04%, respectively. A possible explanation for the discrepancy between the expected and observed mortality rates can be attributed to an underestimation of infectious cases, the specific virus strain of the outbreak, and the various demographic characteristics [49]. Moreover, the necessity for a time-varying mortality rate becomes evident.

All final state estimations are based on the medians and empirical distributions of a total of $M =$ 2000 particles. The utilization of additional particles provided similar results, indicating that selecting 2000 particles is a suitable option. However, it is important to note that higher values increase the computational burden of the filter, while leading to longer execution durations. Table 3 displays the RMSE values for the fitting performance of different models in predicting the SEIRD epidemic over a period of 150 days. The deterministic SEIRD model and four particle filter alternatives were assessed, considering both constant and dynamically estimated transmissibility and mortality rates. We should note that the EI-PF, represents the particle filter with time-dependent transmissibility and mortality rates that is combined with the integration of penalty factors in the weight determination step. With the exception of the EI-PF, all other models lacked penalty factors. We should note that the initializations for the transmissibility and mortality rates of all particle filter alternatives are based on maximum likelihood estimations.

Table 3. RMSE values corresponding to the examination of the fitting efficiency of 5 models regarding new cases and deaths. The last row annotated in bold corresponds to the most efficient modeling methodology.

Model	Cases	Deaths
Deterministic SEIRD	392.0306	3.2385
PF-SEIRD with fixed parameters	414.5950	2.8817
PF-SEIRD with varying ${\beta _t}$ and fixed $\mu$	347.2406	2.8129
PF-SEIRD with varying ${\beta _t}$ and ${\mu _t}$	371.5354	1.7990
EI-PF model	265.7940	1.0494

The final row of the table represents the results for the proposed EI-PF model, which exhibited the most accurate fitting ability. The second-best model was the particle filter with time-varying ${\beta _t}$ and ${\mu _t}$ rates (fourth row), while the next best alternative was the particle filter with fixed $\mu$. Specifically, compared to the latter model, the former showed a 7% increase in RMSE for the new cases, while the latter exhibited a 56.35% increase in RMSE for the deceased cases. In addition, the fitting efficiency of the particle filter with fixed parameters is equivalent to that of the deterministic scheme. The former provides a slightly higher RMSE of 5.76% for the new cases, while the latter shows an increase of 12.38% for the number of deaths.

Furthermore, an additional method of confirming the efficacy of the various particle filters compared in the results section involves evaluating the marginal log-likelihood estimation. In table 4, we present the log-likelihoods for the particle filter with fixed parameters, as well as those for the particle filter with time-varying ${\beta _t}$ and constant $\mu$, the particle filter with time-varying ${\beta _t}$ and ${\mu _t}$, along with the corresponding values for the EI-PF, over a time span of 30, 60, 90, 120, and 150 days of mpox outbreak. Across all five examined periods, the epidemiologically reinforced particle filter consistently exhibits the highest likelihood values. The displayed values validate the effectiveness of the proposed EI-PF, while the worst results are still provided by the particle filter with constant parameters. Figures 2–4 depict the fitting capabilities of the different modeling options analyzed. In relation to the number of deceased cases, we illustrate the evolution of the total number of deaths, to aid in comprehending the importance of penalty factors, which remove unsuitable trajectories regarding the phenomenon under study. Although, we underline that we have employed the number of new deaths in the Bivariate Poisson distribution.

Table 4. Comparison of the marginal log-likelihoods of the 4 tested particle filter alternatives, after 30, 60, 90, 120 and 150 days.

	24 July	23 August	22 September	22 October	21 November
PF-SEIRD with fixed parameters	−1098.366	−3124.187	−4995.891	−6706.178	−9147.702
PF-SEIRD with varying ${\beta _t}$ and fixed $\mu$	−699.253	−1842.930	−3055.894	−4318.325	−5325.463
PF-SEIRD with time-varying ${\beta _t}$ and ${\mu _t}$	−704.247	−1856.259	−3068.571	−4286.329	−5216.908
EI-PF model	−244.794	−507.826	−761.815	−1032.052	−1281.958

Figure 3 shows the fitting capacity of the deterministic SEIRD and the proposed EI-PF. We notice that the proposed model follows efficiently the dynamics of the mpox epidemic, especially after the first 30 days for the reported new cases. Undoubtedly, the algorithm's operation that employs the available information of daily observations utilizing the properties of the likelihood function, enhances significantly the model's ability to adapt the changing dynamics that characterize the studied phenomenon. In contrast, the deterministic SEIRD scheme produced by the abovementioned parametric set, exhibits a significant underestimation of new cases, particularly beyond the first 85 days. Specifically, the deterministic approach predicts zero new cases beyond this timeframe, indicating an early termination of the epidemic outbreak and rendering it impractical for predictive purposes. Additionally, substantial disparities are observed in estimating the total number of deaths, clearly indicating that the conventional deterministic scheme fails to capture the incremental pattern of this variable.

Zoom In Zoom Out Reset image size

Since the proposed EI-PF model provided the best fitting capacity, in figure 4 we also present the final estimations accompanied by the respective 90% confidence intervals. These intervals, highlighted with dark blue color, are generated based on the empirical distribution of the 2000 particles. In figure 5, we present the results of the filter with time-varying ${\beta _t}$ and ${\mu _t},$ but without penalty factors. As we see in the second diagram, corresponding to the number of deceased cases, the estimations of this filter do not follow the increasing monotonic behavior that should accompany the evolution of total deaths. This fact is encountered four times (red circles), around day 54, 70, 76 and 143, where we observe negative jumps in the total number of deaths. On the other hand, this phenomenon is eliminated when we apply the penalty factors (figure 4), while the irregularities that we encounter in the boundaries of the confidence intervals of the second diagram of figure 5, are also satisfactorily smoothed. These observations, in parallel with the RMSE and log-likelihood values of tables 3 and 4, validate the usefulness of the employment of these factors during the weight determination process.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 7, we observe the initial and final ($t = 150\,{\text{d}}$ays) histograms of model's parameters produced by the $M$ particles, providing valuable information about their moments, namely mean, variance, skewness etc. The red plots display the parameters' initial distributions, whose selection was described earlier in the results section (table 2), while the blue ones exhibit the respective final distributions. For parameters $a$ and $\kappa ,$ the final histograms are similar to the initial ones, as these rates are held constant over time, for all particle filter alternatives. On the other hand, the final histogram of ${\beta _t}$ displays a log-normal form—which is quite expected—that is concentrated around lower values compared to the initial distribution. This fact is associated with the attenuation of the examined epidemic wave in the United States. A similar log-normal behavior is shown for the final histogram of ${\mu _t}$, with a much more concentrated density around 0.

Next, we present the evolution of the EI-PF's time-varying parameters for 150 days of mpox (figure 6). The red lines represent the estimations derived from the weighted sum of the filter's particles, accompanied by the corresponding 90% confidence interval. Regarding the evolution of ${\beta _t}$, the filter adjusts fast to lower values, while the most prevalent increase emerge slightly before August where more intensified numbers of new infections occur. Beyond this timeframe, the transmissibility rate gradually decreases, closely tracking the progression of the new cases time series. Moreover, the mortality rate exhibits a relatively stable trend until mid-September, followed by a more fluctuating pattern for the remaining period. The upward spikes are associated with the positive jumps observed, regarding the deceased cases shown in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we extensively investigate the predictive efficiency of the proposed EI-PF using four successive time intervals, regarding new cases and deaths of the mpox outbreak. During this process, we utilize only the transition density provided by the discrete distributions presented in equation (9). For the four following predictive attempts, we fit our model to the data of the first 70, 90, 110 and 140 days, and then we make predictions for the next ten time steps. The time-varying ${\beta _t}$ and ${\mu _t}$ are updated during the fitting process and then they are considered fixed for the prediction period. As a result, predictions are acquired for the periods 3/9/2022–12/9/2022, 23/9/2022–2/10/2022, 13/10/2022–22/10/2022 and 12/11/2022–21/11/2022.

In figure 8, we show the fitting and predictive performance of the EI-PF model. The median estimates obtained through the fitting process are represented by the red lines, while the light blue lines indicate the 10 days predictions. The pink shading highlights each prediction period. Both the fitting and predictive estimations are accompanied by 90% confidence intervals derived from the empirical distribution of the filtering method. For the four examined periods, the RMSE values obtained from the EI-PF model are as follows: 330.3228, 189.0477, 135.6662, and 37.8959 for new cases, and 0, 0, 1.6432, and 1.6733 for deaths, respectively. In comparison, the particle filter without penalty factors yields the following RMSE values: 371.5354, 168.6413, 136.5171, and 40.1995 for new cases, and 1.8140, 1, 1.6432 and 6.4187 for deaths. These results demonstrate the enhanced predictive capability of the EI-PF approach. It is worth noting that the particle filter without penalty factors and with varying ${\beta _t}$ and ${\mu _t}$ was chosen for the above comparison as it constituted the alternative with the second best fitting performance.

Zoom In Zoom Out Reset image size

In the present study, a novel EI-PF is introduced for the modeling of epidemic outbreaks. The presented approach incorporates time-varying transmissibility and mortality rates, which frequently characterize the prevalence of epidemics [19, 20, 25, 50]. Furthermore, two extra parameters were integrated that improve the particle weight distribution and consequently the resampling procedure. By comparing (a) the deterministic model, (b) a particle filtering model with fixed parameters, (c) a particle filter with time-varying ${\beta _t}$, (d) a particle filter with time-varying ${\beta _t}$ and ${\mu _t}$ without penalty factors and (e) the particle filtering model with time-varying ${\beta _t}$ and ${\mu _t}$ and enhanced weighting process, the last model exhibited the best fitting and predictive capacity regarding the mpox data.

Compared to the proposed EI-PF, the deterministic SEIRD approach and the particle filter with fixed parameters produce deteriorated RMSE values with respect to the fitting capacity. The former provides an increase of 47.49% and 208.60% for the new cases and deaths, while the latter produces an increase of 55.98% and 174.60%, respectively. Furthermore, the particle filter with the dynamically estimated ${\beta _t}$ and constant $\mu$ provides 30.64% and 168.05%, while the particle filter with time-varying ${\beta _t}$ and ${\mu _t}$ and no penalty factors generates 39.78% and 71.43% higher RMSE values for the same quantities (table 3). Thus, the utilization of not just the time-varying rates, but also the penalty factors, not only limits particle trajectories that are not compatible with the characteristics of epidemics, but also enhances the model's fitting efficiency both for new infections and deaths.

The effectiveness of the EI-PF model in prediction is extensively examined across four periods, providing estimations for the next 10 days. In comparison to the particle filter with time-varying ${\beta _t}$ and ${\mu _t}$, which serves as the most competitive alternative, the EI-PF demonstrates superior performance, exhibiting lower RMSE values for both deaths and infections in nearly all instances. Consequently, the overall estimations derived from the fitting process (figure 3) and the predictive procedure (figure 7) of the EI-PF approach are highly satisfactory, particularly in capturing the highly irregular behavior of new cases, where several negative jumps representing zero new infections occur. Furthermore, incorporating time-varying transmissibility and mortality parameters is crucial, as the number of new deaths increases when the infection wave attenuates. Models with fixed parameters fail to adapt effectively to this pattern since the onset of most new deceased cases should coincide with the peak of the infection outbreak.

The particle filter approach, even when the parameters of the system are unknown, show significant potential. When applied to emerging communicable illnesses, where there is frequent (e.g. daily) case count reporting but insufficient parameter information, the approach's benefit is likely to be particularly pronounced. The dynamic model's current state and parameter values can be updated adaptively with the aid of the stochastic approach. By regularly updating the model with the latest data, it is feasible to predict the future and then foresee potential trade-offs between public health measures. The particle filtering, by incorporating stochasticity in its subpopulations could be more inclusive for divergent cases or other hidden states that are not accounted for in the initial investigation. Moreover, the time-dependent rates consider the population variability that is associated with demography or different environment conditions.

Aiming to propose a particle filter whose estimations are compatible with the specificities of epidemics, we improve the particle selection procedure by providing additional information derived from the specificities of the studied phenomenon. Since each particle represents a possible trajectory of the epidemic's spread in the population, the elimination of particles presenting fewer deceased or recovered cases as we move to future time steps seems necessary. This behavior derives from the stochasticity that accompanies the resampling step of particle filters. Even if we draw samples from a density which is concentrated on a set where ${D_t} \unicode{x2A7E} {D_{t - 1}}$ and ${R_t} \unicode{x2A7E} {R_{t - 1}}$, the stochasticity of the resampling step may generate more particles that show fewer deceased cases on time step $t$ compared to $t - 1$, leading to imbalances in the final filter's estimation (figure 4, diagram 2). On the other hand, particles that display positive derivatives for the total number of deaths and recoveries with respect to time should be rewarded.

We emphasize that the choice of penalty factors is based on the increasing trend observed in the cumulative number of deaths and recoveries, a standard trend to all epidemic scenarios. As a result, our goal is to present a modification of the particle filter technique applicable to any epidemic situation, regardless of the specific data collected. To ensure the upward trend in recoveries and deaths, the most straightforward approach is integrating the derivatives of these states. Moreover, the rewarding of positive derivatives is proved to be beneficial for our dataset, while it ensures the generation of particle trajectories that satisfy the properties of the examined natural phenomenon. The impact of this rewarding can be easily adjusted through the $\omega$ factors, with respect to the specificities of other epidemic scenarios. Simultaneously, ${\omega _1}$ and ${\omega _3}$, are utilized to balance the influence of derivatives against the expected likelihood levels. Furthermore, ${\omega _2}$ and ${\omega _4}$ provide users with the capability to mitigate the effects of penalty factors beyond the initial timeframe, if required for their particular application.

Since we cannot avoid the implementation of the resampling process, due to its contribution to the confrontation of the degeneracy issue, which is associated with the methodology of particle filters, we culminate in the inclusion of the penalty factors that reduce the possibility of selecting particles that display negative jumps, leading to a smoother estimation process, which is compatible with the laws of the physical phenomenon. We note that this method is quite efficient especially during periods where we observe stability or mild increases in the total recoveries or deaths. Finally, even if the approach is oriented towards the utilization of recoveries and deaths, it presents noteworthy improvements in estimating the number of new cases, too.

In general, we shall note that a stochastic methodology is highly recommended for epidemics, since a model that aims to describe all possible transitions of this event would lead to a highly complex model, while errors in the records are unknown, making the process computationally demanding and the fitting performance unreliable due to overfitting effects. For instance, transitions characterized by minor or negligible rates, like the transitions from the exposed directly to deceased state or an infection rate that derives from the exposed individuals, are integrated in the system's noise. We note that the SEIRD scheme should be considered appropriate for the examination of mpox, since the publicly available datasets only include records for the number of infected and deceased cases, which is a frequently observed phenomenon is other cases of epidemics. Thus, the implementation of a more complex epidemiological structure, can lead to overfitting effects—that accompany both artificial intelligence and statistical models [25, 51, 52]—since any additional states added to the model, are not supported by real-time data.

The proposed method comes with limitations. The assumption of a fixed population with homogenous contacts represents a simplification of the phenomenon, as multiple members could either enter or exit the population during specific time intervals. These assumptions are quite frequent in mathematical modeling of epidemics, although additional information is necessary for the valid calibration of the model. We suggest that the employment of an open system is more suitable in describing long time interval (1 year or above) or when there is reliable information about immigration waves during the studied period. This can be achieved by adopting extra parameters that represent both the inflow of new members in the population due to births or inbound immigrating population and the outflow of cases from the system due to deaths unrelated to the epidemic or outbound migration.

The modeling of contagious disease's evolution plays an important role regarding many aspects of public health and could be used to develop policy guidelines to prevent further transmission. The recent COVID-19 outbreak highlighted the significance of an efficient and timely evaluation of the spread dynamics aiming to reduce the health and socio-economic drawbacks of any epidemic. The recent mpox epidemic, although it had a lower magnitude compared to COVID-19, has alerted medical professionals to quickly evaluate the situation and provide appropriate solutions to suppress the spread.

In the current study, a novel, hybrid particle filtering epidemiological model is proposed, which effectively combines the elements of a deterministic SEIRD model with the inclusion of stochastic parameters, in order to efficiently evaluate the dynamics of the disease. The model is applied to the United States mpox data from 25 June to 21 November 2022, and indicated that particle filtering in combination with dynamic models can be utilized for the efficient evaluation of emerging communicable diseases, especially when the system's parameters are unknown.

Throughout the iterative determination of particle weights, we emphasize the incorporation of additional parameters called penalty factors, which represents one of the main novelties of the present analysis. The necessity of the resampling process to address the degeneracy issue of particle filters, necessitates the use of penalty factors, aiming to eliminate irregular particle trajectories that do not follow the dynamics of epidemics. Finally, the incorporation of these parameters improves the fitting and prediction of new cases and deaths, as it is extensively described in the results section.

The proposed model exhibits resilience when applied to the available data, making it a valuable tool for predicting the progression of emerging infectious diseases and assisting in the development of health policy guidelines. It is particularly useful in situations where real-time observations and system dynamics alone are inadequate for effectively explaining the changing patterns of an epidemic due to uncertainties. This approach can be extended to various existing and future epidemic events, including Ebola, dengue, and hepatitis, as the employed transitions and epidemiological states are broadly applicable. Further research could explore additional stochastic approaches that provide insights into crucial epidemic characteristics [53–56] or investigate the influence of chronic diseases [57] on the prevalence of new epidemic phenomena.

The utilized dataset exists in the 'Our World in Data' Monkeypox GitHub repository and can be accessed through the link https://raw.githubusercontent.com/owid/monkeypox/2088d25b06d19276d07e7182870d771fafbea450/owid-monkeypox-data.csv.

This research was funded by the Special Account for Research Funds of Aristotle University of Thessaloniki (AUTH).

The authors declare that there are no conflicts of interest.