# Hand, Foot, and Mouth Disease in Thailand: A Comprehensive Modelling of Epidemic Dynamics

Authors: Suraj Verma; M. A. Razzaque; U. Sangtongdee; C. Arpnikanondt; B. Tassaneetrithep; D. Arthan; C. Paratthakonkun; N. Soonthornworasiri
Journal: Computational and Mathematical Methods in Medicine
 (2021)

---

## Abstract

Hand, foot, and mouth disease (HFMD) is a highly contagious disease with several outbreaks in Asian-Pacific countries, including Thailand. With such epidemic characteristics and potential economic impact, HFMD is a significant public health issue in Thailand. Generally, contagious/infectious diseases’ transmission dynamics vary across geolocations due to different socioeconomic situations, demography, and lifestyles. Hence, a nationwide comprehensive model of the disease’s epidemic dynamics can provide information to understand better and predict a potential outbreak of this disease and efficiently and effectively manage its impact. However, there is no nationwide and comprehensive (i.e., the inclusion of reinfections in the model) model of HFDM dynamics for Thailand. This paper has endeavoured to promote nationwide comprehensive modelling of HFMD’s epidemic dynamics and comprehend the reinfection cases. We have formulated the SEIRS epidemiological model with dynamic vitals, including reinfections, to explore this disease’s prevalence. We also introduced periodic seasonality to reproduce the seasonal effect. The pattern of spread of this disease is uneven across the provinces in Thailand, so we usedK-means clustering algorithm to cluster those provinces into three groups (i.e., highly, moderately, and least affected levels). We also analysed health records collected from district hospitals, which suggest significant reinfection cases. For example, we found that 11% (approximately) of infectious patients return for repeat treatment within the study period. We also performed sensitivity analysis which indicates that the basic reproduction number (R0) is sensitive to the rate of transmission (β) and the rate at which infected people recover (γ). By fitting the model with HFMD confirmed data for the provinces in each cluster, the basic reproduction number (R0) was estimated to be 2.643, 1.91, and 3.246 which are greater than 1. Based on this high R0, this study recommends that this disease will persist in the coming years under identical cultural and environmental conditions.


---

## Body

## 1. Introduction

HFMD is a paediatric disease that became an epidemic in the Asian-pacific region by the 1990s. This disease has spread in several Asian-Pacific countries, including Australia, Cambodia, China, Japan, Malaysia, Taiwan, Vietnam, and Thailand [1]. In Thailand, the first outbreak was recorded in 2012 with 45,464 infected people which was the largest in a decade. Thailand experienced the second outbreak of this disease in 2016 [2]. With such frequent outbreaks, HFMD has been a critical health issue in Thailand and has a tremendous socioeconomic impact. In addition to children’s health, this disease affects their social life and education as they need to remain isolated from healthy children [3].Though HFMD is a paediatric disease, even parents are affected as they have to take special care of their children during this period. HFMD outbreaks have caused a considerable economic strain in the affected nation’s communities and healthcare systems.Epidemic control agencies in Thailand have encountered many violent outbreaks that resulted in deaths. Despite being a common disease every year, it may not be diagnosed by a physician due to unclear appearances of rash [3, 4]. Many severe symptoms are misdiagnosed with a similar disease—encephalitis, pneumonitis, pulmonary oedema, or myocarditis. The outbreak of this virus in Thailand can be due to indirect transmission by environment infection, asymptomatic infected people, and patients who lack medical care [4]. Although serotyping is an immune response, a person can get HFMD more than once if they are infected with a different serotype [5].Knowing a disease’s (e.g., HFMD) current spread and forecasting the epidemic are significant for healthcare and financial planning. Mathematical models of epidemiology explaining the population dynamics of contagious diseases have played a significant role in the better understanding of epidemiological patterns and control of diseases [6]. The transmission patterns of infectious diseases vary through geolocations due to various socioeconomic conditions, demographics, and people’s lifestyles [7]. Hence, a nationwide comprehensive model is required to understand the transmission of this disease.There is no research that estimates the reinfection rate of HFMD in Thailand to the best of our knowledge. Also, no nationwide comprehensive epidemiological model has been developed to estimate the basic reproduction number (R0) for Thailand. In this context, this research is aimed at developing a nationwide comprehensive epidemiological model of HFMD in Thailand. The objectives of this research are (i) to estimate the reinfection rate of HFDM in Thailand, (ii) to identify and cluster the affected provinces of Thailand based on the severity of infection, and (iii) develop a nationwide comprehensive epidemiological model with dynamic vitals for those clustered provinces. We estimated the recurrence of HFMD from the data collected from three local hospitals in Thailand, and this estimated recurrence rate was considered in the epidemiological model. HFMD has an incubation period in the host before they become infected. The incubation period is when an individual is exposed to the virus but not yet infected. Therefore, we considered the SEIRS (susceptible-exposed-infectious-recovered-susceptible) model in this paper, which addresses both the reinfection cases and the incubation period of the disease. This model describes the spread of HFMD and estimates the infection rate, recovery rate, and the basic reproduction number. Also, periodic seasonality was introduced in the model to address the seasonal outbreaks, which occur every rainy season. The data of hospitalised patients collected from BoE National Disease Surveillance (report 506) [2] was used to fit the proposed model and estimate the basic reproduction number (R0). This paper describes the simulation of the proposed model and performs sensitivity analysis to analyse input parameters’ impact in estimating R0.Section2 provides a brief literature review of related works conducted to estimate the basic reproduction number (R0) to understand HFMD disease transmission. Section 3 presents the materials and methods, including data collection and preprocessing, construction of the theoretical SEIRS model, and sensitivity analysis. The results from the study of the collected data are summarised in Section 4. This section also discussed the results obtained from clustering and fitting the clusters in the model using R2 to check the goodness of fit. In Section 5, we discussed the findings from the analysis and case study performed and the impact of COVID-19 in the spread of HFDM. Finally, we conclude the work with potential future research direction.


