Complex interaction between meteorological factors on the risk of hand, foot, and mouth disease

Hand, foot, and mouth disease (HFMD) is a contagious disease. Large, severe outbreaks affecting thousands of people occurred frequently particularly in the Asian-Pacific regions such as China (Hong et al 2022), Singapore (Min et al 2021), Japan (Gonzalez et al 2019), India (Rao et al 2017), Malaysia (Lee et al 2021), and New Zealand (Hayman et al 2014). Usually, HFMD is self-limiting, but a small proportion of children may experience severe complications such as meningitis, encephalitis, acute flaccid paralysis and neurorespiratory syndrome, leading to death. Therefore, HFMD has become one of the leading public health concerns in the world.

The relationship between meteorological factors and the risk of HFMD has been widely studied. For example, Kim et al found that the HFMD rate increased by 10.3% and 6.6% for every unit rise in average temperature and relative humidity below the cutoffs of 18 °C and 65%, and then decreased as temperature and humidity increased (2016). Gou et al found that the increase in weekly HFMD incidence were associated with every unit increase in average temperature (2.8% increased in HFMD incidence), and relative humidity (2.47%) in Lanzhou (2018). Abdul Wahid et al found that the risk of HFMD increased with rainfall below 60 mm and decreased with rainfall exceeding 60 mm (2021). One of our previous studies also showed an increase in risk up to 23.5 °C temperature and decreased passed this point while risk was increased with increasing relative humidity (Du et al 2019a).

This evidence suggests a potential adverse impact of multiple meteorological factors on the risk of HFMD. The mechanism behind the relationship between meteorological factors and HFMD remains unsure. Several proposed factors, including pathogen infectivity outside the host, human behavior, and immune function fluctuations, could be involved in the risk of HFMD (Fares 2013). For example, during the suitable weather time, individuals are more likely to have outdoor activities, thereby it may increase the chances that susceptible individuals contact with infectious individuals or contaminated environment. In addition, according to the previous studies in the laboratory, the reproduction and survival rates of the enterovirus that cause HFMD are strongly affected by fluctuations in environment (Arita et al 2005).

Although a plenty of evidence on the association between meteorological exposures and the risk of HFMD was available, significant gaps remain in our understanding of interactions among these factors on the occurrence of HFMD. Such knowledge would provide a better understanding of the meteorological exposure periods that lead to the highest risk of HFMD and improve predictions in the temporal trend for HFMD epidemics. However, the existing evidence is limited. Zhang et al assessed the interactive effects between temperature and humidity by including an indicator of high/low temperature/humidity in models but did not find any significant estimates (2016). This study was limited by it is statistical techniques, given that converting two continuous variables into a single categorical variable would increase uncertainty in the results due to the loss of information. In our previous studies, we used classification and regression tree models to evaluate the association of HFMD incidence with temperature, humidity, and air pressure, while allowing the interactions among meteorological factors to be examined (Du et al 2016). Risk of HFMD was higher when the temperature was greater than 24 °C and the relative humidity was less than 61% such as during non-autumn seasons. However, our study was only able to classify HFMD risk by meteorological factors with cutoff values and could not consider continuous interaction effects of meteorological factors (i.e. regression coefficients or RRs). In addition, only a very few interactions between meteorological factors have been considered in the available studies. Given the highly correlated nature of meteorological factors, it is imperative to include multiple interaction terms. Furthermore, health impacts of meteorological events are likely to be nonlinear, leading to an increased complexity of the models.

To address these methodology gaps and add new evidence for the interactive relationship among meteorological factors, we used Bayesian kernel machine regression (BKMR) to investigate interactive effects among the possible meteorological factors on the risk of HFMD while minimizing the uncertainties that arise in traditional models. In addition, we explored the city-specific factors that were attributed to the heterogeneity of interactive effects between cities on HFMD.

2.1. Study area

2.2. Case data

2.3. Meteorological data

2.4. Data of potential effect modifiers

2.5. Statistical analysis

