Exploring the dynamics of monkeypox transmission with data-driven methods and a deterministic model

Introduction: Mpox (formerly monkeypox) is an infectious disease that spreads mostly through direct contact with infected animals or people's blood, bodily fluids, or cutaneous or mucosal lesions. In light of the global outbreak that occurred in 2022–2023, in this paper, we analyzed global Mpox univariate time series data and provided a comprehensive analysis of disease outbreaks across the world, including the USA with Brazil and three continents: North America, South America, and Europe. The novelty of this study is that it delved into the Mpox time series data by implementing the data-driven methods and a mathematical model concurrently—an aspect not typically addressed in the existing literature. The study is also important because implementing these models concurrently improved our predictions' reliability for infectious diseases.

Methods: We proposed a traditional compartmental model and also implemented deep learning models (1D- convolutional neural network (CNN), long-short term memory (LSTM), bidirectional LSTM (BiLSTM), hybrid CNN-LSTM, and CNN-BiLSTM) as well as statistical time series models: autoregressive integrated moving average (ARIMA) and exponential smoothing on the Mpox data. We also employed the least squares method fitting to estimate the essential epidemiological parameters in the proposed deterministic model.

Results: The primary finding of the deterministic model is that vaccination rates can flatten the curve of infected dynamics and influence the basic reproduction number. Through the numerical simulations, we determined that increased vaccination among the susceptible human population is crucial to control disease transmission. Moreover, in case of an outbreak, our model showed the potential for epidemic control by adjusting the key epidemiological parameters, namely the baseline contact rate and the proportion of contacts within the human population. Next, we analyzed data-driven models that contribute to a comprehensive understanding of disease dynamics in different locations. Additionally, we trained models to provide short-term (eight-week) predictions across various geographical locations, and all eight models produced reliable results.

Conclusion: This study utilized a comprehensive framework to investigate univariate time series data to understand the dynamics of Mpox transmission. The prediction showed that Mpox is in its die-out situation as of July 29, 2023. Moreover, the deterministic model showed the importance of the Mpox vaccination in mitigating the Mpox transmission and highlighted the significance of effectively adjusting key epidemiological parameters during outbreaks, particularly the contact rate in high-risk groups.

1 Introduction

In the past 20 years, the world has experienced several epidemics that occurred throughout history, as mentioned in the previous studies by Piret and Boivin (1). Although many diseases, such as H1N1 influenza, COVID-19, and the current Mpox virus, have posed a major public health threat (2), mathematical and data-driven modeling can help to control and prevent the spread of disease by understanding the dynamics of disease. Moreover, a model can predict the potential course of disease progression and offer the best and most practical solutions to prevent pandemics when it incorporates fundamental features (3). In the initial stages of an outbreak, researchers paid widespread attention to mathematical modeling, with the general anticipation that deterministic models would accurately forecast the pandemic’s course, as it helps to understand the exponential increase of infections. This resurgence of disease modeling is not new. However, data-driven modeling opens a new era to explore the dynamics of diseases using real-world data. While numerous data-driven models have interpretability issues, focusing solely on a single data-driven approach may mislead decision-makers. In spite of having some limitations, data-driven techniques have increased their popularity in disease modeling for controlling epidemics using real-world data (4), particularly in epidemiological data. On the other hand, the deterministic model offers a structured and mathematically tractable framework for understanding the spread of infectious diseases, helping disease control, and enabling informed public health decisions. Therefore, integrating mathematical models and data-driven approaches becomes imperative for effective public health data implementation and control strategies, allowing for a more comprehensive and adaptive response to dynamic health challenges.

Mpox is a viral zoonotic infection disease caused by the monkeypox virus, and its transmission is considered one of the threats to human health, which may happen due to increased animal-to-human contact (5), human-to-human transmission (6), and climate change as it influences the environment of their vectors (7). More specifically, transmission through environmental factors can occur when a person or an animal touches surfaces or materials recently contaminated with the virus, whether from infected humans or animals (8). Another zoonotic virus related to Mpox has also been found to be transmitted from animals to humans (9). The virus is a member of the genus Orthopoxvirus (10, 11), and it is related to smallpox (12), but it is generally less severe in humans (13). The disease was first identified in 1958 during an outbreak in an animal facility in Copenhagen, Denmark (14). It was later recognized as a cause of human illness in 1970 in the Democratic Republic of Congo (14–16). The natural reservoir of the virus is believed to be rodents, such as squirrels and prairie dogs, but it has also been identified in other animals, including monkeys, porcupines, and Gambian giant rats (13). However, the ongoing Mpox virus is an emergent human pathogen recently studied by Emily and Sassine (17), Khan et al. (18), and John et al. (19). It has been found that the disease spreads from person to person through close contact with respiratory or other bodily fluids (13). The symptoms of Mpox include fever, headache, muscle aches, and a rash that spreads from the face to other parts of the body and eventually forms a scab before falling off (13). As of July 10, 2023, no specific treatment has been approved for Mpox, but the disease can be prevented through vaccination proposed by Andrea et al. (11), avoiding close contact with infected animals or people, and generally practicing good hygiene (13).

Recently, Kanj et al. (20), Hasan and Saeed (21), and Jeta et al. (22) reviewed nineteen diseases, of which nine employed compartments compartmental (deterministic) models to analyze the transmission dynamics of the Mpox virus in populations comprising both humans and non-humans. These studies explored variations of the classical SIR (Susceptible-Infected-Recovered) vector-borne models (23–25) within the context of human and non-human interaction and utilized them to study Mpox disease transmission. Incorporating vaccination classes (26) is essential for the disease’s accurate representation of progression and incubation. Moreover, the lessons from COVID-19 highlight that vaccinating susceptible populations can address the threats of Mpox and other emerging infectious diseases. Thus, we attempted to introduce a Mpox vaccination class and integrate the smallpox vaccine or dose-1 vaccine into the deterministic model proposed by Mesady et al. (27), Yuan et al. (24), and Esteban et al. (25) to provide human protection strategies along with the model fitting for reported Mpox cases across locations.

The rapid spread of Mpox has disrupted the globe when the COVID-19 disease transmission is declining. During this time, only a few articles have been published on detecting Mpox disease using deep learning on the image data (28). Many articles utilize machine learning techniques to forecast COVID-19 (29, 30) time series such as extreme machine learning (ELM), multilayer perceptron (MLP), long short-term memory (LSTM) networks, gated recurrent unit (GRU), convolution neural network (CNN), and deep neural network (DNN) methods on time series data used to forecast COVID-19 cases and predict a possible ending point of the outbreak. Although LSTM is the most popular neural network model for time series prediction, CNNs can also perform this task effectively. Moreover, CNNs have been shown to outperform RNNs in various tasks involving sequential data (31). CNNs can also learn long-term dependencies in time series data through a combination of multiple 1D convolution layers. In 2022, Ketu and Mishra (19) used CNN-LSTM, a hybrid deep-learning prediction model, to forecast the COVID-19 pandemic across India. While some analytical studies have used the Mpox disease transmission models for understanding the dynamics of this disease (2, 27), less attention has been given to predicting the Mpox virus infection using data-driven modeling (28, 32).

In this study, we forecast Mpox disease transmission for short-term predictions across diverse geographical areas by utilizing concurrently the deterministic modeling and advanced deep learning techniques (including 1D-CNN, LSTM, BiLSTM, hybrid CNN-LSTM, and CNN-BiLSTM) alongside statistical time series models like ARIMA and exponential smoothing. We employed a deterministic model and ran simulations to gain valuable insights into the dynamics of the Mpox outbreak and evaluate the impact of vaccination on disease control. Our analysis incorporated real Mpox data from ourworldindata (33), which helped to identify the disease spread across diverse global locations using our epidemiological model. Furthermore, we estimate critical epidemiological parameters through data fitting within the deterministic model. The consistent outcomes from the eight models utilized in this study indicate reliable results. Moreover, our study highlights the significance of time series models and demonstrates that vaccination rates have the potential to “flatten the curve” of infected cases and influence the basic reproduction number. Increased vaccination among susceptible populations is pivotal in effectively managing disease transmission. While the threat of outbreaks persists, reducing contact rates among high-risk groups can mitigate the disease’s impact in specific regions or communities. Importantly, our eight employed models provide an eight-week short-term prediction showing that Mpox is currently in a declining phase. These findings contribute to a comprehensive understanding of disease dynamics across various locations, guiding targeted intervention strategies for controlling and mitigating infectious disease outbreaks.

The remainder of this paper is divided into six sections. Section 1. provides an account of the datasets and summarizes the global and regional Mpox cases. Section 3. briefly describes deterministic and data-driven models. Section 4. implements the univariate Mpox dataset in the time series models and fits the deterministic model to find the epidemiological parameters. Fitting the data is performed using MATLAB and its optimization toolbox, such as fmincon. Section 5. concludes with a discussion, and the conclusion is drawn in Section 6.

2 Materials

According to the Mpox data ourworldindata (33), the current Mpox virus spread globally, surpassing previous outbreaks and raising severe public health concerns in 2022–2023 by spreading to many regions (112 countries) worldwide. Excluding 27 cases reported before May 1, 2022, our study focuses on ourworldindata (33), which spans from May 1, 2022, to May 31, 2023, comprising 87,875 global new confirmed cases with 143 deaths. Moreover, this study employs various methodologies to predict the trajectory of the 2022–2023 Mpox epidemic in different regions, which is later verified with an additional eight weeks of Mpox cases from June 2023 to July 2023. The current transmission of Mpox has been reported on all continents except for Antarctica (Figure 1), as of May 31, 2023. The continents reporting the most cases are mainly North America (37,058), Europe (25,624), and the South Americas (22,357). Like coronavirus disease 2019 (COVID-19), the United States of America still leads the total confirmed cases at 30,225 as of May 31.

Figure 1

Figure 1. Distribution of the global population sizes, Mpox cases, and incidence: Panel (A) displays the population of each continent, where the world population is 7,795,237,030, and panel (B) shows the distribution of total Mpox cases in seven continents by the pie chart. Panel (C) presents incidence per 100 k individuals across the seven continents, where the maximum incidence occurs in North America (6.26 per 100 k), South America (5.19 per 100 k), and Europe (3.43 per 100 k). Panel (D) shows that Spain has a higher incidence of Mpox cases (greater than 15 per 100 k) compared to the other countries. Specifically, USA, Portugal, and Luxembourg maintain an incidence rate of 8 cases per 100 k, while the incidence rate in the remaining countries is below 8 per 100 k.

2.1 Mpox dataset description

The research materials and methods involve data pre-processing and building predictive models for time series data from May 1, 2022, to May 31, 2023. The first phase involves converting the daily dataset into a weekly time series, as shown in Figure 2, normalizing the data using Equation 1, visualizing mpox cases in global and US regions in Figures 3 and 4, respectively, and splitting the data into training and test sets, as depicted in Figure 5. The second phase involves optimizing and training the models using the best hyperparameters on a train set, which starts on week 0 (May 1, 2022) and ends on week 40 (February 11, 2023). Then, the trained models are evaluated on a test set from week 40 (February 11, 2023), to week 56 (May 31, 2023). Figures 2, 5 provide an overview of the entire procedure. Furthermore, a year is conventionally classified as an epidemic year for a given region if the incidence of Mpox exceeds 100 cases per 100,000 individuals in a given year, for example, a year between May 2022 and May 2023. A non-epidemic year when the incidence of Mpox remains below this threshold like Stolerman et al. (34). In this work, we use this threshold to identify which country has the epidemic year for Mpox and then use the model fitting using the deterministic model in some regions. From Figures 1, 4, we see no locations that may be considered epidemic years for Mpox disease, but the District of Columbia in the USA, is close to this threshold. Figure 3 visualizes global confirmed human Mpox cases by country, continent, and geographic region, and it reveals that Mpox is transmitted more in some regions of the world. It also reveals that the majority of cases are reported in the USA.