## 2. Related Works

In epidemiology, there is a long history of mathematical models. Many epidemiological models like SIR, SEIR, SEIQR, SEIRS, and many other models have been introduced to estimate infection rate of contagious diseases. The majority of these models are compartmental models, where the population is divided into compartments/groups with an assumption that the population transfers from one compartment to another [8]. Several mathematicians and researchers have proposed various epidemiological models to understand the transmission of HFMD better. Chuo et al. [9] used a SIR model to predict the number of infected people during the outbreak of HFMD in Sarawak of Malaysia. However, this model did not consider any seasonality effect and reinfection cases to spread this disease. Later on, a periodic model considered with quarantine effect was proposed by Liu [10]. Du et al. [7] adopted another approach to model this disease using time-series susceptible-infected-recovered model. In this model, a time scale for this disease was considered two weeks based on the characteristics of HFMD where the duration of the transition from infection to recovered and temporary immunity is about 7 to 10 days.Very few studies have been conducted for modelling the outbreak of this disease in Thailand. Chadsuthi and Wichapeng [4] researched the outbreak of HFMD in Bangkok in 2016 and proposed a customised SEIR model. The customised model considered susceptible (S), exposed but not yet infectious (E), infectious but with no symptoms (Ie), infectious with symptoms (I), hospitalised (H), and recovered (R). The basic reproduction number (R0) was predicted to be 1.441. However, the study conducted was during the outbreak of 2016 for Bangkok only, and no seasonal effect was considered in this model.Furthermore, in a study conducted by Huang et al. [11], it is found that multiple viruses, Enterovirus-71 (EV-71) and Coxsackievirus A16 (CA 16), are the causative agents of HFMD. Children infected with one virus are likely to be infected with other viruses that cause the same disease. As per their study conducted in China, 1.9% were reinfected at 12 months, 3.3% at 36 months, and 4.0% at 38.8 months. However, no study is taking care of the cases of reinfection in Thailand. Significantly, to date, no nationwide and comprehensive model has been developed for modelling the transmission of HFDM that considers seasonal effect and reinfection cases. CV-A16 and EV-A71 related to HFMD were analysed with an epidemic model in 2016 in China by [12]. Analysing time series then SIR modelling made this work having a model fit. Seasonality and heterogeneity were key features extracting from the province-level transmission. EV-A71 had a reproduction rate less than CV-A16, which was 26.63 and 27.13, respectively. Although their prediction was concentrated in sensitivity analysis, the primary purpose was to distribute vaccination effectively for herd immunity in a province, not to calculate the reinfection rate as we are proposing.


## 3. Materials and Methods

### 3.1. Data Exploration

For this study, we collected province-wise monthly reported data (from January 2011 to July 2020) of HFMD cases from BoE National Disease Surveillance (report 506) [2]. We also gathered data for HFMD cases for the different age groups of people in Thailand for 2019 from BoE. The population of children under 14 years for years 2019 was collected from the World Bank Group [13]. Moreover, to understand the reinfection cases, patients’ records were collected from three local hospitals: Samut Sakhon, Suan Pueng, and Ranong. These are district hospitals that offer standard medical services for Thai and non-Thai residents. The local datasets were derived from a retrospective study associated with hospital patients’ situation and trend in western borders between 2012 and 2016 [14]. All the personal identifiers were removed from the dataset to maintain the patients’ privacy. The Faculty of Tropical Medicine at Mahidol University gave ethical approval for the retrospective analysis that obtained local datasets.Data collection began with a selection of public-access health facilities. The sample involved the three public hospitals mentioned earlier. Service capabilities are similar to those of medical centres located in urban areas. For the study, we collected the patients’ records admitted to each hospital between 2012 and 2016. As HFMD is a paediatric disease, we separated the under fifteen age patients’ records for the study. The diagnostic statistics data was collected from the outpatient and inpatient medical record database. We used Equation (1) to identify the minimum sample size of the medical record required for our estimation [15].
(1)n=Z2pq/d2where
(i)
n is the minimum sample size needed for the study(ii)
Z is the level of confidence according to the standard normal distribution. For example, for a level of confidence of 95% or 5% level of significance (α=0.05), z=1.96(iii)
p is the prevalence of the disease in children(iv)
q is a subcalculation with 1−p, or percentage of failure(v)
d stands for precision limit or proportion of sampling error which is usually 5% confidence limitFrom the experiment, it was found that the minimum number sample required for the estimation is 384. Patient data were collected from the Medical Record Program for Enterovirus Infectious Disease. The disease codes B08.4, B08.5, and K12.0 were randomly sampled to obtain at least 384 records.


### 3.2. Data Preprocessing

The first step was preprocessing and cleaning the raw data to fit into the model. We performed this operation in two steps. First, we preprocessed the data to get HFMD cases across the provinces in Thailand, as shown in Figure1. For this purpose, we built a data model from raw data of Thailand population, HFMD cases, and Thailand’s provinces geocoordinates. The data was then cleaned by removing irrelevant, null, and empty data, and the cleaned data was transformed into the expected format. Finally, we extracted province-wise infected data and infection-wise health district zones. After getting the data for HFMD cases across different provinces in Thailand, the second step was to get the disease’s reinfection rate in Thailand. For this purpose, we processed the raw data collected from three local hospitals:; Ranong Hospital at the southern border, Samut Sakhon Hospital at the central area, and Rachaburi at the western border in Thailand Figure 2. We compiled all the data and filtered the HFMD cases’ data (ICD-10 code: B084). Then, patients’ details were deidentified to maintain their privacy. To extract the patients who revisited the hospital with HFMD complaint, patients with a record count≥2 were identified. Follow-up records were removed to get the revisit of HFMD cases but not a follow-up treatment.Figure 1
Data collection and preprocessing steps for extracting data to fit into the model.Figure 2
Showing the preprocessing steps for raw data collected from local hospitals in Thailand to estimate the reinfection rate of this disease. Those data were obtained from Ranong Hospital at the southern border, Samut Sakhon Hospital at the central area, and Rachaburi at the western border.


