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.


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)(15)(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)(24)(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 shortterm 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 deeplearning 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.

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.

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.

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 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.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.

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 where y norm t 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.

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.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).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.

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.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.Frontiers in Epidemiology 07 frontiersin.org

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, y as the forecast value, and n as the number of observations.Then the mathematical formulas for these error evaluation measures are as follows: , and nRMSE ¼ RMSE (y max À y min ) : One can also compute nRMSE based on the average ( y) of the actual data values, calculated as RMSE y .

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 Q in susceptible class and death rate m in all classes.Also, u ¼ QN(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), 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): where the group-specific force of infections, l ¼ l h l r , the transmission matrix, T ¼ b hh b hr b rh b rr , and the infectious class, The diagonal and off-diagonal of T represent within-group and between-group interactions, respectively, and we assume that det(T) !0, as it is reasonable to disease transmission.Additionally, the first component h of b 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 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 f 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, t 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).
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: with the initial conditions in a biologically feasible region G h Â G r , where G h and G r are defined by the Equation 6 below:

Model analysis
By summing the first seven equations of group A in the (Equation 5), we have the differential equation for the human population N h is given as follows: Moreover, Figure 12 illustrates the significance of vaccination for the US population.Similarly, for the non-human population, N r , we have , and using u r ¼ Q r N r the corresponding differential is given as follows: 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

Equilibrium state
There are two types of equilibrium points for the model given in (Equation 5).One is the Disease-Free Equilibrium , u r m r , 0, 0, 0), which corresponds to the absence of infectious individuals.The second equilibrium point is called endemic equilibrium, represented by , 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 s in the system (Equation 5), we obtain the endemic equilibrium point as follows:

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 e 0 .Using the technique (50,51), the next-generation matrix (NGM), k ¼ FV À1 at the DFE is defined as where the notations , and 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 9is det(kÀlI) ¼ 0, and this gives the eigenvalues (Equation 10) where, B 1 : , and 3 is well defined and positive, that implies l 3 !l 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: 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 By substituting the values of T i for i ¼ 1, . . ., 6, u h ¼ Q h N h , and u r ¼ Q r N r into (Equation 12), we derive the resulting Remark 2. Increasing the vaccination rate f 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 R 0 , which plays a pivotal role in effectively controlling the transmission of the disease.

Sensitivity analysis of basic reproduction 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 R 0 used in Ngungu et al. (52) and Samuel et al. (53) with respect to parameter p is given by When S R0 p is positive for the parameter p, it indicates an increase in R 0 , while a negative value of S R0 p for the parameter p suggests a decrease in R 0 .Due to complexity of the actual R 0 (Equation 11), we consider the normalized sensitivity index (Equation 14) of R 0 (Equation 13) with respect to the parameter p 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.

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.

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 dth difference is Y t ¼ D d y t .ARIMA is represented as 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 Mpox model sensitivities to its associated parameters for six different locations: Panel (A) and (B) illustrate the sensitivity index of R 0 (Equation 13) for t ¼ 0 and t ¼ 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 b hh , u h , g 1 , m h , t, and f changes when the parameter values change, while the remaining parameters remain unchanged.In other words, the parameters a h , g h , d h , s, 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 R 0 (Equation 13) for six different locations using six different colors and sizes of circles.the time series to its past values y tÀ1 , y tÀ2 , . . ., y tÀp and current residuals e t , while the MA(q) model linearly relates the current value of the time series to its current e t and past residual e tÀ1 , e tÀ2 , . . ., e 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 Furthermore (Equation 15), Y t constitutes a general linear process when expressed as When this equation has only a finite number of non-zero c coefficients, it is referred to as a moving average process.Specifically, the moving average of order q is (Equation 16): Finally, Y t is an autoregressive moving average of orders p and q, written ARMA(p,q), if it can be written where the autoregressive and moving average parameters are f and u, respectively.Note that Equation 17 represents the ARIMA(p,q,d) model with y t , where the d-th order difference Y t ¼ D d y t is substituted into Equation 17.We have used the auto arima function to fit the ARIMA model for the univariate time series according to a provided information criterion (either AIC, BIC, or HQIC).

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 a 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: where, y tþ1jt is the one-step-ahead forecast for the next time period t þ 1, y t is the actual observation at time t, y tjtÀ1 is the forecast for time period t based on the information available up to time t À 1, and a is the smoothing parameter, which takes values between 0 and 1 and controls the weight given to the most recent observation.

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.

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).

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.