Figure 2

Figure 2. Preparing time series data: This figure shows five important phases: data pre-processing steps (panels A and B, top left), weekly time series of reported cases (panel C, top right), weekly time series of reported cases on a normalized scale (panel D bottom right), and the decomposition of the confirmed cases to check seasonality and stationarity (panel E, bottom left). Given the ADF Statistic of $-$2.04679, which is negative, and a $p$-value of 0.26644, exceeding the significance threshold of 0.05, we can infer that the EPI weekly worldwide data exhibits non-stationarity. Panel E also demonstrates the seasonality patterns, as the seasonal decomposition shows the rise and fall, with a period of seven. We observe an increasing trend in the number of reported cases until the EPI week of 14 (between July and August, during the summer), then started decreasing. Panel E displays the additive decomposition, where the seasonal chart scales between $-$100 to 100 and the trend chart scales between $-$1,000 to 1,000.

Figure 3

Figure 3. Visualization of globally reported Mpox cases: Panel (A) displays the confirmed Mpox cases across 112 countries using a treemap visualization. Panel (B) shows the global and regional spread of confirmed Mpox cases on a world map. Among all countries, the USA has recorded the highest number of cases (30,225), followed by Brazil (10,941), Spain (7,555), France (4,146), Colombia (4,090), and Mexico (4,017).

Figure 4

Figure 4. Visualization the spread of Mpox cases in the USA: Panel (A) displays the distribution of Mpox cases (left top) across the US map (35) . Panel (B) shows the time series of confirmed Mpox cases (right top) from week 0 to 51, covering the period from June 4, 2022, to May 27, 2023. Panel (C) illustrates the states with the highest incidence of Mpox per 100,000 (100 k) people in the United States. This panel also shows that the District of Columbia (77.01% per 100 k) has higher Mpox incidence from May 2022 to May 2023. Next, the other states such as New York (18.61% per 100 k), California (14.58% per 100 k), Florida (13.44% per 100 k), Maryland (12.1% per 100 k), Illinois (11.48% per 100 k), Nevada (10.04% per 100 k), and Texas (10.33% per 100 k) are the higher risk of an epidemic. Additionally, Washington (9.06% per 100 k), Arizona (8.28% per 100 k), and New Jersey (25.57% per 100 k) have a higher chance of an epidemic compared to the other states.

Figure 5

Figure 5. Data partitioning into training, validation, and testing sets: This figure illustrates the decomposition of the training, validation, and testing data based on the EPI week for deep learning (1D-CNN, LSTM, BiLSTM, hybrid CNN-LSTM, and CNN-BiLSTM), statistical (ARIMA and exponential smoothing), and deterministic models. Note that data from May 01 to May 31, 2023, was accessed on June 03, 2023. Later, new data for comparing model predictions were accessed on August 01, 2023.

2.2 Processing the time series data

Many real-world disease datasets exhibit temporal patterns, representing the number of infected cases at regular intervals, such as daily or weekly, naturally forming time-series data. Figures 2A,B and A1 (appendix) illustrate common issues stemming from insufficient or low reported data on weekends, particularly around Sundays. To mitigate the impact of reduced weekend reporting, we adopted the epidemiological week (EPI week) or CDC week (where each EPI week begins on Sunday and ends on Saturday) as depicted in Figure 2B. By filtering the data using the EPI week, as shown in Figure 2B, we aimed to reduce the noise inherent in the daily data shown in Figures 2A and A2. This resulted in grouping 395 data points into EPI weeks for each region, providing a dataset spanning 57 weeks.

2.2.1 Data normalization

Modeling time series data presents challenges due to its dynamic and nonlinear nature, as highlighted by Tealab et al. (36) in 2017 and McFarland et al. (37) in 2018. Data normalization techniques, such as the Min-Max formula (Equation 1), are commonly used to address this issue. Notably, data normalization plays a crucial role in ensuring the quality of data prior to model training, and as such, it can significantly impact the performance of any subsequent analysis or modeling effort. The Min-Max formula (Equation 1) described by Wibawa et al. (38) in 2022 is used to normalize data in Figure 2D, resulting in smaller intervals within 0–1. Following (36, 37), we define the data normalization equation as

$y_t^{\mathrm{norm}}=\frac{y_t-y_{min}}{y_{max}-y_{min}},(1)$

where $y_t^{\mathrm{norm}}$ is the result of normalization, and $y_t$ is the data to be normalized, while $y_{min}$ and $y_{max}$ stand for the minimum and maximum value of the data. The deterministic model, ARIMA, and Exponential smoothing models are independent of any other variable; hence, this study does not require normalization for these models. The only normalized data was used in the deep learning models such as CNN, LSTM, BiLSTM, and hybrid CNN-LSTM and CNN-BiLSTM, as normalization can help prevent vanishing and exploding gradients during training.

2.2.2 Split time series data

We define the training, testing, and validation datasets. In 2018, Kuhn and Johnson (39) introduced a section titled “Data Splitting Recommendations,” and later in 2021 Nguyen et al. (40) outlined the limitations of using a sole “test set.” Of course, there are no fixed rules for separation training and testing data sets, nor is there an “ideal” ratio for splitting a dataset, but experts commonly used ratios such as 70:30; 80:20; 65:35; 60:40, etc. Given our current small dataset size, this experiment uses the first 60 percent of the data points in the time series as the training set, 10 percent as validation, and the last 30 percent of the data points as the “test set,” as depicted in Figure 5. Additionally, to validate the predictions made by the models, we consider another eight weeks of cases from June 2023 to July 2023 accessed this data on August 01, 2023.

3 Methods

This section briefly outlines several time series prediction models: a deterministic model and various data-driven methods, including deep learning models such as CNN, LSTM, BiLSTM, hybrid CNN-LSTM, and CNN-BiLSTM along with two classical statistical models: ARIMA and exponential smoothing. We aim to demonstrate how to build and choose the optimal predictive models; Figures 6 and 7 show the building blocks of these methods used to analyze Mpox time series data.

Figure 6

Figure 6. Building Blocks of deterministic and data-driven models: (A) displays different methods and demonstrates how to build the predictive models. (B) identifies the best predictive models by evaluating those models.

Figure 7

Figure 7. Building Blocks of data-driven models: This figure was adapted from Crick (41), illustrating how a model optimizer is generally used in model training. The loss function $L(y,\acute{y})$ is computed by comparing predicted values $(\acute{y})$ and true values $(y)$, and model parameters are adjusted to minimize $L(y,\acute{y})$ that enhances the accuracy of prediction model $f$.

3.1 Performance metrics

In this study, we used three performance criteria (42, 43), namely mean absolute error (MAE), root mean square error (RMSE), and normalized RMSE (nRMSE) to evaluate the predictive accuracy of the developed models, where lower values indicate better performance (44). To define these performance metrics, we set $y$ as the actual data value, $\acute{y}$ as the forecast value, and $n$ as the number of observations. Then the mathematical formulas for these error evaluation measures are as follows:

$\text{MAE}=\frac{\sum _{t=1}^n|y-\acute{y}|}{n},\text{RMSE}=\sqrt{\frac{\sum _{t=1}^n(y-\acute{y}{)}^2}{n}},\text{and}\text{nRMSE}=\frac{\text{RMSE}}{(y_{max}-y_{min})}.(2)$

One can also compute nRMSE based on the average $(\overline{y})$ of the actual data values, calculated as $\frac{\text{RMSE}}{\overline{y}}$.

3.2 Deterministic model

We built a model to analyze the transmission of Mpox within and between human and rodent (non-human) populations. Extending the SEIR-like vector-borne models (24, 25), we introduce a deterministic model and then observe how monkeypox spreads through the human population over time by inputting different values for these variables and running simulations. The models are based on the following assumptions: (1) All individuals are identical, that is, age, sex, social status, and other characteristics do not affect the probability of infection; (2) Human population mixing is homogeneous; (3) There is no inherited immunity; (4) With these assumptions, we consider the population sizes $N(t)$ with a per capita birth rate $\mathrm{\Theta }$ in susceptible class and death rate $\mu$ in all classes. Also, $\theta =\mathrm{\Theta }N(t)$ represents the new birth in the susceptible class. Therefore, the rate of change of a population is proportional to the total population size, i.e. (Equation 3),

$\dot{N}(t)=\mathrm{\Theta }N(t)-\mu N(t)=rN(t),\text{where}r=\mathrm{\Theta }-\mu \text{is the population growth rate.}(3)$

We now propose a deterministic model on the transmission dynamics of monkeypox consisting of two groups: the human population and the rodent (non-human) population. The high-risk group will increase transmission in the low-risk group, and the low-risk group will decrease transmission in the high-risk group. These interactions can be represented using the following form (Equation 4):

$\lambda =\left(\begin{array}{c}{\lambda }_h\\ {\lambda }_r\end{array}\right)=\left(\begin{array}{cc}{\beta }_{hh}& {\beta }_{hr}\\ {\beta }_{rh}& {\beta }_{rr}\end{array}\right)\left(\begin{array}{c}{I}_h\\ I_r\end{array}\right)=\left(\begin{array}{c}{\beta }_{hh}{I}_h+{\beta }_{hr}{I}_r\\ {\beta }_{rh}{I}_h+{\beta }_{rr}{I}_r\end{array}\right),(4)$

where the group-specific force of infections, $\lambda =\left(\begin{array}{c}{\lambda }_h\\ {\lambda }_r\end{array}\right)$, the transmission matrix, $T=\left(\begin{array}{cc}{\beta }_{hh}& {\beta }_{hr}\\ {\beta }_{rh}& {\beta }_{rr}\end{array}\right)$, and the infectious class, $I=\left(\begin{array}{c}{I}_h\\ I_r\end{array}\right)$. The diagonal and off-diagonal of $T$ represent within-group and between-group interactions, respectively, and we assume that $\text{det( T) }\ge 0$, as it is reasonable to disease transmission. Additionally, the first component $h$ of ${\beta }_{hr}$ indicates the infected group , and the second component $r$ indicates the infectious group. In this model, each population can infect the other, but the infection moves through the populations separately. To model the spread of Mpox in a metropolitan area, Yuan et al. (24) suggest a transmission risk and control strategy for Mpox using a matrix of contracts among the groups. With that, we assume the transmission rates are greater within groups than between them, and ${\beta }_{hr}={\beta }_{rh}=c_0(1-k_1),{\beta }_{hh}=c_0{k}_1$, and ${\beta }_{rr}=c_0{k}_2$, where $c_0$ the baseline contact rate among the overall population and $k_1$ and $k_2$ are the proportion of contacts within the human and rodent (non-human) overall contacts, respectively. Moreover, the current outbreak of Mpox suggests that sustained human-to-human transmission is quite feasible, and smallpox vaccination protects against Mpox. As of June 2023, similar to the smallpox vaccination (13), Mpox vaccines such as (45) JYNNEOS (46) and ACAM2000 (47) are recommended for individuals at higher risk during the ongoing Mpox outbreak (45). Therefore, incorporating the vaccination class $V$ with a rate $\varphi$ of susceptible populations into the models Peter et al. (48, 49) is required for the Mpox model, which has been proven effective in preventing Mpox infections. Moreover, the smallpox vaccine and dose-1 Mpox vaccination may not offer complete protection against Mpox; therefore, $\tau$ converts the vaccinated class back to a susceptible class. Additionally, besides the two groups, the human population is further subdivided into seven compartments: susceptible humans $S_h(t)$, exposed humans $E_h(t)$, infected humans $I_h(t)$, hospitalized humans $H(t)$, smallpox or dose-1 Mpox vaccinated humans $V_1(t)$, dose-2 vaccinated humans $V_2(t)$, and recovered humans $R_h(t)$. The rodent population is also subdivided into three compartments: susceptible rodents $S_r(t)$, exposed rodents $E_r(t)$, infected rodents $I_r(t)$, and recovered rodents $R_r(t)$ (Table 1).