### 3.3. Model

We modelled the trajectory of the infection of HFMD disease in Thailand for years 2011 to 2019 using the SEIRS model with vital dynamics, a compartmental model. The incubation period of the disease is considered in this model [16]. This model also addresses the reinfection cases where a fraction of the recovered patients get reinfected with the same disease. The population is divided into four compartments, as shown in Figure 3. The compartments are represented by S, E, I, and R where
(i)
Susceptible (S) denotes people who are not infected but can get infected(ii)
Exposed (E) denotes people who are exposed but not yet infected(iii)
Infected (I) denotes people who are infected and can pass it on(iv)
Removed (R) denotes people who have recovered from this disease or diedFigure 3
Schematic representation of SEIRS epidemiological model with dynamic vitals illustrating the compartmental division of population.Let us assume the total population of children who can be infected with HFMD byN. Then, the total population can be represented as,N=S+E+I+R. The equations for the SEIRS model are as follows:
(2)dSdt=μN−S−βINS+ωR,dEdt=βINS−σ+μE,dIdt=σE−γ+μI,dRdt=γI−μ+ωR.Here,μ represents the birth and mortality rates per capita. Let us assume that the birth rate and the death rate are the same, giving the constant population throughout the year. β represents the contact/transmission rate at which infected individuals contact susceptible people to spread the disease. IfIis infected people, contactSis susceptible people, and the rate at which disease can spread isβ, then βI/N S represents the average number of infections. 1/σ and 1/γ represent the mean latent period and infectious period, respectively. The rate at which recovered patients become susceptible is represented by ω. The basic reproduction number (R0) used to measure the transmission potential of this disease using this model is calculated by Equation (3). At the same time, 1/β+μ+α gives the infectious period and the probability index of becoming infectious rather than dying. σ/σ+μ gives the probability index of becoming infectious while exposed (E).
(3)R0=βσ+μσγ+μ.


### 3.4. Seasonality of Transmission

HFMD in Thailand is a recurring epidemic disease. Historical data demonstrates that the disease outbreaks happen seasonally during the rainy season each year. The transmission rate varies seasonally giving a spike during the rainy season. The seasonal transmission was considered in our model, focusing on the transmission rateβ as β=βt. A sinusoidal trigonometric function was introduced to estimate β as follows:
(4)β=β01+sin2πt+αT,where β0, α, t, and T represent the baseline of transmission rate, a time constant, time, and period of seasonal cycle, respectively [17]. The main advantage of introducing this trigonometric function is to explore the qualitative characteristics of seasonally varying disease transmission and variable transmission rate.


### 3.5. Sensitivity Analysis

To identify the parameters that influence the basic reproduction number (R0), we performed a sensitivity analysis of the model. We used python’s SALib library to perform Sobol sensitivity analysis, a global sensitivity analysis [18]. First, the input parameters (beta, sigma, gamma, and mu) required to estimate the basic reproduction number to calculate the sensitivity were determined. Then, 1000 parameter samples were generated using the satellite library. Finally, the sensitivity index was calculated using the library’s analysis function on the output function using samples. Three different types of sensitivity indices, first-order, second-order, and total-order index, were calculated. The first-order index (S1) measures the contribution of the individual input parameters to R0. The second-order index (S2) measures a contribution of the combination of two parameters to R0. Similarly, the total-order index (ST) measures the contribution of both first- and higher-order indices on R0. Table 1 provides the details of the first-order and total-order sensitivity index of parameters on R0. It is seen that the first-order and total-order indices for beta and gamma are much higher compared to those for sigma and mu, indicating that these parameters are highly influencing parameters for the estimation of R0. Table 2 provides the details of the second-order sensitivity index for the combination of parameters. Here, a combination of beta and gamma has a higher value for the second-order index, indicating that these two parameters are crucial for calculating R0. Figure 4 illustrates the influence of each parameter on the basic reproduction number (R) for 1000 samples. As seen in the figure, we can determine that beta and gamma directly relate to the estimation of R0. For a grouped bar graph, Figure 5 was plotted to analyse the sensitivity of parameters. From this graph, it is seen that beta and gamma are critical parameters for R0. Although there is no cut-off value of sensitivity index, we considered the sensitivity index value>0.05 as an important parameter for estimating R0. In the figure, the total-order index (ST) value greater than the first-order index (S1) value indicates that there are higher-order interactions.Table 1
First-order sensitivity index of input parameters (beta, sigma, gamma, and mu) on the basic reproduction number (R0).
ParamS1S1 confSTST confBeta0.4678410.0584180.6030280.063202Sigma0.0060100.0095480.0139320.006417Gamma0.3550640.0544300.4931390.055598Mu0.0089480.0112500.0184980.004292Table 2
Second-order sensitivity index of input parameters (beta, sigma, gamma, and mu) on the basic reproduction number (R0).
Parameter pairsS2S2 confBeta/sigma0.0016890.073271Beta/gamma0.1318620.109820Beta/mu0.0057060.072474Sigma/gamma0.0038140.013194Sigma/mu0.0047950.012031Gamma/mu0.0082360.071884Figure 4
Sensitivity analysis scatter plots for 1000 sample input parameters.Figure 5
Sobol sensitivity indices for four parameters, beta, sigma, gamma, and mu, are shown in this figure. Total-order sensitivity index (ST) and first-order sensitivity index (S1) are represented in orange and blue colours, respectively. Beta and gamma are critical parameters for the estimation ofR0.