The kernel machine regression (KMR) model (Bobb et al 2018) was used to estimate the risk of HFMD as a function of meteorological exposures. The major advantage of KMR is that it allows the nonlinear impact of continuous exposures and more importantly the interactive effects between any two exposures while holding all other factors constant at a specific value. We specified our model as follows:

$$\begin{align}\log \left( {{\text{cas}}{{\text{e}}_t}} \right) & = h\left( {{\text{facto}}{{\text{r}}_{t1}}, \ldots {\text{facto}}{{\text{r}}_{tM}}} \right) \nonumber\\ &\quad + \beta^*\left( {\cos \left( {\frac{{2\pi t}}{{52}}} \right) + \sin \left( {\frac{{2\pi t}}{{52}}} \right)} \right) + {\varepsilon _t}\end{align} \tag{ 1 }$$

where case_t is the weekly number of HFMD cases, t is the index of week, i.e. 1,...52 × 9 (a total of 9 years), there were weeks with 0 incidence cases, which we replaced with 0.5 to avoid errors in the logarithmic function, h denotes the nonlinear association and interaction between meteorological factors on the risk of HFMD to be estimated, m = 1,...,6 (a total of six meteorological factors), β (β∼ N(0, σ_β ²)) is the corresponding vector of coefficients of the trigonometric functions 'cos((2πt)/52) + sin((2πt)/52)' for modeling seasonality (time trend) as recommended (Meyer and Held 2017), and the residual ∼ N(0, σ ²) is assumed to be independent and identically (iid) normally distributed with a common variance.

We used the Gaussian kernel function to flexibly capture the various forms of h as follows:

$$\begin{align}& h\left( {{\text{facto}}{{\text{r}}_i},{\text{facto}}{{\text{r}}_j}} \right) \nonumber\\ &\quad = {\text{exp}}\left\{ { - {\text{ }}\sum\limits_{m = 1}^M {{r_m}} \left( {{\text{facto}}{{\text{r}}_{im}} - {\text{facto}}{{\text{r}}_{jm}}} \right)} \right\}\end{align} \tag{ 2 }$$

where factor_im and factor_jm represent two different values of the m meteorological factor. r_m denotes the tuning parameter that controls the smoothness of h as a function of factor_m . In particular, we defined univariate exposure as the risk ratios (RRs) of HFMD over the range of a certain exposure while all the other meteorological factors were fixed at their median value. Similarly, we defined bivariate effect as the RRs of HFMD over the range of the first meteorological factor while the second factor was fixed at the 10th, 50th, and 90th percentiles and the rest of factors were fixed at their median value.

We fitted the KMR model in a Bayesian framework using Markov Chain Monte Carlo with 10 000 iterations. The prior distributions for the unknown parameters were specified as follows:

$$\begin{align}\sigma _\beta ^2,\sigma _\varepsilon ^2\sim {\text{Gamma}}\left( {0.001,0.001} \right)\end{align} \tag{ 3 }$$

$$\begin{align}{r_m}\sim {\text{uniform}}\left( {0,0.01} \right).\end{align} \tag{ 4 }$$

To assess the robustness of our findings which could be particularly sensitive to r_m , we performed a sensitivity analysis by specifying a geometric sequence of r_m values to fit separate models to each of the meteorological factors as recommended (Bobb et al 2018). In addition, we assessed model convergence using the trace plots of the variance of the residuals. All the trace plots showed a relatively constant mean and variance, indicating that the model has well converged (figure S11).

We examined the interaction between meteorological factors using the heterogeneity test in meta-analysis. Specifically, for each city, we estimated RRs per each unit increase in the first meteorological factor at its 50th percentile while the second meteorological factor was fixed at its 10th, 50th, and 90th percentile values (therefore, three RRs were obtained), respectively. For each RR, we then pooled city-specific estimates using meta-analysis based on linear mixed-effect model. The Cochran's Q statistic (which follows the Chi-square distribution) and its corresponding P value were obtained from the meta-analysis and were used to determine the significance of the interactive effects. A larger value of Q statistic represents a greater interactive effect. The relationships between the Q statistic and the city-specific factors were examined using the linear regression to identify differences in interaction effects among cities.