Table 1

Table 1. Model’s parameters of Figure A2.

Employing the parameters listed in Table 1 and referring to the compartmental diagram (Figure A2) provided in Appendix A, we can derive the following system of differential equations:

$\begin{array}{r}\begin{array}{c}\mathbf{G}\mathbf{r}\mathbf{o}\mathbf{u}\mathbf{p}\mathbf{A}:\mathbf{H}\mathbf{u}\mathbf{m}\mathbf{a}\mathbf{n}\mathbf{p}\mathbf{o}\mathbf{p}\mathbf{u}\mathbf{l}\mathbf{a}\mathbf{t}\mathbf{i}\mathbf{o}\mathbf{n}\\ & {\dot{S}}_h={\theta }_h-{\lambda }_h\frac{S_h}{N_h}-({\mu }_h+\varphi )S_h+\tau V_1\\ & {\dot{E}}_h={\lambda }_h\frac{S_h}{N_h}-({\mu }_h+{\alpha }_h)E_h\\ & {\dot{I}}_h={\alpha }_h{E}_h-[(1-h){\gamma }_1+h{\gamma }_h+({\mu }_h+{\delta }_h)]I_h\\ & \dot{H}=h{\gamma }_h{I}_h-(d+{\mu }_h+{\delta }_h)H\\ & {\dot{V}}_1=\varphi S_h-(\tau +{\mu }_h+\sigma )V_1\\ & {\dot{V}}_2=\sigma V_1-({\mu }_h+{\sigma }_2)V_2\\ & {\dot{R}}_h={\gamma }_2{V}_2+(1-h){\gamma }_1{I}_h+dH-{\mu }_h{R}_h\\ \mathbf{G}\mathbf{r}\mathbf{o}\mathbf{u}\mathbf{p}\mathbf{B}:\mathbf{R}\mathbf{o}\mathbf{d}\mathbf{e}\mathbf{n}\mathbf{t}\mathbf{p}\mathbf{o}\mathbf{p}\mathbf{u}\mathbf{l}\mathbf{a}\mathbf{t}\mathbf{i}\mathbf{o}\mathbf{n}\\ & {\dot{S}}_r={\theta }_r-{\lambda }_r\frac{S_r}{N_r}-{\mu }_r{S}_r\\ & {\dot{E}}_r={\lambda }_r\frac{S_r}{N_r}-({\mu }_r+{\alpha }_r)E_r\\ & {\dot{I}}_r={\alpha }_r{E}_r-({\mu }_r+{\delta }_r+{\gamma }_r)I_r\\ & {\dot{R}}_r={\gamma }_r{I}_r-{\mu }_r{R}_r\end{array}\end{array}(5)$

with the initial conditions in a biologically feasible region ${\mathrm{\Gamma }}_h\times {\mathrm{\Gamma }}_r$, where ${\mathrm{\Gamma }}_h\mathrm{and}{\mathrm{\Gamma }}_r$ are defined by the Equation 6 below:

$\begin{array}{cc}& {\mathit{\Gamma }}_h=\{(S_h,E_h,I_h,H,V_1,V_2,R_h)\in {\mathcal{R}}_{+}^7:S_h+E_h+I_h+H+V_1+V_2+R_h\le N_h\}\text{and}\\ & {\mathit{\Gamma }}_r=\{(S_r,E_r,I_r,R_r)\in {\mathcal{R}}_{+}^4:S_r+E_r+I_r+R_r\le N_r\}.\end{array}(6)$

3.2.1 Model analysis

By summing the first seven equations of group A in the (Equation 5), we have ${\dot{S}}_h+{\dot{E}}_h+{\dot{I}}_h+\dot{H}+{\dot{V}}_1+{\dot{V}}_2+{\dot{R}}_h={\theta }_h-{\mu }_h(S_h+E_h+I_h+H+$ $V_1+V_2+R_h)-{\delta }_h(I_h+H)={\theta }_h-{\mu }_h{N}_h-{\delta }_h(I_h+H)$. Using ${\theta }_h={\mathrm{\Theta }}_h{N}_h$, the differential equation for the human population $N_h$ is given as follows:

${\dot{N}}_h=({\mathrm{\Theta }}_h-{\mu }_h)N_h-{\delta }_h(I_h+H)(7)$

Moreover, Figure 12 illustrates the significance of vaccination for the US population. Similarly, for the non-human population, $N_r$, we have ${\dot{S}}_r+{\dot{E}}_r+{\dot{I}}_r+{\dot{R}}_r={\theta }_r-{\mu }_r(S_r+E_r+I_r+R_r)-{\delta }_r{I}_r={\theta }_r-{\mu }_r{N}_r-{\delta }_r{I}_r$, and using ${\theta }_r={\mathrm{\Theta }}_r{N}_r$ the corresponding differential is given as follows:

${\dot{N}}_r=({\mathrm{\Theta }}_r-{\mu }_r)N_r-{\delta }_r{I}_r.(8)$

Here, $N_h$ and $N_r$ are sums of state variables for human and rodent populations, respectively, which are not constants in general. However, they are constants when ${\mathrm{\Theta }}_h={\mu }_h,{\delta }_h=0$, and ${\mathrm{\Theta }}_r={\mu }_r,{\delta }_r=0$ as ${\dot{N}}_h=0$, ${\dot{N}}_r=0$, respectively (Equations 7–8).

3.2.2 Equilibrium state

There are two types of equilibrium points for the model given in (Equation 5). One is the Disease-Free Equilibrium(DFE) or Mpox-Free equilibrium point, ${ϵ}^0=(S_h^0,E_h^0,I_h^0,H^0,V_1^0,V_2^0,R_h^0,S_r^0,E_r^0,I_r^0,R_r^0)={\displaystyle (\frac{{\theta }_h({\mu }_h+\sigma +\tau )}{{\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))},0,0,0,{\displaystyle \frac{{\theta }_h\varphi }{{\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma )},{\displaystyle \frac{{\theta }_h\sigma \varphi }{({\gamma }_2+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))},\frac{{\gamma }_2{\theta }_h\sigma \varphi }{{\mu }_h({\gamma }_2+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))},\frac{{\theta }_r}{{\mu }_r},0,0,0)}}}$, which corresponds to the absence of infectious individuals. The second equilibrium point is called endemic equilibrium, represented by $(S_h^{\ast },E_h^{\ast },I_h^{\ast },H^{\ast },V_1^{\ast },V_2^{\ast },R_h^{\ast },S_r^{\ast },E_r^{\ast },I_r^{\ast })$, which can be obtained by equating the right-hand side of the system (Equation 5) with zero. We use the symbolic packages of the software Mathematica for the analytic computations in our deterministic model. Using the notation $T_1={\mu }_h+{\alpha }_h,T_2=(1-h){\gamma }_1,T_3={\mu }_h+{\delta }_h,T_4=\tau +{\mu }_h+\sigma ,T_5={\mu }_r+{\alpha }_r,T_6={\mu }_r+{\delta }_r$, and $T_7={\mu }_h+\sigma$ in the system (Equation 5), we obtain the endemic equilibrium point as follows:

$\begin{array}{rl}& S_h^{\ast }=\frac{N_h{T}_4{\theta }_h}{T_4({\lambda }_h+N_h({\mu }_h+\varphi ))-N_h\tau \varphi },E_h^{\ast }=\frac{T_4{\theta }_h{\lambda }_h}{T_1{T}_4({\lambda }_h+N_h{\mu }_h)+N_h{T}_1(T_4-\tau )\varphi }\\ & I_h^{\ast }=\frac{T_4{\alpha }_h{\theta }_h{\lambda }_h}{T_1(h\ast {\gamma }_h+T_2+T_3)(T_4({\lambda }_h+N_h({\mu }_h+\varphi ))-N_h\tau \varphi )},H^{\ast }=\frac{h{\gamma }_h{T}_4{\alpha }_h{\theta }_h{\lambda }_h}{T_1(d+T_3)(h{\gamma }_h+T_2+T_3)(T_4({\lambda }_h+N_h({\mu }_h+\varphi ))-N_h\tau \varphi )}\\ & V_1^{\ast }=\frac{N_h{\theta }_h\varphi }{T_4({\lambda }_h+N_h({\mu }_h+\varphi ))-N_h\tau \varphi },V_2^{\ast }=\frac{N_h{\theta }_h\sigma \varphi }{({\gamma }_2+{\mu }_h)(T_4({\lambda }_h+N_h({\mu }_h+\varphi ))-N_h\tau \varphi )}\\ & R_h^{\ast }=\frac{(d(h{\gamma }_h+T_2)+T_2{T}_3)T_4{\alpha }_h{\theta }_h{\lambda }_h({\gamma }_2+{\mu }_h)+N_h{T}_1(d+T_3)(h{\gamma }_h+T_2+T_3){\gamma }_2{\theta }_h\sigma \varphi }{T_1(d+T_3)(h{\gamma }_h+T_2+T_3){\mu }_h({\gamma }_2+{\mu }_h)(T_4({\lambda }_h+N_h({\mu }_h+\varphi ))-N_h\tau \varphi )}\\ & S_r^{\ast }=\frac{N_r{\theta }_r}{{\lambda }_r+N_r{\mu }_r},E_r^{\ast }=\frac{{\theta }_r{\lambda }_r}{T_5{\lambda }_r+N_r{T}_5{\mu }_r},I_r^{\ast }=\frac{{\alpha }_r{\theta }_r{\lambda }_r}{T_5(T_6+{\gamma }_r)({\lambda }_r+N_r{\mu }_r)},R_r^{\ast }=\frac{{\alpha }_r{\gamma }_r{\theta }_r{\lambda }_r}{T_5(T_6+{\gamma }_r)({\lambda }_r+N_r{\mu }_r)}.\end{array}$

3.2.3 Basic reproduction number

In our proposed Mpox model (Equation 5), compartments $S_h,V_1,V_2,R_h$, and $S_r$ represent the disease-free states , while compartments $E_h,I_h,H,E_r$, and $I_r$ represent the disease class. To investigate whether an infected individual will promote an epidemic outbreak in a susceptible population, we analyze the number of secondary infections produced by an infected individual during the transmission period, called the basic reproduction number $R_0$, at the DFE ${ϵ}^0$. Using the technique (50, 51), the next-generation matrix (NGM), $\kappa =F{V}^{-1}$ at the DFE is defined as

$\kappa =\left(\begin{array}{ccccc}{\displaystyle \frac{{\alpha }_h{\beta }_{hh}{\theta }_h{T}_4}{N_h{T}_1(h{\gamma }_h+T_2+T_3)({\mu }_h{T}_4+\varphi T_7)}}& {\displaystyle \frac{{\beta }_{hh}{\theta }_h{T}_4}{N_h(h{\gamma }_h+T_2+T_3)(muh\text{T4}+\text{T7}\varphi )}}& 0& {\displaystyle \frac{{\alpha }_r{\beta }_{hr}{\theta }_h{T}_4}{N_h{T}_5{T}_6({\mu }_h{T}_4+\varphi T_7)}}& {\displaystyle \frac{{\beta }_{hr}{\theta }_h{T}_4}{N_h{T}_6({\mu }_h{T}_4+\varphi T_7)}}\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ {\displaystyle \frac{{\alpha }_h{\beta }_{rh}{\theta }_r}{{\mu }_r{N}_r{T}_1(h{\gamma }_h+T_2+T_3)}}& {\displaystyle \frac{{\beta }_{rh}{\theta }_r}{{\mu }_r{N}_r(h{\gamma }_h+T_2+T_3)}}& 0& {\displaystyle \frac{{\alpha }_r{\beta }_{rr}{\theta }_r}{{\mu }_r{N}_r{T}_5{T}_6}}& {\displaystyle \frac{{\beta }_{rr}{\theta }_r}{{\mu }_r{N}_r{T}_6}}\\ 0& 0& 0& 0& 0\end{array}\right),(9)$