### 3.6. Data Fitting and Mathematical Simulation

This study appliesK-means, a self-adaptive weight clustering algorithm, to identify or cluster the provinces that are highly, moderately, and least affected with this disease to better estimate the basic reproduction number for those clustered provinces. The data collected from the BoE for the provinces in each cluster was fitted to simulate the SEIRS model. The initial values for S0 were taken as 3604225, 6563387, and 779287 for provinces in clusters 1, 2, and 3, respectively. Similarly, initial values for E0 and R0 were considered 0. The time of each cycle (T) was considered 12 months. The reinfection rate (ω=0.119) estimated from the three local hospitals’ dataset was also considered in our model. We used SciPy optimiser’s curve fitting tool to fit the model with the data and estimate the parameters. The curve fitting tool uses the least-squares minimisation technique to fit the data in the model. The goodness of fit of the simulated model was estimated using the determination coefficient (R2).


## 3.1. Data Exploration

For this study, we collected province-wise monthly reported data (from January 2011 to July 2020) of HFMD cases from BoE National Disease Surveillance (report 506) [2]. We also gathered data for HFMD cases for the different age groups of people in Thailand for 2019 from BoE. The population of children under 14 years for years 2019 was collected from the World Bank Group [13]. Moreover, to understand the reinfection cases, patients’ records were collected from three local hospitals: Samut Sakhon, Suan Pueng, and Ranong. These are district hospitals that offer standard medical services for Thai and non-Thai residents. The local datasets were derived from a retrospective study associated with hospital patients’ situation and trend in western borders between 2012 and 2016 [14]. All the personal identifiers were removed from the dataset to maintain the patients’ privacy. The Faculty of Tropical Medicine at Mahidol University gave ethical approval for the retrospective analysis that obtained local datasets.Data collection began with a selection of public-access health facilities. The sample involved the three public hospitals mentioned earlier. Service capabilities are similar to those of medical centres located in urban areas. For the study, we collected the patients’ records admitted to each hospital between 2012 and 2016. As HFMD is a paediatric disease, we separated the under fifteen age patients’ records for the study. The diagnostic statistics data was collected from the outpatient and inpatient medical record database. We used Equation (1) to identify the minimum sample size of the medical record required for our estimation [15].
(1)n=Z2pq/d2where
(i)
n is the minimum sample size needed for the study(ii)
Z is the level of confidence according to the standard normal distribution. For example, for a level of confidence of 95% or 5% level of significance (α=0.05), z=1.96(iii)
p is the prevalence of the disease in children(iv)
q is a subcalculation with 1−p, or percentage of failure(v)
d stands for precision limit or proportion of sampling error which is usually 5% confidence limitFrom the experiment, it was found that the minimum number sample required for the estimation is 384. Patient data were collected from the Medical Record Program for Enterovirus Infectious Disease. The disease codes B08.4, B08.5, and K12.0 were randomly sampled to obtain at least 384 records.


## 3.2. Data Preprocessing

The first step was preprocessing and cleaning the raw data to fit into the model. We performed this operation in two steps. First, we preprocessed the data to get HFMD cases across the provinces in Thailand, as shown in Figure1. For this purpose, we built a data model from raw data of Thailand population, HFMD cases, and Thailand’s provinces geocoordinates. The data was then cleaned by removing irrelevant, null, and empty data, and the cleaned data was transformed into the expected format. Finally, we extracted province-wise infected data and infection-wise health district zones. After getting the data for HFMD cases across different provinces in Thailand, the second step was to get the disease’s reinfection rate in Thailand. For this purpose, we processed the raw data collected from three local hospitals:; Ranong Hospital at the southern border, Samut Sakhon Hospital at the central area, and Rachaburi at the western border in Thailand Figure 2. We compiled all the data and filtered the HFMD cases’ data (ICD-10 code: B084). Then, patients’ details were deidentified to maintain their privacy. To extract the patients who revisited the hospital with HFMD complaint, patients with a record count≥2 were identified. Follow-up records were removed to get the revisit of HFMD cases but not a follow-up treatment.Figure 1
Data collection and preprocessing steps for extracting data to fit into the model.Figure 2
Showing the preprocessing steps for raw data collected from local hospitals in Thailand to estimate the reinfection rate of this disease. Those data were obtained from Ranong Hospital at the southern border, Samut Sakhon Hospital at the central area, and Rachaburi at the western border.


## 3.3. Model

We modelled the trajectory of the infection of HFMD disease in Thailand for years 2011 to 2019 using the SEIRS model with vital dynamics, a compartmental model. The incubation period of the disease is considered in this model [16]. This model also addresses the reinfection cases where a fraction of the recovered patients get reinfected with the same disease. The population is divided into four compartments, as shown in Figure 3. The compartments are represented by S, E, I, and R where
(i)
Susceptible (S) denotes people who are not infected but can get infected(ii)
Exposed (E) denotes people who are exposed but not yet infected(iii)
Infected (I) denotes people who are infected and can pass it on(iv)
Removed (R) denotes people who have recovered from this disease or diedFigure 3
Schematic representation of SEIRS epidemiological model with dynamic vitals illustrating the compartmental division of population.Let us assume the total population of children who can be infected with HFMD byN. Then, the total population can be represented as,N=S+E+I+R. The equations for the SEIRS model are as follows:
(2)dSdt=μN−S−βINS+ωR,dEdt=βINS−σ+μE,dIdt=σE−γ+μI,dRdt=γI−μ+ωR.Here,μ represents the birth and mortality rates per capita. Let us assume that the birth rate and the death rate are the same, giving the constant population throughout the year. β represents the contact/transmission rate at which infected individuals contact susceptible people to spread the disease. IfIis infected people, contactSis susceptible people, and the rate at which disease can spread isβ, then βI/N S represents the average number of infections. 1/σ and 1/γ represent the mean latent period and infectious period, respectively. The rate at which recovered patients become susceptible is represented by ω. The basic reproduction number (R0) used to measure the transmission potential of this disease using this model is calculated by Equation (3). At the same time, 1/β+μ+α gives the infectious period and the probability index of becoming infectious rather than dying. σ/σ+μ gives the probability index of becoming infectious while exposed (E).
(3)R0=βσ+μσγ+μ.