All statistical analyses were performed using R version 4.1.3 (the R Foundation for Statistical Computing, Vienna, Austria), including the packages of bkmr, metafor, and ggplot2. There was no missing value in these data.

Figure S1 shows the long-term time trends of HFMD cases and meteorological factors in Guangzhou, in a time-series within a 1 year period. An intuitive relationship can be observed given that the peaks of temperature, humidity, rainfall, and sunshine duration occur at approximately the same time as the peaks in the HFMD cases. However, the peak occurrence times of pressure and wind speed were the inverse of the HFMD cases. Similar results were also observed in Shantou, Zhanjiang, and Meizhou (figures S2–S4).

Figure 1 shows the univariate effect of meteorological factors on the risk of HFMD in Guangzhou. An inverted V-shaped relationship was observed between the risk of HFMD and temperature, with the peak at about 62.5th percentile (22.0 °C). An inverted U-Shaped relationship was observed between the risk of HFMD and wind speed, with the peaks at about 62.5th percentile (3.2 m s⁻¹). An almost linear relationship was observed between the risk of HFMD and pressure, rainfall, and sunshine duration. At the 50th percentile, the RRs were 1.773(95%CI: 1.458, 2.156) for temperature, 1.055(0.874, 1.274) for humidity, 1.013(0.739, 1.391) for rainfall, 1.017(0.871, 1.183) for pressure, 1.238(0.981, 1.558) for wind speed, and 1.016(0.878, 1.174) for sunshine duration, respectively (table 1). Similar results were observed in Meizhou and Shantou (tables S2, S3, figures S5 and S6), while Zhanjiang had its own characteristics (table S4 and figure S7).

Zoom In Zoom Out Reset image size

Figure 2 shows a bivariate relationship between two meteorological factors on the risk of HFMD in Guangzhou. As displayed in the first row of panels, with increased levels of humidity [10th (red curve) vs 50th (green curve) vs 90th (blue curve)], the risk of HFMD associated with the other meteorological factors including pressure, rainfall, sunshine duration, temperature, and wind speed also increased, with a curve for a greater percentile of humidity generally being higher than that for a smaller one. Similarly, we found that the pressure and rainfall had less interactive effects with other meteorological factors as limited disparities between the curves were observed. Based on our estimates, the sunshine duration and wind speed tended to have similar quantitative interaction as humidity with other meteorological factors. In contrast, temperature tended to modify the association between the risk of HFMD and other meteorological factors but not monotonously. These nonmonotonic interactions for meteorological factors indicated a higher level of effect on the risk of HFMD with the median temperature than the other temperature. Similar results were observed in Meizhou and Shantou (figures S8 and S9), while Zhanjiang had weaker interactions (figure S10).

Zoom In Zoom Out Reset image size

Table 2 presents the RRs of HFMD in Guangzhou for each meteorological factor fixed at its median value while considering an additional factor at its 10th, 50th, and 90th percentiles. In the case of temperature and humidity, the RRs were 1.621(95%CI: 1.226, 2.141), 2.638(2.169, 3.208) and 3.734(2.908, 4.792) at median temperature when the humidity was at its 10th, 50th, 90th percentiles, respectively (P_interaction< 0.001). Alternatively, the RRs were 0.901(95%CI: 0.592, 1.369), 2.026(1.679, 2.448) and 0.884(0.632, 1.238) at median humidity when the temperature was at its 10th, 50th, 90th percentiles, respectively (P_interaction< 0.001). For humidity and wind speed, the RRs for humidity held at its 50th percentile value were 0.621(95%CI: 0.491, 0.787), 0.854(0.708, 1.032) and 1.309(1.015, 1.688) when the wind speed was at its 10th, 50th, 90th percentiles, respectively (P_interaction< 0.001). The RRs of wind speed when held at its 50th percentile value were 1.162(0.860, 1.563), 1.929(1.529, 2.423), and 2.766(2.017, 3.785) when the humidity was at its 10th, 50th, 90th percentiles, respectively (P_interaction< 0.001). Similar results were observed in Meizhou and Shantou (tables S5 and S6), while Zhanjiang was observed with fewer significant interactions (table S7).