where the notations $T_1={\mu }_h+{\alpha }_h,T_2=(1-h){\gamma }_1,T_3={\mu }_h+{\delta }_h,T_4=\tau +{\mu }_h+\sigma ,T_5={\mu }_r+{\alpha }_r,T_6={\mu }_r+{\delta }_r$, and $T_7={\mu }_h+\sigma$ are used in the matrices

$\begin{array}{rl}& F=\left(\begin{array}{ccccc}0& {\displaystyle \frac{\beta \text{hh}\theta \text{h}(\mu \text{h}+\sigma +\tau )}{\text{Nh}(\mu \text{h}(\mu \text{h}+\sigma +\tau )+\varphi (\mu \text{h}+\sigma ))}}& 0& 0& {\displaystyle \frac{\beta \text{hr}\theta \text{h}(\mu \text{h}+\sigma +\tau )}{\text{Nh}(\mu \text{h}(\mu \text{h}+\sigma +\tau )+\varphi (\mu \text{h}+\sigma ))}}\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& {\displaystyle \frac{\beta \text{rh}\theta \text{r}}{\mu \text{r}\text{Nr}}}& 0& 0& {\displaystyle \frac{\beta \text{rr}\theta \text{r}}{\mu \text{r}\text{Nr}}}\\ 0& 0& 0& 0& 0\end{array}\right),\\ & V=\left(\begin{array}{ccccc}\alpha \text{h}+\mu \text{h}& 0& 0& 0& 0\\ -{\alpha }_h& {\delta }_h+{\gamma }_1(1-h)+h{\gamma }_h+{\mu }_h& 0& 0& 0\\ 0& -h{\gamma }_h& d+{\delta }_h+{\mu }_h& 0& 0\\ 0& 0& 0& {\alpha }_r+{\mu }_r& 0\\ 0& 0& 0& -{\alpha }_r& {\delta }_r+{\mu }_r\end{array}\right),\text{and}\\ & V^{-1}=\left(\begin{array}{ccccc}{\displaystyle \frac{1}{{\alpha }_h+{\mu }_h}}& 0& 0& 0& 0\\ {\displaystyle \frac{{\alpha }_h}{({\alpha }_h+{\mu }_h)({\delta }_h+{\gamma }_1(1-h)+h{\gamma }_h+{\mu }_h)}}& {\displaystyle \frac{1}{{\delta }_h+{\gamma }_1(1-h)+h{\gamma }_h+{\mu }_h}}& 0& 0& 0\\ {\displaystyle \frac{h{\alpha }_h{\gamma }_h}{({\alpha }_h+{\mu }_h)(d+{\delta }_h+{\mu }_h)({\delta }_h+{\gamma }_1(1-h)+h{\gamma }_h+{\mu }_h)}}& {\displaystyle \frac{h{\gamma }_h}{((d+{\delta }_h+{\mu }_h)({\delta }_h+{\gamma }_1(1-h)+h{\gamma }_h+{\mu }_h)}}& {\displaystyle \frac{1}{d+{\delta }_h+{\mu }_h}}& 0& 0\\ 0& 0& 0& {\displaystyle \frac{1}{{\alpha }_r+{\mu }_r}}& 0\\ 0& 0& 0& {\displaystyle \frac{{\alpha }_r}{({\alpha }_r+{\mu }_r)({\delta }_r+{\mu }_r)}}& \frac{1}{{\delta }_r+{\mu }_r}\end{array}\right).\end{array}$

The NGM (Equation 9) is a non-negative matrix, and, as such, we expect that there will be a single, unique eigenvalue that is positive, real, and strictly greater than all the others. The corresponding characteristic equation of Equation 9 is $det(\kappa -\lambda I)=0$, and this gives the eigenvalues (Equation 10)

${\lambda }_1=0,0,0,{\lambda }_2=\frac{B_1-\sqrt{B_2+B_3^2}}{B_4},{\lambda }_3=\frac{B_1+\sqrt{B_2+B_3^2}}{B_4},(10)$

where, $B_1:=({\alpha }_r{\beta }_{rr}{\theta }_r{N}_h{T}_1(h{\gamma }_h+T_2+T_3)(\varphi ({\mu }_h+\sigma )+{\mu }_h{T}_4)+{\alpha }_h{\beta }_{hh}{\theta }_h{\mu }_r{N}_r{T}_4{T}_5{T}_6$,

$B_2=4{\alpha }_h{\alpha }_r{\theta }_h{\theta }_r{\mu }_r{N}_h{N}_r{T}_1{T}_4{T}_5{T}_6(h{\gamma }_h+T_2+T_3)({\beta }_{hr}{\beta }_{rh}-{\beta }_{hh}{\beta }_{rr})(\varphi ({\mu }_h+\sigma )+{\mu }_h{T}_4)$,

$B_3=({\alpha }_h{\beta }_{rr}{\theta }_r{N}_h{T}_1({\gamma }_{h}h+T_2+T_3)(\varphi ({\mu }_h+\sigma )+{\mu }_h{T}_4)+{\alpha }_h{\beta }_{hh}{\theta }_h{\mu }_r{N}_r{T}_4{T}_5{T}_6)$, and

$B_4=(2{\mu }_r{N}_h{N}_r{T}_1{T}_5{T}_6(h{\gamma }_h+T_2+T_3)(\varphi ({\mu }_h+\sigma )+{\mu }_h{T}_4)$.

Since $det(T)\ge 0$, $B_2+B_3^2$ is well defined and positive, that implies ${\lambda }_3\ge {\lambda }_2$. Now, the basic reproduction number is defined as the largest eigenvalue (spectral radius) of the NGM (Equation 9), and it can be obtained as follows:

${\mathcal{R}}_0=\frac{B_1+\sqrt{B_2+B_3^2}}{B_4}.(11)$

Remark 1.From an epidemiological perspective, an epidemic outbreak will occur if and only if $R_0>1$. Moreover if $det(T)=0$, then $B_2=0$, and the new basic reproduction number becomes

${\mathcal{R}}_0=\frac{B_1+B_3}{B_4}=\frac{{\alpha }_h{\beta }_{hh}{\theta }_h{T}_4}{N_h{T}_1({\gamma }_{h}h+T_2+T_3)(\varphi ({\mu }_h+\sigma )+{\mu }_h{T}_4)}+\frac{{\alpha }_r{\beta }_{rr}{\theta }_r}{{\mu }_r{N}_r{T}_5{T}_6}.(12)$

By substituting the values of $T_i$ for $i=1,\dots ,6$, ${\theta }_h={\mathrm{\Theta }}_h{N}_h$, and ${\theta }_r={\mathrm{\Theta }}_r{N}_r$ into (Equation 12), we derive the resulting

${\mathcal{R}}_0=\frac{{\alpha }_h{\beta }_{hh}{\mathrm{\Theta }}_h({\mu }_h+\sigma +\tau )}{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\delta }_h+{\gamma }_1(1-h)+{\gamma }_{h}h+{\mu }_h)}+\frac{{\alpha }_r{\beta }_{rr}{\mathrm{\Theta }}_r}{{\mu }_r({\alpha }_r+{\mu }_r)({\delta }_r+{\mu }_r)}.(13)$

Remark 2.Increasing the vaccination rate $\varphi$ among susceptible individuals $S_h$ leads to an increase in $B_4$. This increase in vaccination is particularly significant as it results in a reduction of ${\mathcal{R}}_0$, which plays a pivotal role in effectively controlling the transmission of the disease.

3.2.4 Sensitivity analysis of basic reproduction ${\mathcal{R}}_0$

As epidemiological models often involve estimated or fitted parameters, the inherent uncertainty lies in the values of these parameters when making conclusions about the underlying epidemic. Therefore, we carried out a sensitivity analysis of this model parameter concerning the basic reproduction number. The normalized sensitivity index of ${\mathcal{R}}_0$ used in Ngungu et al. (52) and Samuel et al. (53) with respect to parameter $p$ is given by

$S_p^{{\mathcal{R}}_0}=\frac{p}{{\mathcal{R}}_0}\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }p}.(14)$

When $S_p^{{\mathcal{R}}_0}$ is positive for the parameter $p$, it indicates an increase in ${\mathcal{R}}_0$, while a negative value of $S_p^{{\mathcal{R}}_0}$ for the parameter $p$ suggests a decrease in ${\mathcal{R}}_0$. Due to complexity of the actual ${\mathcal{R}}_0$ (Equation 11), we consider the normalized sensitivity index (Equation 14) of ${\mathcal{R}}_0$ (Equation 13) with respect to the parameter $p=({\beta }_{hh},{\alpha }_h,{\mathrm{\Theta }}_h,{\mu }_h,{\delta }_h,h,{\gamma }_h,{\gamma }_1,\sigma ,\tau ,\varphi ,{\beta }_{rr},{\alpha }_r,{\mathrm{\Theta }}_r,{\mu }_r,{\delta }_r)$ as follows:

$\begin{array}{rl}& S_{{\alpha }_h}^{{\mathcal{R}}_0}=\frac{{\alpha }_h}{{\mathcal{R}}_0}\frac{{\beta }_{hh}{\mathrm{\Theta }}_h({\mu }_h+\sigma +\tau )}{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\gamma }_1+{\delta }_h-h{\gamma }_1+{\gamma }_{h}h+{\mu }_h)}\\ & S_{{\beta }_{hh}}^{{\mathcal{R}}_0}=\frac{{\beta }_{hh}}{{\mathcal{R}}_0}\frac{{\alpha }_h{\mathrm{\Theta }}_h({\mu }_h+\sigma +\tau )}{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\gamma }_1+{\delta }_h-h{\gamma }_1+{\gamma }_{h}h+{\mu }_h)}\\ & S_{{\mathrm{\Theta }}_h}^{{\mathcal{R}}_0}=\frac{{\mathrm{\Theta }}_h}{{\mathcal{R}}_0}\frac{{\alpha }_h{\beta }_{hh}({\mu }_h+\sigma +\tau )}{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\gamma }_1+{\delta }_h-h{\gamma }_1+{\gamma }_{h}h+{\mu }_h)}\\ & S_{{\mu }_h}^{{\mathcal{R}}_0}=\frac{{\mu }_h}{{\mathcal{R}}_0}\frac{A}{({\alpha }_h+{\mu }_h{)}^2({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ){)}^2({\gamma }_1+{\delta }_h-h{\gamma }_1+h{\gamma }_h+{\mu }_h{)}^2}\\ & S_{{\delta }_h}^{{\mathcal{R}}_0}=-\frac{{\delta }_h}{{\mathcal{R}}_0}\frac{{\alpha }_h{\beta }_{hh}{\mathrm{\Theta }}_h({\mu }_h+\sigma +\tau )}{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\gamma }_1+{\delta }_h-h{\gamma }_1+{\gamma }_h+{\mu }_h{)}^2}\\ & S_h^{{\mathcal{R}}_0}=\frac{h}{{\mathcal{R}}_0}\frac{{\alpha }_h{\beta }_{hh}{\theta }_h({\gamma }_1-{\gamma }_h)({\mu }_h+\sigma +\tau )}{N_h({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\gamma }_1+{\delta }_h-h{\gamma }_1+h{\gamma }_h+{\mu }_h{)}^2}\\ & S_{{\gamma }_h}^{{\mathcal{R}}_0}=-\frac{{\gamma }_h}{{\mathcal{R}}_0}\frac{{\alpha }_h{\beta }_{hh}h{\mathrm{\Theta }}_h({\mu }_h+\sigma +\tau )}{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\gamma }_1+{\delta }_h-h{\gamma }_1+h{\gamma }_h+{\mu }_h{)}^2}\\ & S_{{\gamma }_1}^{{\mathcal{R}}_0}=-\frac{{\gamma }_1}{{\mathcal{R}}_0}\frac{{\alpha }_h{\beta }_{hh}(1-h){\mathrm{\Theta }}_h({\mu }_h+\sigma +\tau )}{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\gamma }_1+{\delta }_h-h{\gamma }_1+h{\gamma }_h+{\mu }_h{)}^2}\\ & S_{\sigma }^{{\mathcal{R}}_0}=-\frac{{\delta }_h}{{\mathcal{R}}_0}\frac{{\alpha }_h{\beta }_{hh}{\mathrm{\Theta }}_h\tau \varphi }{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ){)}^2({\gamma }_1+{\delta }_h-h{\gamma }_1+h{\gamma }_h+{\mu }_h)}\\ & S_{\tau }^{{\mathcal{R}}_0}=\frac{\tau }{{\mathcal{R}}_0}\frac{{\alpha }_h{\beta }_{hh}{\mathrm{\Theta }}_h\varphi ({\mu }_h+\sigma )}{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ){)}^2({\gamma }_1+{\delta }_h-h{\gamma }_1+h{\gamma }_h+{\mu }_h)}\\ & S_{\varphi }^{{\mathcal{R}}_0}=-\frac{\varphi }{{\mathcal{R}}_0}\frac{{\alpha }_h{\beta }_{hh}{\mathrm{\Theta }}_h({\mu }_h+\sigma )({\mu }_h+\sigma +\tau )}{({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ){)}^2({\gamma }_1+{\delta }_h-h{\gamma }_1+h{\gamma }_h+{\mu }_h)}\\ & S_{{\beta }_{rr}}^{{\mathcal{R}}_0}=\frac{{\beta }_{rr}}{{\mathcal{R}}_0}\frac{{\alpha }_r{\mathrm{\Theta }}_r}{{\mu }_r({\alpha }_r+{\mu }_r)({\delta }_r+{\mu }_r)}\\ & S_{{\alpha }_r}^{{\mathcal{R}}_0}=\frac{{\alpha }_r}{{\mathcal{R}}_0}\frac{{\beta }_{rr}{\theta }_r}{N_r({\alpha }_r+{\mu }_r{)}^2({\delta }_r+{\mu }_r)}\\ & S_{{\mathrm{\Theta }}_r}^{{\mathcal{R}}_0}=\frac{{\mathrm{\Theta }}_r}{{\mathcal{R}}_0}\frac{{\alpha }_r{\beta }_{rr}}{{\mu }_r({\alpha }_r+{\mu }_r)({\delta }_r+{\mu }_r)}\\ & S_{{\mu }_r}^{{\mathcal{R}}_0}=-\frac{{\mu }_r}{{\mathcal{R}}_0}\frac{{\alpha }_r{\beta }_{rr}{\mathrm{\Theta }}_r(2{\mu }_r({\alpha }_r+{\delta }_r)+{\alpha }_r{\delta }_r+3{\mu }_r^2)}{{\mu }_r^2({\alpha }_r+\mu -r{)}^2({\delta }_r+{\mu }_r{)}^2}\\ & S_{{\delta }_r}^{{\mathcal{R}}_0}=-\frac{{\delta }_r}{{\mathcal{R}}_0}\frac{{\alpha }_r{\beta }_{rr}{\mathrm{\Theta }}_r}{{\mu }_r({\alpha }_r+{\mu }_r)({\delta }_r+{\mu }_r{)}^2}\\ & A={\alpha }_h{\beta }_{hh}{\mathrm{\Theta }}_h(-({\alpha }_h+{\mu }_h)({\mu }_h+\sigma +\tau )({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))\\ & -({\alpha }_h+{\mu }_h)({\mu }_h+\sigma +\tau )(2{\mu }_h+\sigma +\tau +\varphi )({\gamma }_1+{\delta }_h-h{\gamma }_1+{\gamma }_{h}h+{\mu }_h)\\ & +({\alpha }_h+{\mu }_h)({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\gamma }_1+{\delta }_h-h{\gamma }_1+h{\gamma }_h+{\mu }_h)\\ & -({\mu }_h+\sigma +\tau )({\mu }_h({\mu }_h+\sigma +\tau )+\varphi ({\mu }_h+\sigma ))({\gamma }_1+{\delta }_h-h{\gamma }_1+h{\gamma }_h+{\mu }_h))\end{array}$

Figure 8

Figure 8. Mpox model sensitivities to its associated parameters for six different locations: Panel (A) and (B) illustrate the sensitivity index of ${\mathcal{R}}_0$ (Equation 13) for $\tau =0$ and $\tau =0.01$, respectively. These are computed with the fitted parameters described in Table 3, variable parameters listed in Table 4, and remaining parameters given in Table 2. The figure demonstrates that the sensitivity index of the parameters ${\beta }_{hh},{\theta }_h,{\gamma }_1,{\mu }_h,\tau$, and $\varphi$ changes when the parameter values change, while the remaining parameters remain unchanged. In other words, the parameters ${\alpha }_h,{\gamma }_h,{\delta }_h,\sigma$, and $h$ are almost constant and do not depend on another parameter, but the remaining parameters show the dependence on a second parameter, affecting the varying sensitivity index. The simulation depicts the sensitivity index of ${\mathcal{R}}_0$ (Equation 13) for six different locations using six different colors and sizes of circles.

Table 2

Table 2. Description of the parameters with the values used in the deterministic modeling of Mpox transmission.

Table 3

Table 3. Population size and fitted parameters in deterministic model fitting.

Table 4

Table 4. Variable parameters and errors in deterministic model fitting.

For our Mpox model (Equations 5 and 6), we analyze the sensitivity index to its associated parameter for six locations to find the most sensitivity parameters in Figure 8.

3.3 Data-driven models

We have now introduced two statistical time series models (ARIMA and exponential smoothing) and five deep learning models (Figure 9). The primary objective of these models is to analyze Mpox data/patterns and make predictions.

3.3.1 ARIMA model

ARIMA is a widely used time series model that is suitable for all kinds of data, including changing trends, seasonality, periodic changes, and random disturbances, i.e., it deals with non-stationary time series. Suppose $y_t$ is the actual observation and $d$th difference is $Y_t={\mathrm{\Delta }}^d{y}_t$. ARIMA is represented as $\text{ARIMA}(p,d,q)$, where $p$ is the autoregression order, $d$ is the degree of difference, and $q$ is the moving average order. The ARIMA model addresses non-stationary time series by modeling the difference stationary time series using an ARMA model (44). It can function as an ARMA, AR, I, or MA model. The AR(p) model linearly relates the current value or observed value $y_t$ at time t of the time series to its past values $y_{t-1},y_{t-2},\dots ,y_{t-p}$ and current residuals ${ϵ}_t$, while the MA(q) model linearly relates the current value of the time series to its current ${ϵ}_t$ and past residual ${ϵ}_{t-1},{ϵ}_{t-2},\dots ,{ϵ}_{t-q}$ values. The ARIMA model is basically an ARMA model fitted on $d$-th order differenced time series such that the final differenced time series is stationary. $Y_t$ is p’th order autoregressive process, written AR(p) and is defined as