## 3.4. Seasonality of Transmission

HFMD in Thailand is a recurring epidemic disease. Historical data demonstrates that the disease outbreaks happen seasonally during the rainy season each year. The transmission rate varies seasonally giving a spike during the rainy season. The seasonal transmission was considered in our model, focusing on the transmission rateβ as β=βt. A sinusoidal trigonometric function was introduced to estimate β as follows:
(4)β=β01+sin2πt+αT,where β0, α, t, and T represent the baseline of transmission rate, a time constant, time, and period of seasonal cycle, respectively [17]. The main advantage of introducing this trigonometric function is to explore the qualitative characteristics of seasonally varying disease transmission and variable transmission rate.


## 3.5. Sensitivity Analysis

To identify the parameters that influence the basic reproduction number (R0), we performed a sensitivity analysis of the model. We used python’s SALib library to perform Sobol sensitivity analysis, a global sensitivity analysis [18]. First, the input parameters (beta, sigma, gamma, and mu) required to estimate the basic reproduction number to calculate the sensitivity were determined. Then, 1000 parameter samples were generated using the satellite library. Finally, the sensitivity index was calculated using the library’s analysis function on the output function using samples. Three different types of sensitivity indices, first-order, second-order, and total-order index, were calculated. The first-order index (S1) measures the contribution of the individual input parameters to R0. The second-order index (S2) measures a contribution of the combination of two parameters to R0. Similarly, the total-order index (ST) measures the contribution of both first- and higher-order indices on R0. Table 1 provides the details of the first-order and total-order sensitivity index of parameters on R0. It is seen that the first-order and total-order indices for beta and gamma are much higher compared to those for sigma and mu, indicating that these parameters are highly influencing parameters for the estimation of R0. Table 2 provides the details of the second-order sensitivity index for the combination of parameters. Here, a combination of beta and gamma has a higher value for the second-order index, indicating that these two parameters are crucial for calculating R0. Figure 4 illustrates the influence of each parameter on the basic reproduction number (R) for 1000 samples. As seen in the figure, we can determine that beta and gamma directly relate to the estimation of R0. For a grouped bar graph, Figure 5 was plotted to analyse the sensitivity of parameters. From this graph, it is seen that beta and gamma are critical parameters for R0. Although there is no cut-off value of sensitivity index, we considered the sensitivity index value>0.05 as an important parameter for estimating R0. In the figure, the total-order index (ST) value greater than the first-order index (S1) value indicates that there are higher-order interactions.Table 1
First-order sensitivity index of input parameters (beta, sigma, gamma, and mu) on the basic reproduction number (R0).
ParamS1S1 confSTST confBeta0.4678410.0584180.6030280.063202Sigma0.0060100.0095480.0139320.006417Gamma0.3550640.0544300.4931390.055598Mu0.0089480.0112500.0184980.004292Table 2
Second-order sensitivity index of input parameters (beta, sigma, gamma, and mu) on the basic reproduction number (R0).
Parameter pairsS2S2 confBeta/sigma0.0016890.073271Beta/gamma0.1318620.109820Beta/mu0.0057060.072474Sigma/gamma0.0038140.013194Sigma/mu0.0047950.012031Gamma/mu0.0082360.071884Figure 4
Sensitivity analysis scatter plots for 1000 sample input parameters.Figure 5
Sobol sensitivity indices for four parameters, beta, sigma, gamma, and mu, are shown in this figure. Total-order sensitivity index (ST) and first-order sensitivity index (S1) are represented in orange and blue colours, respectively. Beta and gamma are critical parameters for the estimation ofR0.


## 3.6. Data Fitting and Mathematical Simulation

This study appliesK-means, a self-adaptive weight clustering algorithm, to identify or cluster the provinces that are highly, moderately, and least affected with this disease to better estimate the basic reproduction number for those clustered provinces. The data collected from the BoE for the provinces in each cluster was fitted to simulate the SEIRS model. The initial values for S0 were taken as 3604225, 6563387, and 779287 for provinces in clusters 1, 2, and 3, respectively. Similarly, initial values for E0 and R0 were considered 0. The time of each cycle (T) was considered 12 months. The reinfection rate (ω=0.119) estimated from the three local hospitals’ dataset was also considered in our model. We used SciPy optimiser’s curve fitting tool to fit the model with the data and estimate the parameters. The curve fitting tool uses the least-squares minimisation technique to fit the data in the model. The goodness of fit of the simulated model was estimated using the determination coefficient (R2).


## 4. Results