The interactive effects between humidity and other meteorological factors were greater in Guangzhou than the other three cities (table S8). For temperature and humidity, the heterogeneity statistics (Cochran's Q) were 19.132, 0.265, 1.741, and 2.118 for Guangzhou, Shantou, Meizhou, and Zhanjiang, respectively. GDP per capita was positively associated with the interactive effects between humidity and other meteorological factors including sunshine duration (P= 0.013), rainfall (P= 0.017), temperature (P= 0.021), win speed (P= 0.011), and pressure (P= 0.013) (tables 3 and S9). In addition, the number of per capita passenger traffic were also positively associated with the interactive effects between humidity and other meteorological factors including sunshine duration (P< 0.001), rainfall (P= 0.001), temperature (P= 0.002), wind speed (P< 0.001), and pressure (P< 0.001). Furthermore, we found increased per capital area of parks and green land was also associated with greater interactive effects between humidity and other meteorological factors including sunshine duration (P= 0.042), rainfall (P= 0.035), temperature (P= 0.031), wind speed (P= 0.045), and pressure (P= 0.042).

We found greater interactive effects of temperature with other meteorological factors in Shantou than the other three cities. For temperature and humidity, the Q statistics for the interaction between humidity and temperature were 24.876, 69.877, 31.785, and 24.217 for Guangzhou, Shantou, Meizhou, and Zhanjiang, respectively. The annual sunshine duration was positively associated with the interactive effects between temperature and other meteorological factors including rainfall (P = 0.013), humidity (P = 0.022), and wind speed (P = 0.024). While increased rate of sewage treatment was associated with greater interactive effects of temperature with other meteorological factors including sunshine duration (P = 0.014) and pressure (P = 0.015).

We also conducted multiple sensitivity analyses. First, to verify the reliability of the results of the inter-city heterogeneity analysis, we doubled the number of cities for sensitivity analysis. The results of the sensitivity analysis found that sunshine duration, GDP per capita, rate of sewage treatment, and per capita area of parks and green land in city were dominant predictors underlying the inter-city heterogeneity, which was consistent with the results of the main analysis (table S10). Second, we further used the multivariable linear regression instead of the KMR to explored the interactions between meteorological factors in Guangzhou. However, only 3 out of 15 interaction terms were statistically significant (table S11).

We found that increased humidity may be associated with a greater effect of temperature, rainfall, pressure, wind speed, and sunshine duration on the risk of HFMD. For humidity and temperature, one of our previous studies reported a relatively high risk of HFMD when the temperature was greater than 24.85 °C and the relative humidity was above 80.59% (Du et al 2016). Similar results were also observed for other infectious diseases. For example, Liu et al reported an interaction between temperature and relative humidity on the risk of influenza in Shaoyang, China. The cold temperatures and lower humid conditions were associated with the increased risks (Liu et al 2019). The potential mechanisms underlying the interaction may be that humid conditions improve the viability and spread efficiency of the virus by providing suitable circumstances for transmission (Lowen and Palese 2009). Particularly, the humid environment was more suitable for the reproduction and spread of enteroviruses which is a major pathogen of HFMD.

We found that a moderate rather than low or high temperature strengthened the effect of humidity, rainfall, pressure, wind speed, and sunshine duration on the risk of HFMD. Although there were few studies on temperature modifying the effects of other meteorological factors, most studies on temperature and HFMD reported that moderate levels of temperature correspond to the highest risk. Consistent with the findings of many previous studies, the current findings and our previous studies have found an inverted V-shaped temperature-HFMD relationship (Du et al 2019a). Qi et al also found an inverted V-shaped relationship between temperature and HFMD (2018). There is a potential that at the peak of the inverted V-shaped relationship between temperature and HFMD, the median temperatures, which is the most suitable for outdoor activities, increased the infectivity of enterovirus (Averett 2016). These temperatures along with other meteorological conditions, such as humidity, increase the likelihood of contact between infected and susceptible people by increasing their social network (Du et al 2019b).

In addition to the humidity and temperature, we also found that increased wind speed could strengthen the effect of sunshine duration, rainfall, humidity, temperature, and pressure on the risk of HFMD. Without considering the interaction effects, Ma et al found that wind speed was positively associated with HFMD consultation rates (2010). Liao et al reported a positive correlation between HFMD incidence in children and wind speed (2015). The potential reason for an increase in HFMD with increased wind speed could be that wind can transport pollutants and water droplets containing these viruses, accelerating the propagation velocity of droplets, and therefore spreading the HFMD virus. This mechanism can also explain the role of wind speed in promoting the risk on HFMD associated with other meteorological factors.

We found that the cities with greater per capita GDP and area of parks and green land were suffering from greater interactive effects between humidity and other meteorological factors (e.g. temperature). Without assessing the relationship between meteorological factors, Gao et al (2021) and Li et al (2018) also found that GDP was positively associated with increased HFMD incidence. This association was potentially due to susceptible people in cities with a higher GDP having more access to resources (e.g. raincoats or umbrella) and likely had more options for outdoor spaces at different humidity levels. The per capital area of parks and green land may increase risk of transmission due to shared space among susceptible populations, and prior studies by our group found that the green space would increase the contact rate and increased the risk of HFMD (Du et al 2022).

Cities with higher annual hours of sunshine duration and lower rate of sewage treatment were found to have statistically significant interactive effects of temperature on other meteorological factors (e.g. humidity). Again, without considering the relationship between other meteorological factors, van Pham et al also found that the risk of HFMD increased with the increase in hours of sunshine (2019). The potential reason the number of hours of sunshine could impact HFMD is that susceptible people may enjoy more outdoor activities in days with favorable temperature and increased hours of sunshine. In regard to the rate of sewage treatment, Ji et al found the HFMD viral pathogens were highly prevalent in the raw sewage (86%) and treated wastewater (29%) (2012). Thus, the lower rate of sewage treatment would result in high background levels of environmental exposure, thereby increasing the risk of infection in susceptible populations.

While our study had many strengths, including the large sample size (a total of 468 weeks in 9 years) and the new insight for complex interactions between meteorological factors, there were also some limitations. First, as with other ecological studies, we cannot determine a causal relationship between meteorological factors and the risk of HFMD because the time series study is essentially an ecological design. Second, HFMD infections may have mild symptoms or be asymptomatic resulting self-management of health rather than seeking medical care. Although the cases reported to the surveillance system could be assumed to cover the almost entirely infections, underestimation is unavoidable. Third, the possible higher-order interactions were not analyzed in this study due to limitation of the existing methodologies. Future studies would be needed to clarify such higher-order interactions.

Our study contributed further understanding of complex interactions between meteorological factors on the risk of HFMD. It appears that good weather leads to higher infection rates of HFMD. This finding may be related to improved viability of the virus under these conditions as well as increased social interactions between infected and uninfected people leading to viral spread. Our findings provide epidemiological evidence for meteorological interactions on HFMD, which may guide health authorities to integrate various meteorological conditions to precisely identify high-risk days and thus make timely and effective preventive and control strategies.

This work was supported by the National Natural Science Foundation of China (Nos. 82103947 and 81973150) and the Natural Science Foundation of Guangdong, China (No. 2022A1515010823).

The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.

This study was approved by the Institutional Review Board at the School of Public Health, Sun Yat-sen University (Approval number: 2021-018). All patient information of HFMD surveillance data were collected from China Center for Disease Control and Prevention (China CDC). No confidential information was included and all methods were carried out in accordance with relevant guidelines and regulations. All data employed in this study were deidentified prior to analysis and anonymized confidentiality, thus informed consent was not required.

The authors declare no competing interests.