Results 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, t 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;  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. 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.

Fitting deterministic model
Suppose we are given the data {(t 1 , y 1 ), . . .(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 hI(t) of model estimated parameter q, where h is a "reporting rate."Then the sum of squared errors (SSE) between hI(t) and the data can be measured by In the least-square fitting, we find the value q of the model parameter q such that SSE(q) is the minimum.We use the Matlab function ode45 to simulate our proposed model, and Matlab's minimization-constrained function fmincon takes the leastsquares error function SSE(q) 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 q ¼ [c 0 , k 1 , k 2 , a h , a r , t] using the reporting rate h ¼ 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 t is greater than or equal to 2, represented as max t[t I(t) ! 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., max t[t 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.

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.

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 where N r ¼ 8,000,000.We run simulations over the time range t ¼ [0, 65] with step size Dt ¼ 0:1.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.
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.

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.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 , where N r ¼ 8,000,000.We run simulations over the time domain t ¼ [0, 65] with step size Dt ¼ 0:01.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.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.

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.

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 7 Building
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, y) is computed by comparing predicted values ( y) and true values (y), and model parameters are adjusted to minimize L(y, y) that enhances the accuracy of prediction model f.

FIGURE 9
FIGURE 9Neural 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 (TablesA4-A9).

FIGURE 10 Data
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, Table2provides the initial data used for the fitting, while Table3contains the corresponding fitting parameters for each location.Additionally, Table4shows various variable parameters using the testing errors defined in Equation2.

FIGURE 11 Importance
FIGURE 11 Importance of Vaccinations: This figure illustrates the significant impact of vaccination on the USA population size N h ¼ 331,002,651.Panel (A) depicts the percentage of epidemic sizes for f 0 [ [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 f 0 increases.Overall, raising the vaccination rate f 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 b hh ¼ 5:5, b rr ¼ 2:5, b hr ¼ b rh ¼ 2:0251.For each simulation, we also consider the initial conditionsI 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,000,000.We run simulations over the time range t ¼ [0, 65] with step size Dt ¼ 0:1.

FIGURE 12 Epidemic
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 Dc 0 ¼ Dk 1 ¼ 0:1.The epidemic region is depicted in yellow in the figure, while the nonepidemic parametric region is represented in blue.For these simulations, we use the parameter values given in Table2, fitted parameter in Table3, and variable parameters value in Table4.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 Ã(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,000,000.We run simulations over the time domain t ¼ [0, 65] with step size Dt ¼ 0:01.

FIGURE 13
FIGURE 13Deep 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 Figure14and subsequently used to predict the next eight weeks from 2023-06-10 to 2023-07-29 (week 57 to week 64) in Figure15.Similarly, we optimized the model for other locations as provided in the Appendix C.

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

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

TABLE 4
Variable parameters and errors in deterministic model fitting.

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

TABLE A3 Optimized
Exp Smoothing model's initial seasonal coefficients.
TABLE A5 Hyperparameter tuning on North America data using deep learning approach.
TABLE A6 Hyperparameter tuning on South America data using deep learning approach.
TABLE A7 Hyperparameter tuning on Europe data using deep learning approach.
TABLE A8 Hyperparameter tuning on Brazil data using deep learning approach.
TABLE A9 Hyperparameter tuning on USA data using deep learning approach.Model building on the training data (June 4, 2022, to February 11, 2023) with 10% validation of training data Test : Feb. 18-May 27 (2023) TABLE A10 Prediction errors for different models.