We analysed the collected data to understand the spread of HFMD in different provinces throughout the year. We found that the disease outbreak generally occurs in June, July, and August of each year throughout Thailand. Two outbreaks were seen in 2012 and 2016, where the number of infected cases during the rainy season was very high compared to other seasons (Figure6). We plotted a heat map for HFMD cases to visualise the spread of this disease across provinces during the years 2011 to 2019. From the heat map, we found that, according to the severity of infections, the provinces could be grouped into three clusters (Figure 7). Thus, the K-means algorithm was applied to cluster the provinces based on the number of infection cases from 2011 to 2019. Based on the severity of infection of the disease as observed in the heat map and K-means clustering result, we identified three different clusters of provinces in Thailand. The provinces are clustered into highly, moderately, and least affected provinces, as shown in Figure 8 and Table 3.Figure 6
Illustration of HFMD confirmed cases from the year 2011 to 2019 throughout Thailand. This figure shows that there is a sharp peak in the number of HFMD cases during the rainy season of each year.Figure 7
Heat map to illustrate HFMD cases from 2011 to 2019 across all the provinces in Thailand. From the figure, we can visualise that the provinces can be grouped into three groups as per the spread of this disease. Bangkok was seen to be the most affected province.Figure 8
Geomapping with clustering of provinces based on the number of infections.Table 3
The provinces of Thailand are clustered in 3 clusters based on the number and severity of infections.
ClustersInfection severityProvincesCluster 1LowLampang, Lamphun, Mae Hong Son, Nan, Phayao, Phrae, Phetchabun, Sukhothai, Tak, Uttaradit, Chai Nat, Kamphaeng Phet, Phichit, Uthai Thani, Angthong, Nakhon Nayok, Nonthaburi, P.Nakhon S.Ayutthaya, Pathum Thani, Sing Buri, Kanchanaburi, Nakhon Pathom, Phetchaburi, Prachuap Khiri Khan, Ratchaburi, Samut Sakhon, Samut Songkhram, Suphan Buri, Chachoengsao, Chanthaburi, Prachin Buri, Sa Kaeo, Trat, Kalasin, Mahasarakham, Roi Et, Bungkan, Loei, Nakhon Phanom, Nong BuaLam Phu, Nong Khai, Sakon Nakhon, Udon Thani, Chaiyaphum, Amnat Charoen, Mukdahan, Yasothon, Chumphon, Krabi, Nakhon Si Thammarat, Phangnga, Phuket, Ranong, Narathiwat, Pattani, Phatthalung, Satun, Songkhla, Trang, YalaCluster 2HighBangkokCluster 3AverageChiang Mai, Chiang Rai, Phitsanulok, Nakhon Sawan, Lop Buri, Saraburi, Chon Buri, Rayong, Samut Prakan, Khon Kaen, Buri Ram, Nakhon Ratchasima, Surin, Si Sa Ket, Ubon Ratchathani, Surat ThaniTo analyse the age group of infected people in Thailand, we plotted a bar graph for 2019. As expected, we observed that children under 14 years are more infected than older people (Figure9).Figure 9
Age group of HFMD infected 2019.The severity of infection of this disease across the provinces in Thailand was visualized by plotting a heat map over the provincial map of Thailand (Figure10). We used “Paintmaps” online tool to plot the heat map over Thailand’s provincial map and visualise the geographical spread of this disease for the year 2019. This geographical map illustrates that HFMD has spread in most of the provinces, and Bangkok is the most affected province. From the case study performed at three different local hospitals in Thailand, we found that a significant number of patients who visited the hospital were reinfected with HFMD. For the selected three hospitals (2012-2016), we found different cumulative reinfection cases or rates which are not insignificant, and it is increasing with the inclusion of an additional year (Table 4). For example, the reinfection rate or percentage for the year 2012-2015 is 9.39%, which has changed to 11.92%, adding the year 2016.Figure 10
Shade maps of Thailand entirely with 77 provinces. This map illustrates the severity of HFMD infection spreading and reported over the country in 2019.Table 4
Reinfection cases were observed at three local hospitals in Thailand.
YearReinfection casesPatientsPercentage2012143833.662012 to 2013426966.032012 to 20149111517.902012 to 201514715659.392012 to 2016249208911.92


### 4.1. Fitting the Model

SEIRS model was fitted using curve fitting, and the goodness of fit of the model to the recorded data for each cluster is illustrated in Figures11–13 for provinces in clusters 1, 2, and 3, respectively. The model fitted the HFMD models for the years 2011 to 2019. For better fitting, the model was fitted for each year, and values for the parameters α, β0, σ, μ, γ were estimated.Figure 11
SEIRS model for provinces in cluster 1 for the years 2011-2019.Figure 12
SEIRS model for provinces in cluster 2 for the years 2011-2019.Figure 13
SEIRS model for provinces in cluster 3 for years 2011-2019.


### 4.2. Estimation of Parameters

For the model simulation, the value ofN was considered the population of children under 14 years. The initial infected population, I0, was taken from the BoE. The parameters α, β0, σ, μ, and γ were estimated by fitting the recorded data in the model using curve fit. Tables 5–7 illustrate the estimated parameters for the provinces in clusters 1, 2, and 3, respectively. The tables also show the calculated basic reproduction number R0.Table 5
Estimation of parameters for provinces in cluster 1.
Parameters201120122013201420152016201720182019N6563521I0134E00————————T12ω0.119α-0.41-0.5770.0540.3371.578-0.2560.483-0.083-0.575β03.0611.8351.9642.3791.5392.0072.0762.2042.186σ0.73940.0182.1753.0150.7111.8182.8363.854.62γ3.0830.9190.4690.631-0.0240.9280.3840.7160.762μ-0.3311.2610.811.0670.3981.1510.9961.0331.065R015.2613.8771.8292.9840.3693.8022.1213.0383.246Table 6
Estimation of parameters for provinces in cluster 2.
Parameters201120122013201420152016201720182019N779326I039E00————————T12ω0.119α-0.512-0.7810.8790.8131.438-0.0011.1150.998-0.317β04.7494.8211.4881.2160.681.3881.5221.6731.812σ1.10374.4380.7614.6834.25329.9752.1341.4367.488γ5.8613.539-0.623-0.052-0.8250.552-0.111-0.0990.624μ-0.4772.980.9470.8461.0660.8770.7970.7171.036R045.0730.2210.2150.8180.1311.9270.7610.6892.643Table 7
Estimation of parameters for provinces in cluster 3.
Parameters201120122013201420152016201720182019N3604332I0107E00————————T12ω0.119α-0.6330.2640.9360.271-0.435-0.1941.104-0.177-0.092β02.8020.6340.8251.8962.0221.6321.5722.161.907σ0.62283.552.9075.7281.75223.5592.2082.9163.074γ2.2790.34-0.5780.6160.3680.674-0.3820.730.384μ-0.1740.1310.8911.0330.8661.0541.0550.840.917R08.1890.2980.1982.6491.672.6990.7152.6321.91