$Y_t={\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+\dots +{\varphi }_p{Y}_{t-p}+{ϵ}_t.(15)$

Furthermore (Equation 15), $Y_t$ constitutes a general linear process when expressed as $Y_t={ϵ}_t+{\psi }_1{ϵ}_{t-1}++{\psi }_2{ϵ}_{t-2}+\cdots$, with the condition that $\sum _{i=1}^{\mathrm{\infty }}{\psi }_i^2<\mathrm{\infty }$. When this equation has only a finite number of non-zero $\psi$ coefficients, it is referred to as a moving average process. Specifically, the moving average of order $q$ is (Equation 16):

$Y_t=-{\theta }_1{ϵ}_{t-1}-{\theta }_2{ϵ}_{t-2}-\dots -{\theta }_q{ϵ}_{t-q}+{ϵ}_t.(16)$

Finally, $Y_t$ is an autoregressive moving average of orders $p$ and $q$, written ARMA(p,q), if it can be written

$Y_t={\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+\dots +{\varphi }_p{Y}_{t-p}+{ϵ}_t-{\theta }_1{ϵ}_{t-1}-{\theta }_2{ϵ}_{t-2}+\dots -{\theta }_q{ϵ}_{t-q},(17)$

where the autoregressive and moving average parameters are $\varphi$ and $\theta$, respectively. Note that Equation 17 represents the ARIMA(p,q,d) model with $y_t$, where the $d$-th order difference $Y_t={\mathrm{\Delta }}^d{y}_t$ is substituted into Equation 17. We have used the auto$\mathrm{\_}$arima function to fit the ARIMA model for the univariate time series according to a provided information criterion (either AIC, BIC, or HQIC).

3.3.2 Exponential smoothing model

The Exponential Smoothing model is a simple but effective method for forecasting time series data. It depends more on the most recent observations for predicting future time series values compared to older observations. The smoothing parameter $\alpha$ controls the weight given to the most recent observation vs. the previous forecasted value, and it is typically chosen by optimization techniques to minimize the forecast error. The basic Exponential Smoothing model, as denoted by Equation 18, is presented as follows:

${\acute{y}}_{t+1|t}=\alpha y_t+(1-\alpha ){\acute{y}}_{t|t-1},\text{,for}t=1,\cdots ,T,(18)$

where, ${\acute{y}}_{t+1|t}$ is the one-step-ahead forecast for the next time period $t+1$, $y_t$ is the actual observation at time $t$, ${\acute{y}}_{t|t-1}$ is the forecast for time period $t$ based on the information available up to time $t-1$, and $\alpha$ is the smoothing parameter, which takes values between 0 and 1 and controls the weight given to the most recent observation.

3.3.3 CNN, LSTM, and bidirectional LSTM models

CNN is short for Convolutional Neural Network, which is a deep learning algorithm primarily used for image and video processing, as well as other types of data. For time series forecasting, a 1D-CNN architecture is used, consisting of an input layer, convolutional layer, pooling layer, flattened layer, fully connected layer, and output layer. The convolutional layers apply filters to the input data to identify features, while the pooling layers reduce the data’s dimensions to reduce computation. A flattened layer converts the reduced-size feature map into a one-dimensional array. Finally, the fully connected layers perform classification or regression based on the input data’s features. We have used ReLU as the activation function in the convolutional layer.

The Long Short-Term Memory (LSTM) model architecture is a type of recurrent neural network (RNN) designed to handle the vanishing gradient problem often encountered in traditional RNNs. LSTM networks have a complex structure that allows them to capture long-term dependencies in sequential data. Briefly, the basic structure of the LSTM unit has a memory cell, and three primary gates: an input gate, an output gate, and a forget gate. These gates determine whether to allow information to pass through the cell or forget it. The input gate decides which information from the current input and the previous hidden state should be passed to the current cell state, while the forget gate determines which information from the previous cell state should be retained or forgotten. Finally, the output gate determines the amount of output that should be passed to the next LSTM cell or the final output. The internal architecture of the LSTM model can be observed in Esmail et al. (54).

Unlike standard LSTM, the bidirectional LSTM (BiLSTM) is a recurrent neural network (RNN) architecture that processes sequential data in both forward and backward directions. It combines two LSTMs, one processing the input sequence in the forward direction and the other in the reverse direction. By considering each element’s past and future context in the sequence, BiLSTM can capture more comprehensive information. The main advantage of using BiLSTM is that they can capture both past and future context for each element in a sequence, allowing the model to make more informed predictions. However, BiLSTMs also have some limitations. They require processing the entire input sequence in both forward and backward directions, which can be computationally expensive and time-consuming. Overall, bidirectional LSTMs provide a powerful tool for capturing context and dependencies in sequential data by considering both past and future information.

3.3.4 Hybrid deep learning models (CNN-LSTM, CNN-bidirectional LSTM)

The hybrid deep learning model, CNN-LSTM, combines the strengths of CNNs and LSTM networks for sequence classification, including natural language processing and time series analysis. Similarly, CNN-BiLSTM is a neural network architecture that combines CNNs with BiLSTM layers. The combination of CNN and BiLSTM layers in this architecture allows the model to effectively capture both spatial and sequential dependencies in the input data. The CNN layer extracts relevant features and reduces dimensionality, while the LSTM layer captures temporal dependencies and models sequence data. The model takes input data, passes it through the CNN layer, and feeds the output into the LSTM or BiLSTM layers to produce the final output using the dense layer. Various applications such as speech recognition, sentiment analysis, and weather forecasting have successfully utilized the CNN-LSTM model. It is a powerful tool for analyzing complex sequential data requiring feature extraction and sequence modeling (Figure 9).

Figure 9

Figure 9. Neural network architecture: CNN, LSTM, BiLSTM, CNN-LSTM, and CNN-BiLSTM architectures involve neural network layers for feature extraction on input data. Model descriptions for time series prediction problems from sequences of data are provided in the Appendix (Tables A4–A9).

4 Results

We implement five deep learning models: CNN, LSTM, BiLSTM, hybrid CNN-LSTM, and CNN-BiLSTM prediction model to forecast the Mpox virus transmission. A 60-10-30 split in training, validation, and test data will be applied. We evaluate the performance of the data-driven models on test data. We also compare five deep learning models with the statistical time series ARIMA and Exponential Smoothing models. Finally, we evaluate model predictions with the real Mpox data.

4.1 Results on deterministic model

Vaccination for all individuals globally is really challenging, but vaccination is one of the effective strategies to prevent any epidemic like COVID-19. Note that the smallpox vaccine is 85% effective against Mpox due to its genetic similarity to the smallpox virus (55, 56). Since no vaccine is 100% effective (57), and people who have been Mpox vaccinated can still get Mpox, so we assume the failure rate, $\tau$ in our proposed deterministic model given in Equation (5). However, CDC (58) reported that approximately 1,233,453 people in the USA received the Mpox vaccine as of June 20, 2023, and only 23% of the population at risk has been fully vaccinated in the USA. Moreover, there are limited vaccinations in many areas of the world; for example, the Africa Centre for Disease Control stated that the continent, with a population of over 1.2 billion people, had no access to Mpox vaccines (55). Moreover, Figure 11 illustrates the significance of vaccination for the US population.

We obtained values from the literature for most biological parameters, which in this model are either constant (Table 2) or data fitting with our model (Table 3). To simulate our model, we initialize the week one infected individual by one, i.e., $I_h(0)=1,I_r(0)=1$ by ignoring the scenario of whether there were any reported cases or not. We then fit our Mpox model with the reported data and uncover the fitted parameters in Table 3 and the variable parameters in Table 4.

4.1.1 Fitting deterministic model

Suppose we are given the data $\{(t_1,y_1),\dots (t_n,y_n)\}$, where $y_n$ is the $n$-th observed data point, and $n$ is the total number of data points. We want to fit $\eta I(t)$ of model estimated parameter $\vartheta$, where $\eta$ is a “reporting rate.” Then the sum of squared errors (SSE) between $\eta I(t)$ and the data can be measured by

$L:=\text{SSE}(\vartheta )=\sum _{j=1}^n\Vert y_j-\eta I(t_j){\Vert }^2.(19)$

In the least-square fitting, we find the value $\hat{\vartheta }$ of the model parameter $\vartheta$ such that $\text{SSE}(\vartheta )$ is the minimum. We use the Matlab function ode45 to simulate our proposed model, and Matlab’s minimization-constrained function fmincon takes the least-squares error function $\text{SSE}(\vartheta )$ and uses a direct search routine to find a minimum value of least squares error. In this fitting, we fix most of the parameters given in Table 2 and estimate the parameters $\vartheta =[c_0,k_1,k_2,{\alpha }_h,{\alpha }_r,\tau ]$ using the reporting rate $\eta ={10}^{-1}$ on the model incidence of $I$. More importantly, Figure 10 presents the Mpox deterministic model fitting in various locations, and its algorithmic steps are given in Appendix A (Algorithm 1).

To investigate the particular impact of disease transmission of Mpox via inter-group contacts, we now generate colormaps based on the baseline contact rate $c_0$ within the overall population and the proportion of human contacts, $k_1$. To incorporate these numerical simulations, an outbreak is defined as in (61, 62) when the maximum value of the function I(t) within a specified time interval $\tau$ is greater than or equal to 2, represented as $\underset{t\in \tau }{max}I(t)\ge 2$. On the other hand, the absence of an outbreak is determined when the maximum value of infected cases, I(t), within the same time frame, remains below 2, i.e., $\underset{t\in \tau }{max}I(t)<2$, where minor fluctuations in the total count of infected individuals may arise due to numerical instability. These conditions produce the epidemic (yellow-colored) vs. non-epidemic (blue-colored) parametric regions in Figure 12.

4.2 Results on data-driven models

We now implement the deep learning model to the EPI weekly time series data and found that the hybrid CNN-LSTM network model performed best in both test and prediction. We choose different window sizes to group the confirmed cases in the time series data for training and testing the five neural networks: CNN, LSTM, BiLSTN, and hybrid CNN-LSTM and CNN-BiLSTM. Then, the optimal window size and hyperparameters were determined by analyzing their errors’ (MAE, RMSE, nRMSE) effect on the network’s performance. Several studies, such as the one by Sun et al. (63) in 2020, have demonstrated the practical effectiveness of the Adam optimizer and its favorable performance compared to other adaptive learning rate algorithms. As a result, we have implemented the Adam optimizer in all of our machine-learning models.

5 Discussion

While recent deep learning studies on the Mpox approach stayed limited to classifying or diagnosing Mpox, for example, (64) and many more, our study modeled to explore the dynamics of disease transmission. In our modeling framework, we have done model fitting in Figure 10 and made predictions in Figure 15 and found important epidemiological parameters in Table 3. More importantly, our findings of the deterministic model (Equation 5) demonstrate that with an increasing vaccination rate, the percentage of the epidemic size decreases, as depicted in Figure 11. This result also highlights the importance of vaccination for disease control and prevention. The findings presented in Figure 12 indicate that effective control of the epidemic can be achieved by making strategic modifications to the baseline contact rate $c_0$ and the proportion of contacts within the human population $k_1$, as depicted in the accompanying Figure. Additionally, we have constructed the deep learning models and presented a comprehensive visualization of the entire time series dynamics in Figure 13. Next, Figure 14 depicts the disease dynamics observed on the test set for the data-driven models across various regions: (A) World, (B) North America, (C) South America, (D) Europe, (E) Brazil, and (F) USA. This figure also reveals that the deep learning models, namely CNN, LSTM, BiLSTM, CNN-LSTM, and CNN-BiLSTM, closely align with the reported incidence dynamics, capturing the trends effectively. However, the ARIMA and Exponential smoothing models deviate from the real incidence on test data, particularly on the World data.

Figure 10

Figure 10. Data fitting for the Mpox deterministic model in various locations: The presented figure illustrates the Mpox data fitting of the deterministic model in various locations. As a mathematical model of the Mpox spread, we use the deterministic model shown in Figure A2. Also, Table 2 provides the initial data used for the fitting, while Table 3 contains the corresponding fitting parameters for each location. Additionally, Table 4 shows various variable parameters using the testing errors defined in Equation 2.

Figure 11

Figure 11. Importance of Vaccinations: This figure illustrates the significant impact of vaccination on the USA population size $N_h=331\text{,}002\text{,}651$. Panel (A) depicts the percentage of epidemic sizes for ${\varphi }_0\in [0,0.25]$. Panel (B) demonstrates a decreasing fraction of incidence with increasing vaccination rates. Panels (C) and (D) depict the cumulative incidence and the prevalence, respectively, highlighting a substantial reduction when the vaccination rate ${\varphi }_0$ increases. Overall, raising the vaccination rate $\varphi$ of Mpox offers a solution to mitigate the risks posed by Mpox as well as other emerging infectious diseases. For these simulations, we use the parameter values given in Table 2 except the variable parameter values ${\beta }_{hh}=5.5,{\beta }_{rr}=2.5,{\beta }_{hr}={\beta }_{rh}=2.0251$. For each simulation, we also consider the initial conditions $I_h(0)=I_r(0)=1,S_h(0)=N_h-I(0)-V_1(0),E_h(0)=H(0)=V_1(0)=V_2(0)=R_h(0)=E_r(0)=R_r=0,S_r(0)=N_r-I_r(0)$, where $N_r=8\text{,}000\text{,}000$. We run simulations over the time range $\tau =[0,65]$ with step size $\mathrm{\Delta }\tau =0.1$.

Figure 12

Figure 12. Epidemic (yellow colored) vs. non-epidemic (blue colored) parametric regions depending on $c_0$ and $k_1$: The figure shows the colormaps for the range of $c_0=[0,1]$, and $k_1=[0,1]$ with step size $\mathrm{\Delta }{c}_0=\mathrm{\Delta }{k}_1=0.1$. The epidemic region is depicted in yellow in the figure, while the non-epidemic parametric region is represented in blue. For these simulations, we use the parameter values given in Table 2, fitted parameter in Table 3, and variable parameters value in Table 4. This figure also illustrates that with these epidemiological parameters, it is possible to control the epidemic by adjusting the baseline contact rate $c_0$ and the proportion of contacts within the human population $k_1$ in overall contacts. For each simulation, we also consider the initial conditions $I_h(0)=I_r(0)=1,V_1(0)=N_h\ast (1/100),S_h(0)=N_h-I(0)-V_1(0),E_h(0)=H(0)=V_2(0)=R_h(0)=E_r(0)=R_r=0,S_r(0)=N_r-I_r(0)$, where $N_r=8\text{,}000\text{,}000$. We run simulations over the time domain $\tau =[0,65]$ with step size $\mathrm{\Delta }\tau =0.01$.

Figure 13

Figure 13. Deep learning model training for Mpox incidence in different locations: This figure shows the predicted vs. actual incidence per week in six different locations worldwide using five different deep learning models. We use 60% for the training dataset with 10% for validation during the model building (details given in the appendix) and then utilize the trained model to predict the incidence (number of new cases) per week for the entire time series data spanning from 2023-02-11 to 2023-06-03. For example, in the World and Europe time series data, the model was trained on the data from 2022-05-07 to 2023-02-04 (week 0 to week 39), and the predictions were then extended to cover the entire time whole period from 2022-05-07 to 2023-06-03 (week 0 to week 56) and compared with the actual data. The actual data are displayed in red in the panel, while the predictions from various deep learning models, including CNN, LSTM, BiLSTM, CNN-LSTM, and CNN-BiLSTM, are presented in magenta, green, blue, black, and cyan, respectively. These models were also evaluated using test data from 2023-02-11 to 2023-06-03 (week 40 to week 56), with the model predictions given in Figure 14 and subsequently used to predict the next eight weeks from 2023-06-10 to 2023-07-29 (week 57 to week 64) in Figure 15. Similarly, we optimized the model for other locations as provided in the Appendix C.

Figure 14

Figure 14. Data driven model on test data for Mpox incidence across different locations.

Figure 15

Figure 15. Mpox predictions of different locations using data-driven and deterministic models. Appendix C contains the data-driven model construction process, along with error measurements on the test data. Additional details, i.e., model parameters and model construction are provided in Tables A1, A2, A3, A4, A5, A6, A7, A8, and A9. Moreover, the deterministic model used the same parameter values described in Figure 10.

In the context of the prediction shown in Figure 15, the hybrid deep learning models, CNN-LSTM and CNN-BiLSTM, exhibit good performance compared to other models across all test dynamics, demonstrating their potential for accurate disease prediction and capturing the complexities of the Mpox outbreak in different regions. Notably, most of the predictive models indicate that the Mpox disease is currently in a decline phase. The predictive trend remarkably mirrors the real data. However, it is worth noting that some models, for example, CNN, LSTM, BiLSTM, and CNN-BiLSTM indicate the possibility of another peak in the disease’s trajectory in selected geographical regions. In terms of overall performance, it is evident that both CNN-LSTM and deterministic models exhibit better predictive capabilities. Their accuracy (Table A10) and effectiveness in forecasting Mpox disease dynamics appear to stand out among the various models considered in this analysis.

5.1 Modeling limitation

The deterministic model is designed to incorporate explainability and account for real-world scenarios, although it does rely on specific assumptions. In contrast, data-driven modeling offers less explainability in their predictions but demonstrates good predictive capability when applied to time series data. We normalized the complete time series dataset using a data normalization approach, which involved selecting max and min values from each location’s entire time series. Note that this normalization data is only used in deep learning modeling. Another limitation of the deep learning models is that they are trained on random initial weights during the training.

6 Conclusion

This paper conducted a comprehensive analysis of global Mpox univariate time series data in diverse geographical locations. We proposed a deterministic model and utilized advanced deep learning techniques such as 1D-CNN, LSTM, BiLSTM, hybrid CNN-LSTM, and CNN-BiLSTM, alongside statistical time series models like ARIMA and exponential smoothing, to gain deeper insights into the Mpox disease dynamics. Moreover, the spatial pattern analysis of global Mpox data from May 1, 2022, to May 31, 2023, offered insights into the geographical distribution of the disease, helping public health authorities and policymakers focus on areas with higher risks. Our deterministic model highlighted the critical role of vaccination rates in flattening the curve of infection dynamics and influencing the basic reproduction number. It underscored the importance of increasing vaccination among susceptible populations to control disease transmission effectively. Through data fitting, we estimated crucial epidemiological parameters within our proposed deterministic model. The results showed the importance of reducing contact rates in high-risk groups to mitigate the disease outbreaks. Additionally, these findings contributed to a comprehensive understanding of disease dynamics across diverse locations and informed targeted intervention strategies for controlling and mitigating infectious disease outbreaks. Furthermore, our deterministic and data-driven models extended their utility by providing short-term (eight weeks) predictions across various geographical locations, including the World, USA, Brazil, and three continents: North America, South America, and Europe. The findings, verified by real data, suggested that Mpox is approaching its decline phase as of July 29, 2023.

7 Computer configuration

We performed all experiments in this study using Python programming language on Google Colab with the TensorFlow and Keras libraries through the Google Chrome browser and MATLAB. Computer info: MacBook Pro; Chip: AppleM2 Pro; Memory: 16 GB.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

HKD: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Software, Supervision, Validation, Visualization, Writing – original draft, Project administration, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article.

This research was partially supported by Graduate College (OSU) and GPSGA research material grant (OSU).

Acknowledgments

The author is grateful to Dr. Md Sakhawat Hossain (Clemson University) for carefully reading the manuscript. The author would like to acknowledge Dr. Lucas M. Stolerman (OSU), Dr. Chanaka Kottegoda (OSU), S M Rakibur Rahman (OSU), Dr. Pedro Maia (UTA), Dr. Sathyanarayanan N. Aakur (Auburn University), and Sanjoy Kundu (Auburn University) for valuable conversations. Thanks to the organizers of the 15th SISMID, Dept. of Biostatistics, University of Washington, for the financial support, which facilitated much of the work conducted during the training period. The author also thanks Dr. Micaela Martinez (WE ACT for Environmental Justice) and Dr. Matt Ferrari (Pennsylvania State University) for meaningful interactions at the 2023 SISMID.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Piret J, Boivin G. Pandemics throughout history. Front Microbiol. (2021) 11:1–16. doi: 10.3389/fmicb.2020.631736

Crossref Full Text | Google Scholar

2. Bunge EM, Hoet B, Chen L, Lienert F, Weidenthaler H, Baer LR, et al. The changing epidemiology of human monkeypox—a potential threat? a systematic review. PLoS Negl Trop Dis. (2022) 16:e0010141. doi: 10.1371/journal.pntd.0010141

PubMed Abstract | Crossref Full Text | Google Scholar

3. Hsu J. Here’s how computer models simulate the future spread of new coronavirus. Sci Am. (2020). Available online at: https://www.scientificamerican.com/article/heres-how-computer-models-simulate-the-futurespread-of-new-coronavirus/.34276078

PubMed Abstract | Google Scholar

4. Kuhl E. Data-driven modeling of COVID-19-lessons learned. Extreme Mech Lett. (2020) 40:100921. doi: 10.1016/j.eml.2020.100921

PubMed Abstract | Crossref Full Text | Google Scholar

5. Wolfe ND, Dunavan CP, Diamond J. Origins of major human infectious diseases. Nature. (2007) 447:279–83. doi: 10.1038/nature05775

PubMed Abstract | Crossref Full Text | Google Scholar

6. Olaniyi S, Chuma F. Lyapunov stability and economic analysis of monkeypox dynamics with vertical transmission and vaccination. Int J Appl Comput Math. (2023) 9:679–91. doi: 10.1007/s40819-023-01572-w

Crossref Full Text | Google Scholar

7. Caminade C, McIntyre KM, Jones AE. Impact of recent and future climate change on vector-borne diseases. Ann N Y Acad Sci. (2019) 1436:157–73. doi: 10.1111/nyas.13950

PubMed Abstract | Crossref Full Text | Google Scholar

8. Madubueze CE, Onwubuya IO, Nkem GN, Chazuka Z. The transmission dynamics of the monkeypox virus in the presence of environmental transmission. Front Appl Math Stat. (2022) 8:1–21. doi: 10.3389/fams.2022.1061546

Crossref Full Text | Google Scholar

9. Giulio DBD, Eckburg PB. Human monkeypox: an emerging zoonosis. Lancet Infect Dis. (2004) 4:15–25. doi: 10.1016/s1473-3099(03)00856-9

PubMed Abstract | Crossref Full Text | Google Scholar

10. Chen Y, Li M, Fan H. The monkeypox outbreak in 2022: adaptive evolution associated with APOBEC3 may account for. Signal Transduct Target Ther. (2022) 7:323. doi: 10.1038/s41392-022-01181-x

PubMed Abstract | Crossref Full Text | Google Scholar

11. McCollum AM, Damon IK. Human monkeypox. Clin Infect Dis. (2014) 58:260–7. doi: 10.1093/cid/cit703

PubMed Abstract | Crossref Full Text | Google Scholar

12. Ladnyj ID, Ziegler P, Kima E. A human infection caused by monkeypox virus in Basankusu Territory. Bull World Health Organ. (1972) 46:593–7. Available online at: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2480792/.4340218

PubMed Abstract | Google Scholar

13. Centers for Disease Control and Prevention (CDC). Monkeypox (2021). Available online at: https://www.cdc.gov/monkeypox/index.html (Accessed July 29, 2023).

14. Parker S, Buller RM. A review of experimental and natural infections of animals with monkeypox virus between 1958 and 2012. Future Virol. (2013) 8:129–57. doi: 10.2217/fvl.12.130

PubMed Abstract | Crossref Full Text | Google Scholar

15. Reynolds MG, Yorita KL, Kuehnert MJ, Davidson WB, Huhn GD, Holman RC, et al. Clinical manifestations of human monkeypox influenced by route of infection. J Infect Dis. (2006) 194:773–80. doi: 10.1086/505880

PubMed Abstract | Crossref Full Text | Google Scholar

16. Reynolds MG, Davidson WB, Curns AT, Conover CS, Huhn G, Davis JP, et al. Spectrum of infection and risk factors for human monkeypox. Emerging Infect Dis. (2007) 13:1332–9. doi: 10.3201/eid1309.070175

Crossref Full Text | Google Scholar

17. Siegrist EA, Sassine J. Antivirals with activity against mpox: a clinically oriented review. Clin Infect Dis. (2023) 76:155–64. doi: 10.1093/cid/ciac622

PubMed Abstract | Crossref Full Text | Google Scholar

18. Khan MI, Qureshi H, Bae SJ, Awan UA, Saadia Z, Khattakb AA. Predicting monkeypox incidence: fear is not over! Journal of Infection. (2023) 86:256–308. doi: 10.1016/j.jinf.2022.12.021

PubMed Abstract | Crossref Full Text | Google Scholar

19. Thornhill JP, Barkati S, Walmsley S, Rockstroh J, Antinori A, Harrison LB, et al. Monkeypox virus infection in humans across 16 countries—April–June 2022. N Engl J Med. (2022) 387:679–91. doi: 10.1056/NEJMoa2207323

PubMed Abstract | Crossref Full Text | Google Scholar

20. Kanj El Eid R, Allaw F, Haddad SF, Kanj SS. Human monkeypox: a review of the literature. PLoS Pathog. (2022) 18:e1011008. doi: 10.1371/journal.ppat.1010768

PubMed Abstract | Crossref Full Text | Google Scholar

21. Hasan S, Saeed S. Monkeypox disease: an emerging public health concern in the shadow of COVID-19 pandemic—an update. Trop Med Infect Dis. (2022) 7:283. doi: 10.3390/tropicalmed7100283

PubMed Abstract | Crossref Full Text | Google Scholar

22. Jeta M, Sekkak I, Ortiz AM, Moyles I, Nasri B. Mathematical modeling of monkeypox: a scoping review. One Health. (2023) 16:1–8. doi: 10.1016/j.onehlt.2023.100540

Crossref Full Text | Google Scholar

23. Tang S, Gao S, Zhang F, Liu Y. Role of vector resistance and grafting infection in huanglongbing control models. Infect Dis Model. (2023) 8(2):491–513. doi: 10.1016/j.idm.2023.04.006

PubMed Abstract | Crossref Full Text | Google Scholar

24. Yuan P, Tan Y, Yang L, Aruffo E, Ogden NH, Bélair J, et al. Assessing transmission risks and control strategy for monkeypox as an emerging zoonosis in a metropolitan area. J Med Virol. (2023) 95(1):1–15. doi: 10.1002/jmv.28137

Crossref Full Text | Google Scholar

25. Esteban EP, Almodovar-Abreu L. Assessing the impact of vaccination in a COVID-19 compartmental model. Inform Med Unlocked. (2021) 27:1–10. doi: 10.1016/j.imu.2021.100795

Crossref Full Text | Google Scholar

26. Dong S, Xu L. Application of a time-delay sir model with vaccination in COVID-19 prediction and its optimal control strategy. Res Sq. (2022) 111:10677–92. doi: 10.1007/s11071-023-08308-x

PubMed Abstract | Crossref Full Text | Google Scholar

27. El-Mesady A, Elsonbaty A, Adel W. On nonlinear dynamics of a fractional order monkeypox virus model. Chaos Solitons Fract. (2022) 164:112716. doi: 10.1016/j.chaos.2022.112716

Crossref Full Text | Google Scholar

28. Sitaula C, Shahi TB. Monkeypox virus detection using pre-trained deep learning-based approaches. J Med Syst. (2022) 46:78. doi: 10.1007/s10916-022-01868-2

PubMed Abstract | Crossref Full Text | Google Scholar

29. Chimmula V, Zhang L. Time series forecasting of COVID-19 transmission in Canada using LSTM networks. Chaos Solitons Fract. (2020) 135:109864. doi: 10.1016/j.chaos.2020.109864

Crossref Full Text | Google Scholar

30. Oshinubi K, Amakor A, Peter OJ, Rachdi M, Demongeot J. Approach to COVID-19 time series data using deep learning and spectral analysis methods. AIMS Bioeng. (2022) 9:1. doi: 10.3934/bioeng.2022001

Crossref Full Text | Google Scholar

31. Markova M. Convolutional neural networks for forex time series forecasting. AIP Conf Proc. (2022) 2459:030024. doi: 10.1063/5.0083533

Crossref Full Text | Google Scholar

32. McAndrew T, Majumder MS, Lover AA, Venkatramanan S, Bocchini P, Besiroglu T, et al. Early human judgment forecasts of human monkeypox. Lancet Digit Health. (2022) 4:e569–71. doi: 10.1016/S2589-7500(22)00127-3

PubMed Abstract | Crossref Full Text | Google Scholar

33. https://ourworldindata.org/monkeypox (Accessed June 03, 2023 and July 29, 2023).

34. Stolerman LM, Maia PD, Kutz JN. Forecasting dengue fever in Brazil: an assessment of climate conditions. PLoS One. (2019) 14(8):1–16. doi: 10.1371/journal.pone.0220106

PubMed Abstract | Crossref Full Text | Google Scholar

35. https://www.worldometers.info/ (Accessed July 29, 2023).

36. Tealab A, Hefny H, Badr A. Forecasting of nonlinear time series using ANN. Future Comput Inform J. (2017) 2:39–47. doi: 10.1016/j.fcij.2017.05.001

Crossref Full Text | Google Scholar

37. McFarland JM, Ho ZV, Kugener G, Dempster JM, Montgomery PG, Bryan JG, et al. Improved estimation of cancer dependencies from large-scale RNAi screens using model-based normalization and data integration. Nat Commun. (2018) 9:1–13. doi: 10.1038/s41467-018-06916-5

PubMed Abstract | Crossref Full Text | Google Scholar

38. Wibawa AP, Utama ABP, Elmunsyah H, Pujianto U, Dwiyanto FA, Hernandez L. Time-series analysis with smoothed convolutional neural network. J Big Data. (2022) 9:1–18. doi: 10.1186/s40537-022-00599-y

Crossref Full Text | Google Scholar

39. Kuhn M, Johnson K. Applied Predictive Modeling. Vol. 9. 2nd ed. Groton, CT: Pfizer Global Research and Development (2018). p. 44.

Google Scholar

40. Nguyen QH, Ly HB, Ho LS, Al-Ansari N, Le HV, Tran VQ, et al. Influence of data splitting on performance of machine learning models in prediction of shear strength of soil. Math Probl Eng. (2021) 2021:1–15. doi: 10.1155/2021/4832864

Crossref Full Text | Google Scholar

41. Crick C. Machine Learning Class on Loss Functions. Oklahoma State University (OSU) (2023).

Google Scholar

42. Olaniyi S, Obabiyi OS, Okosun KO, Oladipo AT, Adewale SO. Mathematical modelling and optimal cost-effective control of COVID-19 transmission dynamics. Eur Phys J Plus. (2020) 28(11):1–20. doi: 10.1140/epjp/s13360-020-00954-z

Crossref Full Text | Google Scholar

43. Bikia V, Rovas G, Stergiopulos N. Cardiac output estimated from an uncalibrated radial blood pressure waveform: validation in an in-silico- generated population. Front Bioeng Biotechnol. (2023) 11:1199726. doi: 10.3389/fbioe.2023.1199726

PubMed Abstract | Crossref Full Text | Google Scholar

44. Ceylan Z. Estimation of COVID-19 prevalence in Italy, Spain, and France. Sci Total Environ. (2020) 729:138817. doi: 10.1016/j.scitotenv.2020.138817

PubMed Abstract | Crossref Full Text | Google Scholar

45. https://vaccinateyourfamily.org/ (Accessed June 20, 2023).

46. https://www.fda.gov/media/131078/download (Accessed June 20, 2023).

47. (2022). https://www.fda.gov/vaccines-blood-biologics/vaccines/acam2000 (Accessed June 20, 2023).

48. Peter OJ, Kumar S, Kumari N, Oguntolu FA, Oshinubi K, Musa R. Transmission dynamics of monkeypox virus: a mathematical modelling approach. Model Earth Syst Environ. (2022) 8:3423–34. doi: 10.1007/s40808-021-01313-2

PubMed Abstract | Crossref Full Text | Google Scholar

49. Peter OJ, Oguntolu FA, Ojo MM, Oyeniyi AO, Jan R, Khan I. Fractional order mathematical model of monkeypox transmission dynamics. Phys Scr. (2022) 97:1–25. doi: 10.1088/1402-4896/ac7ebc

Crossref Full Text | Google Scholar

50. Diekmann O, Heesterbeek JA. Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis, and Interpretation. Wiley Series in Mathematical and Computational Biology. (2000).

Google Scholar

51. Diekmann O, Heesterbeek JAP, Metz JAJ. On the definition and computation of the basic reproduction ratio ${\mathcal{R}}_0$ in models for infectious diseases in heterogeneous populations. J Math Biol. (1990) 28:365–82. doi: 10.1007/BF00178324

PubMed Abstract | Crossref Full Text | Google Scholar

52. Ngungu M, Addai E, Adeniji A, Adam UM, Oshinubi K. Mathematical epidemiological modeling and analysis of monkeypox dynamism with non-pharmaceutical intervention using real data from United Kingdom. Front Public Health. (2023) 11:1–16. doi: 10.3389/fpubh.2023.1101436

Crossref Full Text | Google Scholar

53. Akinyemi ST, Idisi IO, Rabiu M, Okeowo VI, Iheonu N, Dansu EJ, et al. A tale of two countries: optimal control and cost-effectiveness analysis of monkeypox disease in Germany and Nigeria. Healthc Anal. (2023) 4:1–26. doi: 10.1016/j.health.2023.100258

Crossref Full Text | Google Scholar

54. Esmail W, Stockmanns T, Ritman J. Machine learning for track finding at panda. In: Proceedings of the CTD/WIT (2019). https://arxiv.org/abs/1910.07191

55. Ogunkola IO, Abiodun OE, Bale BI, Elebesunu EE, Ujam SB, Umeh IC, et al. Monkeypox vaccination in the global south: fighting a war without a weapon. Clin Epidemiol Glob Health. (2023) 22:1–4. doi: 10.1016/j.cegh.2023.101313

Crossref Full Text | Google Scholar

56. Turner M, Mandia J, Keltner C, Haynes R, Faestel P, Mease L. Monkeypox in patient immunized with ACAM2000 smallpox vaccine during 2022 outbreak. Emerging Infect Dis. (2022) 18:2336–8. doi: 10.3201/eid2811.221215

Crossref Full Text | Google Scholar

57. Farrar JL, Lewis NM, Houck K, Canning M, Fothergill A, Payne AB. Demographic and clinical characteristics of mpox in persons who had previously received 1 dose of JYNNEOS vaccine and in unvaccinated persons - 29 U.S. jurisdictions, May 22–September 3, 2022. Am J Transplant. (2022) 23:298–303. doi: 10.1016/j.ajt.2023.01.003

Crossref Full Text | Google Scholar

58. Vaccines data for 2023 response. Available online at: https://www.cdc.gov/poxvirus/mpox/response/2022/vaccines_data.html (Accessed June 20, 2023).

59. Miura F, Backer JA, Rijckevorsel GV, Bavalia R, Raven S, Petrignani M. Time scales of human mpox transmission in the Netherlands. J Infect Dis. (2023) 229(3):800–4. doi: 10.1093/infdis/jiad091

PubMed Abstract | Crossref Full Text | Google Scholar

60. DeWitt ME, Polk C, Williamson J, Shetty AK, Passaretti CL, McNeil CJ. Global monkeypox case hospitalisation rates: a rapid systematic review and meta-analysis. eClinicalMedicine. (2022) 54:1–16. doi: 10.1016/j.eclinm.2022.101710

Crossref Full Text | Google Scholar

61. Stolerman LM, Coombs D, Boatto S. Sir-network model and its application to dengue fever. SIAM J Appl Math. (2015) 75:2581–609. doi: 10.1137/140996148

Crossref Full Text | Google Scholar

62. Das HK, Stolerman LM. Epidemic thresholds and disease dynamics in metapopulations: the role of network geometry and human mobility. arXiv. arXiv:2308.01998 X, 1–40 (2023).

63. Sun S, Cao Z, Zhu H, Zhao J. A survey of optimization methods from a machine learning perspective. IEEE Trans Cybern. (2020) 50:3668–81. doi: 10.1109/TCYB.2019.2950779

PubMed Abstract | Crossref Full Text | Google Scholar

64. Eliwa EHI., El Koshiry AME, El-Hafeez TA, Farghaly HM. Utilizing convolutional neural networks to classify monkeypox skin lesions. Sci Rep. (2023) 13:1–20. doi: 10.1038/s41598-023-41545-z

Crossref Full Text | Google Scholar

Appendix

The purpose of this appendix is to present a few results that are outlined in the main manuscript.

A: Mpox daily incidence and the schematic representation of the deterministic model

Figure A1

Figure A1. Dynamics of daily Mpox incidence: This figure illustrates the daily incidence of Mpox in various locations from May 2022 to May 2023. The figure exhibits seasonality and includes missing reported data. In contrast, the weekly data in Figure 2C captures the proper disease dynamics, which can be a pragmatic choice.

Figure A2

Figure A2. Schematic representation of the deterministic model: This figure shows a modified SEIR model with two-phase vaccinations (smallpox or dose one vaccine and dose two vaccines) for mpox virus transmissions.

B. Algorithm for data fitting via least squares estimation in the deterministic model

Algorithm 1

Algorithm 1. Parameter Estimation for a Deterministic Model

C. Parameters and hyperparameters for data-driven modeling

For data-driven modeling, we use different libraries in Python. First, we utilize auto_arima to train the ARIMA model from the pmdarima library. Next, we fit the Exponential Smoothing model to the data using the ExponentialSmoothing function. Then, we employ the Keras deep learning library to develop the neural network models.

Table A1

Table A1. Optimized ARIMA model selection based on location and error on test data.

Table A2

Table A2. Optimized Exp Smoothing model selection based on location and performance analysis.

Table A3

Table A3. Optimized Exp Smoothing model’s initial seasonal coefficients.

Table A4

Table A4. Hyperparameter tuning on World data using deep learning approach.

Table A5

Table A5. Hyperparameter tuning on North America data using deep learning approach.

Table A6

Table A6. Hyperparameter tuning on South America data using deep learning approach.

Table A7

Table A7. Hyperparameter tuning on Europe data using deep learning approach.

Table A8

Table A8. Hyperparameter tuning on Brazil data using deep learning approach.

Table A9

Table A9. Hyperparameter tuning on USA data using deep learning approach.

Table A10

Table A10. Prediction errors for different models.