### 4.3. Validity of Model

The goodness of fitting of the model was tested by regression analysis. The determination coefficientR2 was calculated using the following formulae (5), (6), and (7):
(5)SStot=∑iy1−y∧2,(6)SSres=∑iy1−f12,(7)R2=1−SSresSStotThe standard of model fitting has been evaluated for each cluster independently by calculatingR2 using Equation (7). The result for R2 obtained for each cluster model is illustrated in Table 8. The values are greater than 0.6 which illustrates that at least 60% of the recorded data fitted well in the model. Out of three clusters, for cluster 2, the model fitted the recorded data very well with R2 77.7%.Table 8
Calculation ofR2.
ClustersR2Cluster 10.63Cluster 20.777Cluster 30.756


## 4.1. Fitting the Model

SEIRS model was fitted using curve fitting, and the goodness of fit of the model to the recorded data for each cluster is illustrated in Figures11–13 for provinces in clusters 1, 2, and 3, respectively. The model fitted the HFMD models for the years 2011 to 2019. For better fitting, the model was fitted for each year, and values for the parameters α, β0, σ, μ, γ were estimated.Figure 11
SEIRS model for provinces in cluster 1 for the years 2011-2019.Figure 12
SEIRS model for provinces in cluster 2 for the years 2011-2019.Figure 13
SEIRS model for provinces in cluster 3 for years 2011-2019.


## 4.2. Estimation of Parameters

For the model simulation, the value ofN was considered the population of children under 14 years. The initial infected population, I0, was taken from the BoE. The parameters α, β0, σ, μ, and γ were estimated by fitting the recorded data in the model using curve fit. Tables 5–7 illustrate the estimated parameters for the provinces in clusters 1, 2, and 3, respectively. The tables also show the calculated basic reproduction number R0.Table 5
Estimation of parameters for provinces in cluster 1.
Parameters201120122013201420152016201720182019N6563521I0134E00————————T12ω0.119α-0.41-0.5770.0540.3371.578-0.2560.483-0.083-0.575β03.0611.8351.9642.3791.5392.0072.0762.2042.186σ0.73940.0182.1753.0150.7111.8182.8363.854.62γ3.0830.9190.4690.631-0.0240.9280.3840.7160.762μ-0.3311.2610.811.0670.3981.1510.9961.0331.065R015.2613.8771.8292.9840.3693.8022.1213.0383.246Table 6
Estimation of parameters for provinces in cluster 2.
Parameters201120122013201420152016201720182019N779326I039E00————————T12ω0.119α-0.512-0.7810.8790.8131.438-0.0011.1150.998-0.317β04.7494.8211.4881.2160.681.3881.5221.6731.812σ1.10374.4380.7614.6834.25329.9752.1341.4367.488γ5.8613.539-0.623-0.052-0.8250.552-0.111-0.0990.624μ-0.4772.980.9470.8461.0660.8770.7970.7171.036R045.0730.2210.2150.8180.1311.9270.7610.6892.643Table 7
Estimation of parameters for provinces in cluster 3.
Parameters201120122013201420152016201720182019N3604332I0107E00————————T12ω0.119α-0.6330.2640.9360.271-0.435-0.1941.104-0.177-0.092β02.8020.6340.8251.8962.0221.6321.5722.161.907σ0.62283.552.9075.7281.75223.5592.2082.9163.074γ2.2790.34-0.5780.6160.3680.674-0.3820.730.384μ-0.1740.1310.8911.0330.8661.0541.0550.840.917R08.1890.2980.1982.6491.672.6990.7152.6321.91


## 4.3. Validity of Model

The goodness of fitting of the model was tested by regression analysis. The determination coefficientR2 was calculated using the following formulae (5), (6), and (7):
(5)SStot=∑iy1−y∧2,(6)SSres=∑iy1−f12,(7)R2=1−SSresSStotThe standard of model fitting has been evaluated for each cluster independently by calculatingR2 using Equation (7). The result for R2 obtained for each cluster model is illustrated in Table 8. The values are greater than 0.6 which illustrates that at least 60% of the recorded data fitted well in the model. Out of three clusters, for cluster 2, the model fitted the recorded data very well with R2 77.7%.Table 8
Calculation ofR2.
ClustersR2Cluster 10.63Cluster 20.777Cluster 30.756


## 5. Discussion

The major finding of this study is the reinfection rate of HFMD in Thailand. HFMD-recovered patients become immune to the virus, and the chances of getting reinfected with the same virus is scarce. However, as the disease is caused by multiple viruses (Enterovirus-71 (EV-71) and Coxsackievirus A16 (CA 16)) [19], there is a chance of getting infected with other virus causing the same disease. From the result of the case study performed at three local hospitals in Thailand, it is found that 11.92% (approximately) of the recorded cases are reinfection cases. The case of reinfection is significantly high, especially during the rainy seasons.HFMD is most prevalent spreading in younger ages compared to the older groups of patients across the country. During the mid of each year, the total number of patients increases consistently in a similar previous pattern. One of the reasons is that the rainy season in Thailand starts from early May to late July and this disease transmits quite rapidly during the rainy season. Moreover, primary and secondary schools begin the first semester in the middle of May [20], which leads to an increase in physical contact between infected and susceptible children. At the beginning of the school term, children are not provided with full support, especially wellness and health checks and physical education. These may cause the infection to reach a peak during the mid of each year. The reinfection is also the supporting cause of outbreak during rainy seasons each year.Figure6 illustrates that the number of infected cases is very high during rainy seasons and the outbreak of this disease occurred in 2012 and 2016. The heat map (Figure 7) and geographical map (Figure 10) help us to visualise the spread of this disease across all the provinces. We found that all provinces are not equally affected. Bangkok is observed to be the most affected province. Since this disease’s transmission is uneven across the provinces, it is essential to identify the provinces with similar disease severity. Hence, we utilised the K-means algorithm to cluster the provinces based on the severity of the infection cases. The resulting cluster grouped the provinces into three categories, highly, moderately, and least affected provinces. The main advantage of clustering the provinces based on the number of infections is that positively affected provinces can get more attention and the healthcare facilities and human resources can take appropriate actions to control the disease.Mathematical models have always been an essential tool for understanding the spread of epidemiological diseases. No comprehensive model has been developed to address this disease’s nationwide transmission with reinfection cases in Thailand. To understand and predict the spread of this disease, we have built a comprehensive compartmental SEIRS model with vital dynamics by incorporating the BoE surveillance data and reinfection rate estimated from the three hospitals to estimate theR0 of HFMD. We have also introduced periodic seasonality to address the seasonal effect in preventing this disease. Generally, the time from initial infection to the occurrence of signs and symptoms (incubation period) of HFDM is 3 to 7 days [19]. During this period, infected people are not yet infectious.According to the proposed model,E (exposed) addresses the incubation (latent) period. This model also addresses the rate of reinfection where a fraction of the recovered patients become susceptible because of different virus type. ω is the rate at which recovered patients become susceptible. The proposed SEIRS model was fitted with recorded HFMD cases based upon each cluster’s provinces to better estimate the parameters. Figures 11–13 illustrate the goodness of fit of recorded data for HFMD cases in provinces in clusters 1, 2, and 3, respectively. The model fitted well the outbreaks of this disease seen in the rainy season of each year. While fitting the model, α, β0, σ, μ, and γ parameters were estimated which are illustrated in Tables 5–7 for provinces in three clusters. The basic reproduction number (R0) obtained from the simulation of three clustered provinces is greater than 1, indicating that this disease will persist in the coming years under the same social and environmental condition.HFMD is transmitted primarily through direct contact with contaminated discharge, saliva, or stool from infected patients or contaminated objects. Moreover, attending kindergarten or child-care centres or schools increases the risk factor for the transmission of this disease [21]. Avoiding physical interaction with infected people, such as kissing, embracing, and sharing cups and eating utensils, plays a significant role in minimising the spread of this disease. Because of the COVID-19 pandemic, a nationwide lockdown and quarantine of infected patients were imposed. The schools were closed, and no physical contact was made with infected patients which played significant role in reducing HFMD cases in 2020. For each month of this year, the number of infected cases was decreasing. The predicted number of HFMD-infected patients and the recorded number of cases for the year 2020 for each cluster are illustrated in Figures 14(a)–14(c). These figures illustrate a significant decrease in the number of infected cases. Even during the rainy season, the cases were decreasing because of the impact of COVID-19 lockdown imposed in the year 2020.Figure 14
These figures illustrate the predicted and actual HFMD cases in 2020. It is seen that the actual HFMD cases in each clustered province are significantly low compared to the prediction.
(a)
Provinces in cluster 1(b)
Provinces in cluster 2(c)
Provinces in cluster 3This study estimated the HFMD reinfection rate from limited data collected from three local hospitals in Thailand. This estimate could have been improved by collecting more and random samples from hospitals from different regions of the country or the Bureau of Epidemiology (BoE). However, due to the COVID pandemic, we could not collect sufficient data. Our model considers the reinfection rate for the estimation of the transmission rate (β), the recovery rate (γ), and the basic reproduction number (R0), so an improved reinfection rate estimate can improve the estimation of β, γ, and R0. Also, HFMD-infected patients may not show any symptoms but can still spread the disease. Such infections are asymptomatic infections. Many HFMD-infected patients are hospitalised, which reduces the physical contact between infected and susceptible patients. Thus, asymptomatic infections and hospitalisation cases are essential factors for estimating the basic reproduction number. However, due to the limited available data, these factors have not been considered in our model.


## 6. Conclusions

This paper has estimated the recurrence rate of HFMD disease in Thailand. We have applied the estimation to develop a nationwide and comprehensive model of epidemic dynamics using the SEIRS model. With dynamic vitals, including reinfections, our results indicate that this disease’s transmission is most prevalent in midyear. From the experiment performed using the model and study conducted at three provincial-level hospitals, we can conclude that this seasonal disease is a critical health issue.The spread of this disease is uneven throughout the provinces across the country. Hence, we clustered the provinces into three clusters based on the number of infections from the years 2011 to 2019 and the basic reproduction number (R0) obtained from the SEIRS model for all these clusters is greater than 1, which indicates that under the same social and environmental condition, this disease will persist in the coming years.Moreover, we learned that the reinfection rate, especially during the rainy season, is remarkably high in Thailand. Such circumstance is one of the main reasons for higher infection rates or cases during the rainy season. However, the number of infected patients is significantly low for 2020 due to the COVID-19 pandemic. The lockdown imposed on account of the pandemic minimises direct or indirect physical contact with infected people and spread of COVID-19. This measure also helped to minimise the spread of HFMD. So, taking proper and informed actions could improve the management of HFMD-like diseases in Thailand and globally. Informed actions need real-time or quasi-real-time data about infected people. Hence, our future endeavours will focus on smartphone-based digital diagnosis and reporting of HFMD.



---

*Source: 6697522-2021-03-05.xml*