Analysis of HFMD Transmissibility Among the Whole Population and Age Groups in a Large City of China

Background: Hand-Foot-and-Mouth-Disease (HFMD) has been widely spread in Asia, and has result in a high disease burden for children in many countries. However, the dissemination characteristics intergroup and between different age groups are still not clear. In this study, we aim to analyze the differences in the transmissibility of HFMD, in the whole population and among age groups in Shenzhen city, by utilizing mathematical models.

Methods: A database that reports HFMD cases in Shenzhen city from January 2010 to December 2017 was collected. In the first stage, a Susceptive-Infected-Recovered (SIR) model was built to fit data of Shenzhen city and its districts, and R_eff was used to assess transmissibility in each district. In the second stage, a cross-age groups SIR model was constructed to calculate the difference in transmissibility of reported cases among three age groups of EV71 virus: 0–3 years, 3–5 years, and over 5 years which was denoted as age group 1, 2, and 3, respectively.

Results: From 2010 to 2017, 345,807 cases of HFMD were reported in Shenzhen city, with peak incidence in spring and autumn in Shenzhen city and most of its districts each year. Analysis of the EV71 incidence data by age group revealed that age Group 1 have the highest incidence (3.13 ×10⁻⁷–2.31 ×10⁻⁴) while age group 3 had the lowest incidence (0–3.54 ×10⁻⁵). The differences in weekly incidence of EV71 between age groups were statistically significant (t₁₂ = 7.563, P < 0.0001; t₂₃ = 12.420, P < 0.0001; t₁₃ = 16.996, P < 0.0001). The R² of the SIR model Shenzhen city population-wide HFMD fit for each region was >0.5, and P < 0.001. R_eff values were >1 for the vast majority of time and regions, indicating that the HFMD virus has the ability to spread in Shenzhen city over the long-term. Differences in R_eff values between regions were judged by using analysis of variance (ANOVA) (F = 0.541, P = 0.744). S_iI_iR_i-S_jI_jR_j models between age groups had R² over 0.7 for all age groups and P <0.001. The R_eff values between groups show that the 0–2 years old group had the strongest transmissibility (median: 2.881, range: 0.017–9.897), followed by the over 5 years old group (median: 1.758, range: 1.005–5.279), while the 3–5 years old group (median: 1.300, range: 0.005–1.005) had the weakest transmissibility of the three groups. Intra-group transmissibility was strongest in the 0–2 years age group (median: 1.787, range: 0–9.146), followed by Group 1 to Group 2 (median: 0.287, range: 0–1.988) and finally Group 1 to Group 3 (median: 0.287, range: 0–1.988).

Conclusion: The incidence rate of HFMD is high in Shenzhen city. In the data on the incidence of EV71 in each age group, the highest incidence was in the 0–2 years age group, and the lowest incidence was in the over 5 years age group. The differences in weekly incidence rate of EV71 among age groups were statistically significant. Children with the age of 0–2 years had the highest transmissibility.

Background

Hand-Foot-and-Mouth Disease (HFMD) has been a health concern in Asia since the late 1990s (1), with outbreaks reported in Malaysia, Japan, Singapore, Vietnam, and Cambodia (2–4). HFMD has a high disease burden in Asia, with 358,764 cases of HFMD in Japan each year (5). Also, it remains a major health problem for children in China, affecting more than 2 million children each year. Especially, Guangdong Province has been one of the most severely affected provinces (6). The annual incidence of HFMD in Guangdong Province exceeds 30 cases in every 10,000 people and case number has accounted for approximately 15% of the total number of cases in China in recent years (7). In addition, Shenzhen city had the highest number of HFMD cases among 143 cities in mainland China between 2009 and 2014 (8).

HFMD is a viral disease that has received great attention in the last two decades because of its high incidence rate in the pediatric population. The clinical presentation of HFMD is characterized by fever and a blistering rash, mostly on the hands, feet and oral mucosa (9). The disease is generally mild and self-limiting, but neurologic and cardiopulmonary-related complications may occur in disease outbreaks (10). Among all common febrile and rash illnesses (11), HFMD has remained the one that most frequently affects young children (12, 13). A series of analyses have shown that children younger than 5 years old are more likely to develop HFMD and that it is more likely to accumulate in younger children (14–16).

The main pathogens of HFMD infection are enterovirus 71 (EV71) and coxsackievirus A16 (CVA16) (14, 17), and EV71 infection is particularly the main viral subtype causing severe and fatal cases (12).

However, the mechanism by which EV71 causes severe central nervous system complications is still unclear. It is suggested that this may be due to a combination of pathologic immune response and direct viral action, but there is no effective treatment and no biomarker that can be used as an early warning of severe HFMD. In terms of vaccines, Asian countries that have experienced a history of HFMD infection and pandemics are actively developing vaccines as a preventive measure. Among them, a vaccine for EV-A71 in China has been available in 2016 (9, 10), and a bivalent EV-A71 / CV-A16 vaccine should enter clinical trials in the near future (13). The research of Joseph T. Wu et al. showed that the cost of routine EV71 vaccination is cost-effective in China (15). However, it is promoted as a class II vaccine in China, causing some limitations in terms of immunization coverage.

Currently, many mathematical models have also been applied to the study of HFMD, such as the application of Bayesian spatio-temporal models to analyze the factors influencing HFMD and the effects of interventions (16, 18, 19), the application of lagged nonlinear models to assess the effects of meteorological factors on incidence (20, 21). SIR models (22–24), SEIR models (25), SEILR models (26), and other models (27, 28) are also used in the study of HMFD. In addition, some scholars have also focused on the differences in HMFD between different age groups (29, 30). It indicates that the existing studies related to HFMD still focus mostly on spatial and temporal variation (11, 31), the expression of overall epidemiological characteristics and the relationship between meteorological factors and the incidence of HFMD, while there is a lack of studies targeting the analysis of EV71 virus by age groups.

In this study, we assumed that there is variability in the transmissibility among the major pathogens of HFMD and variability and possible interaction in the transmission of EV71 virus among different age groups. Based on this hypothesis, we collected data on reported HFMD cases in a city with a high prevalence of HFMD in China (Shenzhen city), grouped the cases by age according to their basic information, and then used a seasonally adjusted SIR model to calculate transmissibility of different pathogens among different age groups, and further investigated the magnitude of all-virus transmissibility of HFMD in Shenzhen city. Besides, in this paper, we studied for the first time the characteristics of the transmission pattern of HFMD among different age groups in the whole population, established a transmission model of EV71 virus transmission among different age groups, and analyzed the reasons for the differences between groups.

Materials and Methods

Data Collection

In this study, we collected HFMD cases from January 2010 to December 2017 in Shenzhen city, Guangdong Province, including the number of reported cases and deaths per day, age, sex, exposure history, date of onset, and severity of disease, and obtained basic information on the population of Shenzhen city, including the year-end resident, birth rate and death rate, by searching the Shenzhen city Statistical Yearbook 2010–2017. Between 2010 and 2017, the population exceeded 12 million, and the average birth rate was 18.63 per 1,000 (range: 16.41/1,000–21.46/1,000) and the median death rate was 7.32 per 1,000 (range: 6.19/1,000–9.72/1,000). When we come to the Population Composition of Shenzhen city, the number of people in each age group increased linearly from 2014 to 2017, but the composition ratio remained stable: the 0–2 years old group remained stable around 2.611% (range: 2.609–2.613%), the 3–5 years old group remained stable around 2.163% (range: 2.161–2.165%) and the over 5 years old group remained stable around 95.225% (range: 95.222–95.230%).

The population of the 0–2 years old group was 281.240 thousand, 297.156 thousand, 310.928 thousand and 327.383 thousand from 2014 to 2017 respectively. The population of the 3–5 years old group from 2014 to 2017 was 232.886 thousand, 246.202 thousand, 271.202 thousand and 271.202 thousand respectively. The population of the over 5 years old group from 2014 to 2017 was 1024.789 thousand, 10,835.555 thousand, 11,929.736 thousand and 11,929.736 thousand respectively.

Shenzhen city is a large city located in the south of the coastal area of Guangdong Province, China. Shenzhen city has a southern subtropical monsoon climate with long summers and short winters, mild climate, abundant sunshine and rainfall, and an annual average temperature of 23.0°C. The climate is strongly influenced by the monsoon, with “warm temps” and humid weather in spring, high temperatures and rain in summer, little rain in fall with drought and often typhoons, and short and dry winters with little rain (32). Overall, it is a city with long summers and short winters, strongly influenced by the monsoon climate.

Model Introduction

SIR Model

Based on the transmission mechanism of HFMD and the type of case report data, the Susceptible -Infectious -Removed (SIR) model was used (Figure 1).

FIGURE 1

Figure 1. The Susceptive-Infected-Recovered model.

In this model, S denotes susceptible individuals, I denotes infectious individuals, R denotes recovered individuals and N denotes the total population size. The parameters br, dr, f , β and γ refer to the natural birth rate of the population, the mortality rate of the population, the mortality rate of HFMD, the relative transmission rate and the relative recovery rate, respectively.

$\begin{array}{l}\frac{dS}{dt}=brN-\beta SI-drS\\ \frac{dI}{dt}=\beta SI-\gamma I-(dr+f)I\\ \frac{dR}{dt}=\gamma I-drR\\ N=S+I+R\end{array}$

S_iI_iR_i-S_jI_jR_j Model

Considering the interactions between age groups, we created the age group SIR model (Figure 2).

FIGURE 2

Figure 2. The cross-age groups Susceptive-Infected-Recovered model.

We set 0–2 years old as age group 1; 3–5 years old as age group 2; and over 5 years old as age group 3.

$\begin{array}{l}\frac{d{S}_1}{dt}=brN-{\beta }_{11}{S}_1{I}_1-{\beta }_{21}{S}_1{I}_2-{\beta }_{31}{S}_1{I}_3\\ \frac{d{I}_1}{dt}={\beta }_{11}{S}_1{I}_1+{\beta }_{21}{S}_1{I}_2+{\beta }_{31}{S}_1{I}_3-{\gamma }^{\prime }{I}_1\\ \frac{d{R}_1}{dt}={\gamma }^{\prime }{I}_1\end{array}$

$\begin{array}{l}\frac{d{S}_2}{dt}=-{\beta }_{22}{S}_2{I}_2-{\beta }_{12}{S}_2{I}_1-{\beta }_{32}{S}_2{I}_3\\ \frac{d{I}_2}{dt}={\beta }_{22}{S}_2{I}_2+{\beta }_{12}{S}_2{I}_1+{\beta }_{32}{S}_2{I}_3-{\gamma }^{\prime }{I}_2\\ \frac{d{R}_2}{dt}={\gamma }^{\prime }{I}_2\end{array}$

$\begin{array}{l}\frac{d{S}_3}{dt}=-{\beta }_{33}{S}_3{I}_3-{\beta }_{13}{S}_3{I}_1-{\beta }_{23}{S}_3{I}_2-dr{S}_3\\ \frac{d{I}_3}{dt}={\beta }_{33}{S}_3{I}_3+{\beta }_{13}{S}_3{I}_1+{\beta }_{23}{S}_3{I}_2-{\gamma }^{\prime }{I}_3-dr{I}_3\\ \frac{d{R}_3}{dt}={\gamma }^{\prime }{I}_3-dr{R}_3\end{array}$

$\begin{array}{l}N=S_1+S_2+S_3+I_1+I_2+I_3+R_1+R_2{+R}_3\end{array}$

Parameter Estimation

As shown in Table 1, there are eight parameters (β₁, γ, br, dr, and f ).

1) The propagation coefficient among individuals in the SIR model is set to β. β = 0.340609 is obtained by fitting the reported data of Shenzhen city.

2) The parameters br, dr, and f (33) were calculated based on finding the Shenzhen city Yearbook and the collected data. The annual br and dr values were collected; the weekly values of the two parameters were calculated; and the weekly br and dr values were 0.000352 (range: 0.000330–0.000383) and 0.0000129 (range: 0.0000127–0.0000187). The references give the susceptibility of the population and the likelihood of the population and the infectiousness of the population were 0.67, 0.4, and 0.433, respectively.

3) The age group SIR model in which f₁, f₂, and f₃ denote the disease and death rates in the three age groups 0–2 years, 3–5 years, and over 5 years, respectively, were taken as 0.0003. β₁₁ denotes the transmission coefficient among individuals in age group 1; β₂₂ denotes the transmission coefficient among individuals in age group 2 minus; β₁₂ denotes the transmission coefficient from individuals in age group 1 to individuals in age group 2; β₁₃ denotes the transmission coefficient from individuals in age group 1 to individuals in age group 2; and β₁₄ denotes the transmission coefficient from individuals in age group 1 to individuals in age group 2. Propagation coefficient from individuals within age group 1 to individuals in age group 3; β₂₁ denotes the propagation coefficient from individuals within age group 2 to individuals in age group 1; β₂₃ denotes the propagation coefficient from individuals within age group 2 to individuals in age group 3; β₃₁ denotes the propagation coefficient from individuals within age group 3 to individuals in age group 1; and β₃₂ denotes the propagation coefficient from individuals within age group 3 to individuals in age group 2. The R_eff subscript corresponds to the same age group object.

4) According to the seasonal cyclic pattern of the incidence data, they were divided into 18 segments, and β were fitted separately within and between each age group.

TABLE 1

Table 1. Parameter estimation.

Indicators of Transmission Capacity

The basic reproduction number is an important parameter to determine whether the disease is epidemic or not, which refers to the number of new cases expected to be directly transmitted by one infectious agent in the susceptible population during its transmission period. Effective reproduction number (R_eff) was set as the evaluation index, i.e., the basic reproduction number after intervention, to evaluate the impact of the intervention on the relative transmission capacity of HFMD among the population. The calculation process of the R_eff equation can be seen whereupon:

Step 1: Divide the derivatives of all infected compartments into two parts: the first part represents F new infections and the other part ${V}$ represents the transformation between non-new compartments.

$\begin{array}{l}\frac{d}{dt}\left(\begin{array}{c}{I}_1\\ I_2\\ I_3\end{array}\right)=\left(\begin{array}{c}{S}_1=\left(I_1{\beta }_{11}+I_2{\beta }_{21}+I_3{\beta }_{31}\right)\\ S_2=\left(I_1{\beta }_{12}+I_2{\beta }_{22}+I_3{\beta }_{32}\right)\\ S_3=\left(I_1{\beta }_{13}+I_2{\beta }_{23}+I_3{\beta }_{33}\right)\end{array}\right)-\left(\begin{array}{c}{I}_1\gamma \\ I_{2}dr+I_2\gamma \\ I_{3}dr+I_3\gamma \end{array}\right)\\ \stackrel{def}{=}{F}-{V}\end{array}$

Step 2: Taking the derivatives of the vectors ${F}$ and ${V}$ with respect to the variable [I₁ I₂ I₃], respectively, yields the corresponding Jacobi matrices ${F}$ and ${V}$ :

$\begin{array}{l}{F}=\left(\begin{array}{ccc}{S}_1{\beta }_{11}& S_1{\beta }_{21}& S_1{\beta }_{31}\\ S_2{\beta }_{12}& S_2{\beta }_{22}& S_2{\beta }_{32}\\ S_3{\beta }_{13}& S_3{\beta }_{23}& S_3{\beta }_{33}\\ ‏\end{array}\right)\\ {V}=\left(\begin{array}{ccc}\gamma & 0& 0\\ 0& dr+\gamma & 0\\ 0& 0& dr+\gamma \\ ‏\end{array}\right)\end{array}$

and further calculate the ${V}$ ⁻¹:

$\begin{array}{l}{{V}}^{-1}=\left(\begin{array}{ccc}\frac{1}{\gamma }& 0& 0\\ 0& \frac{1}{dr+\gamma }& 0\\ 0& 0& \frac{1}{dr+\gamma }\\ ‏\end{array}\right)\end{array}$

Step 3: Calculate the next generation matrix M = ${F}$ ${V}$ ⁻¹. The elements in row j and column i of the matrix M are the R_{eff, ij.} of the intergroup interactions. The result of the calculation of R_ij can be obtained by replacing S_i=N_i into the equation.

$\begin{array}{l}M=\left(\begin{array}{ccc}\frac{S_1{\beta }_{11}}{\gamma }& \frac{S_1{\beta }_{21}}{dr+\gamma }& \frac{S_1{\beta }_{31}}{dr+\gamma }\\ \frac{S_2{\beta }_{12}}{\gamma }& \frac{S_2{\beta }_{22}}{dr+\gamma }& \frac{S_2{\beta }_{32}}{dr+\gamma }\\ \frac{S_3{\beta }_{13}}{\gamma }& \frac{S_3{\beta }_{23}}{dr+\gamma }& \frac{S_3{\beta }_{33}}{dr+\gamma }\\ ‏\end{array}\right)\end{array}$

The elements in row j and column i of the matrix M are the R_{eff, ij.} of the intergroup interactions. The result of the calculation of R_ij can be obtained by replacing S_i=N_i into the equation.

Step 4: R₀ is defined as the maximum eigenvalue of M, whereupon:

$\begin{array}{l}{R}_0={\lambda }_{max}(M)\end{array}$

The analytic expressions for the eigenvalues of the matrix M are particularly complex. The following approach is taken for this purpose.

1) First calculate the value of R_ij and then the maximum eigenvalue of the numerical matrix M.

2) Use the definition method to calculate the overall R₀.

$\begin{array}{l}{R}_0={\displaystyle \sum _{i=1}^n}\left[P(x_0\in N_i){\displaystyle \sum _{i=1}^n}{R}_{ij}\right]\end{array}$

Software Introduction

Berkeley Madonna 8.3.18 (developed by Robert Macey and George Oster of the University of California at Berkeley. Copyright ©1993–2001 Robert I. Macey & George F. Oster) ran model coefficients estimated by curve fitting to the data and simulated the effects of the intervention. SPSS 13.0 (IBM Corp., Armonk, NY, USA) was used for statistical testing, with t-tests used to calculate differences between age groups and R² used to assess curve fitting.

Results

Epidemiological Features

From 2010 to 2017, 345,807 cases of HFMD were reported in Shenzhen city, including 129,812 cases in Baoan District, 26,377 cases in Luohu District, 24,008 cases in Futian District, 23,009 cases in Nanshan District, 136,832 cases in Longgang District, and 5,769 cases in Yantian District. Among them, 233,895 cases were 0–2 years old, 98,911 cases were 3–5 years old, and 16,071 cases were over 5 years old. Among them, 16 cases died, and the mortality rate was 0.0059%.

Analyzing the incidence data of Shenzhen city from 2010 to 2017, the median incidence rate was 5.26/100,000 people. The lowest incidence rate was 0.11/100,000 people in week 4 of 2011; the highest incidence rate was 45.10/100,000 people in week 39 of 2017. From the spatio-temporal distribution map (Figure 3), it can be seen that the burden of HMFD in Shenzhen city gradually increased from 2010 to 2017, and the increase in incidence rate was mainly concentrated in Baoan and Longgang districts, among which Longgang district had a more obvious trend of rising first and then decreasing. In addition, after 2013, the burden of disease in Nanshan District, Futian District, Luohu District, and Yantian District differed significantly from Baoan District and Longgang District, and there were also slight fluctuations during the period.

FIGURE 3

Figure 3. The Map of HMFD burden distribution in Shenzhen city.

By analyzing the weekly reported case data, there are approximately two epidemic cycles per year, alternating seasonally from spring to summer and from summer to autumn. Except for 2013, Shenzhen city and most of its districts had peak incidence in both spring and autumn, especially in Luohu and Yantian districts, where the peak incidence in spring tended to be higher than that in autumn (Figure 4). Among the districts, Baoan, Longgang and Yantian had higher incidence rates, with Longgang having the highest incidence rate (1.28 × 10⁻²-7.53 × 10⁻⁴) and Futian having the lowest incidence rate (0–2.75 × 10⁻⁴). In 2017, Baoan district had the highest incidence rate in Shenzhen city, with an incidence rate of 87.36/100,000 people. The differences in incidence rates among districts were statistically significant as calculated by ANOVA (F = 65.006, P < 0.0001). And the results of multiple comparisons in different districts were showed in Table 2.

FIGURE 4

Figure 4. Fitting results for HMFD incidence rates in Shenzhen city and its districts.

TABLE 2

Table 2. The P values of multiple comparisons for different districts (LSD).

Analysis of the incidence data by age group also revealed a peak in spring and autumn, with a significantly higher peak in spring than in autumn. Incidence rates were the highest (3.13 × 10–7 to 2.31 × 10–4) in group 1 (the 0–2 years age group) and the lowest (0–3.54 × 10–5) in group 3 (the over 5 years age group), with the incidence in group 1 being 10.02 times higher than that in the group 3 (Figure 5). The differences in weekly EV71 prevalence between age groups were statistically significant (t₁₂ = 7.563, P < 0.0001; t₂₃ =12.420, P < 0.0001; t₁₃ =16.996, P < 0.0001;). The composition ratio of the three age groups did not change significantly during 2014–2017, and the median composition ratio of group 2 was 65.267% (59.654–74.922%), the median composition ratio of group 2 (the 3–5 years group) was 29.875% (20.385–33.895%), and the median composition ratio of group 3 was 5.197% (4.155–6.452%), all fluctuating within a certain range.

FIGURE 5

Figure 5. Fitting results for EV71 incidence rates by age group.

Model Fitting Effect

SIR Model

First, we fitted the SIR model for the weekly incidence of population-wide HFMD in Shenzhen city and each sub-district. In the fitting process, segments were fitted to the data. According to the seasonal characteristics of HFMD, each peak period was divided into a segment, with an average of 2–3 segments per year in Shenzhen city. The results of fitting for population-level HFMD in each district of Shenzhen city were good, with R² >0.5 and p < 0.001 for all districts (Figure 4). The quartiles of reported incidence in Shenzhen city were P₂₅ = 1.78 × 10⁻⁵, P₅₀ = 5.26 × 10⁻⁵, P₇₅ = 1.19 × 10⁻⁴. The total incidence reported in Shenzhen city had an increasing trend year by year, with the highest peak in 2017 (4.51 × 10⁻⁴), while the remaining districts showed a trend of increasing and then decreasing trend, except for Futian and Nanshan districts.

S_iI_iR_i-S_jI_jR_j Model

We also developed the S_iI_iR_i-S_jI_jR_j model to fit the HFMD incidence data between age groups in Shenzhen city for 0–2 years, 3-5 years, and over 5 years of age for transmissibility. As can be seen from the fitted curves (Figure 5), the fit was good among the age groups. The fitted R² values for each group are listed in the listed fit effect table, and it can be seen that R² is over 0.7 for all age groups, and p < 0.001 for all (Table 3). The quartiles for age group 1 were P₂₅ = 1.08 × 10⁻⁵, P₅₀ = 3.17 × 10⁻⁵, P₇₅ = 6.58 × 10⁻⁵, age group 2 were P₂₅ = 4.98 × 10⁻⁶, P₅₀ = 1.23 × 10⁻⁵, P₇₅ = 2.62 × 10⁻⁵, and age group 3 were P₂₅ = 8.26 × 10⁻⁷, P₅₀ = 2.26 × 10⁻⁶, P₇₅ = 4.97 × 10⁻⁶. The transmission coefficients β for each age group at each time period were obtained by fitting, and the fitting results are shown in Table 4.

TABLE 3

Table 3. The P values and fitted R² values for each group.

TABLE 4

Table 4. The weekly transmission coefficients β for each age group.

Assessment of Transmissibility

The calculated R_eff values of Shenzhen city and each district are shown in Figure 6, and it can be seen that the R_eff values for most of the time are >1, indicating that the HFMD virus were likely to continue spreading in the long term.

FIGURE 6

Figure 6. The transimissibility of the first stage.

The median R_eff of Shenzhen city was 4.001 (Interquartile range [IQR]: 0.021–8.929), with R_eff values <1 in the second cycle of 2011; the median R_eff for Baoan was 4.069 (IQR: 2.096–10.139), with R_eff values all over 1; the median R_eff of Longgang District was 4.126 (IQR: 2.935–8.922), with R_eff values all over 1 and a gradual upward trend from 2010 to 2016; the median R_eff of Luohu District had a median R_eff of 4.298 (IQR: 1.730- 9.314), with R_eff values over 1; the median R_eff of Nanshan District was 3.939 (IQR: 1.937- 8.002), with R_eff values all over 1; the median R_eff of Futian District was 4.777 (IQR: 2.924- 11.105), with a decreasing trend of R_eff over 1; the median R_eff in Yantian District was 5.131 (IQR: 1.730- 9.314), with R_eff over 1. R_eff values were slightly <1, indicating that HFMD will continue to be prevalent in Shenzhen city and its districts but will not cause a major outbreak. The differences in R_eff values between districts were not statistically significant (F = 0.541, P = 0.744) by Analysis of Variance (ANOVA).

The results of the cross-age groups SIR model calculations showed that the 0–2 years age group had the strongest transmissibility (Median[M]: 2.881, IQR: 0.017–9.897), followed by the over 5 years age group (M: 1.758, IQR: 1.005–5.279), while the 3–5 years age group (M:1.300, IQR:0.005–1.005) had the weakest transmissibility among the three groups.

The R_eff for the 0–2 years age group had the largest value in the first half of 2010, followed by a peak in late 2013 and then a gradual decrease; the R_eff for the 3–5 years age group fluctuated within a certain range and had the largest value in the second half of 2017. R_eff values for the over 5 years age showed an increasing and then decreasing trend, peaking at the end of 2013 (Figure 7A).

FIGURE 7

Figure 7. The transimissibility of the second stage (A) Group1 indicates the R_eff value of 0–2years old age group; Group 2 indicates the R_eff value of 3–5 years old age group; Group 3 indicates the R_eff value of over 5 years old age group; (B) Rij denotes the transmission coefficient from individuals in age group i to individuals in age group j, which both i and j take the value 1, 2, 3.

Intra-group transmissibility was strongest in the 0–2 years age group (R_eff: M = 1.787, IQR = 0–9.146), followed by that of Group 1 to Group 2 (R_eff12: M = 0.287, IQR = 0–1.988) and finally Group 1 to Group 3 (R_eff: M = 0.287, IQR = 0–1.988 median: 0.287, range: 0–1.988). The R_eff for this intra-group showed a significant decreasing trend during the survey years, and Group 1 to Group 2 showed a trend of increasing before decreasing R_eff and peaked in 2014.

The transmission in the 3–5 years age group was the highest in the intra-group (median: 0.627, range: 0.005–2.928) and Group 2 to Group 1(median: 0.497, range: 0-4.617), and the lowest in the Group 2 to Group 3(median: 0.060, range: 0-0.590); where the intra-group transmissibility of this age group showed a decreasing trend in R_eff year by year.

In the over 5 years age group, the R_eff from Group 3 to Group 1 was the largest and perennially over1 (median: 1.217, range: 1.000–2.619), followed by Group 3 to Group 2(median: 0.063, range: 0–3.412), while intra-group transmission (median: 0.104, range: 0–1.260) was the smallest; where the R_eff from Group 3 to Group 2 showed a trend of rising and then falling and peaked in the second half of 2015, and the intra-group transmission of this age group rose yearly until 2014 peaked and then fell and tended to zero (Figure 7B).

Discussion

Analysis and observation of the reported incidence of HFMD in Shenzhen city revealed that the incidence rate of HFMD in Shenzhen city is higher than the national-level incidence rate (34, 35), which may be due to the fact that Shenzhen city is a coastal city with a hot and humid climate and a high population density, conditions that favor the reproduction and spread of the virus. In addition, the incidence of HFMD in Asian subtropical countries is characterized by distinct peaks in spring and autumn (36), and the incidence data reported in Shenzhen city are consistent with the double peak at this latitudinal. In this study, we found that the number of reported cases in Shenzhen city increased over time, and the increased incidence suggests that the HFMD pathogen is still spreading in the city. In addition, we found that the four districts with the lowest incidence rates (Nanshan District, Futian District, Luohu District, and Yantian District) also had the top four GDP per capita in Shenzhen city, while the two districts with higher incidence rates had relatively lower GDP per capita in the city's ranking (33). The economic level affects the local sanitation and, therefore, HFMD, as it is an intestinal infectious disease.

There are some limitations in this paper. In terms of model selection, in the past, we also chose the SEIAR model for simulation when cases were reported in a daily reported data (23, 37). However, it is rather unfortunate that exposed and asymptomatic populations were not introduced into the model because of the inability to obtain more precise units of incidence data due to the reporting of HFMD month-based data in Shenzhen. Meanwhile, the SIR model had stronger stability and avoided the problem of parameter identifiability than SEIAR (38). In addition, for the setting of susceptible individuals, we set all patients before 2010 to be susceptible. This is because HMFD was only incorporated into the infectious disease management system in China in 2008, and thus data on the recovered population are not available. The proportion of patients before that time is very small compared to the 10 million population in Shenzhen, so we believe that the effect is weak. The omission of this population may have an impact on the results of the age group component of the study, which needs to be discussed in further studies in future.

The results of the goodness-of-fit tests revealed that our model fitted the reported data well, which means that our modeling procedure has good validity and the model can be used to calculate the transmission rate of each pathogen.

In terms of the overall pattern, R_eff was >1 in Shenzhen city most of the time and in all administrative district, indicating that HFMD is still prevalent in the city and the disease burden remains serious, requiring attention and prevention and control concerns, whilst R_eff in Longgang district showed a gradual upward trend from 2010 to 2016. This may be due to the fact that Longgang, as a newly developed administrative region, has experienced rapid growth in population size and density in recent years, whereas the development of health care resources has failed to match the rate of population growth, leading to a rise in the transmissibility of HFMD. In contrast, the R_eff values in other districts experienced ups and downs between 1 and 11 with no significant trend.

When comparing the R_eff of the three age groups, it is clear that EV71 is most transmissible in the 0–2 years old group, while the R_eff of the 3–5 years old group is perennially lower than the other two groups. This may be due to the relatively high concentration of activities in the community among infants aged 0–2 years. Therefore, the proportion of infants and young children in the exposed population is also higher, leading to more disease transmission in this age group, and also to a greater R_eff of EV71 in the 3–5 years and the over 5 years group.

The 3–5 years age group should show a greater transmissibility due to increased socialization, while the smaller R_eff may be due to certain hygienic habits developed after school. Also, today's society shows concern for the disease and the parents and kindergarten teachers give more hygiene protection. The higher R_eff in the over 5 years old group compared to the previous group may be related to the older age group of EV71 in recent years, leading to an increased exposure rate. It can be seen that the transmission coefficient of HFMD for children in the 3–5 years group started to increase year by year in 2016. The greater contribution of transmission of the 3–5 years group to the transmission of the 0–3 years group may also be related to the comprehensive two-child policy launched in 2015.

After the launch of the HFMD vaccine in 2016, it is observed that the R_eff also declined in the city and in all districts except Yantian District. The results of the TSIR model developed by JT. Wu et al. similarly suggest that a mass EV-A71 vaccination program for infants and young children might significantly reduce the overall burden of HFMD (22). In our study, differences in the transmission capacity of different age groups for around this time point could be observed. However, the peak in the number of cases after the launch of the vaccine in the autumn of 2017 was significantly higher than in previous years, which may be due to the significant increase in the number of children after the introduction of the two-child policy, resulting in an increase in the number of HFMD cases, and may also account for the higher values of β₂₁ in the second-stage fitting results in 2017.

In summary, we estimate that the disease burden of HFMD in Shenzhen city is severe with varying disease trends across districts. There is variability in transmission in different age groups, all of which are slightly weaker, respectively, relative to the overall picture. In addition, vaccination efforts have been effective, but prevention and control of HFMD caused by enterovirus such as Cox A16 and others need to be further strengthened.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

PL, JR, TC, and YN made substantial contributions to conception and design. PL, QL, JR, TC, and YN collected the data and conceived the experiments. PL, QL, YW, FX, JR, ZL, and YN conducted the experiments and analyzed the results. JH, PL, JR, BD, LL, CL, WL, FX, TY, and YW involved in the visualization of the results. PL, JR, SY, CL, YN, and YW wrote the manuscript. QL, JR, PL, TC, and YN revised it critically for important intellectual content. All authors approved the final manuscript and agreed to be accountable for all aspects of the work.

Funding

This study was partly supported by the Bill & Melinda Gates Foundation (INV-005834).

Conflict of Interest

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

Publisher's Note

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

References

1. Hii YL, Rocklov J, Ng N. Short term effects of weather on hand, foot and mouth disease. PLoS One. (2011) 6:e16796. doi: 10.1371/journal.pone.0016796

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Chua KB, Kasri AR. Hand foot and mouth disease due to enterovirus 71 in Malaysia. Virol Sin. (2011) 26:221–8. doi: 10.1007/s12250-011-3195-8

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Fujimoto T, Iizuka S, Enomoto M, Abe K, Yamashita K, Hanaoka N, et al. Hand, foot, and mouth disease caused by coxsackievirus A6, Japan, (2011). Emerg Infect Dis. (2012) 18:337–9. doi: 10.3201/eid1802.111147

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ang LW, Koh BK, Chan KP, Chua LT, James L, Goh KT. Epidemiology and control of hand, foot and mouth disease in Singapore, 2001-2007. Ann Acad Med Singap. (2009) 38:106–12.

PubMed Abstract | Google Scholar

5. Hand Foot Mouth Disease Situation Update 2018. (2018). Available online at: https://www.who.int/westernpacific/health-topics/hand-foot-and-mouth-disease-(hfmd)

Google Scholar

6. Overview of the National Epidemic of Statutory Infectious Diseases in 2017. (2017). Available online at: http://www.nhc.gov.cn/jkj/s3578/201802/de926bdb046749abb7b0a8e23d929104.shtml

Google Scholar

7. Health GPDo. Prevalence of National Legal Infectious Diseases in Guangdong Province. (2017).

Google Scholar

8. Xiao X, Gasparrini A, Huang J, Liao Q, Liu F, Yin F, et al. The exposure-response relationship between temperature and childhood hand, foot and mouth disease: A multicity study from mainland China. Environ Int. (2017) 100:102–9. doi: 10.1016/j.envint.2016.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Li XW, Ni X, Qian SY, Wang Q, Jiang RM, Xu WB, et al. Chinese guidelines for the diagnosis and treatment of hand, foot and mouth disease (2018 edition). World J Pediatr. (2018) 14:437–47. doi: 10.1007/s12519-018-0189-8

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Abdennebi B, Bendisari K, Mansouri B, Abada M. Intracranial tumor developing from the vitelline sac. Apropos of a case] Neurochirurgie. (1988) 34:133–6.

PubMed Abstract | Google Scholar

11. Zou JJ, Jiang GF, Xie XX, Huang J, Yang XB. Application of a combined model with seasonal autoregressive integrated moving average and support vector regression in forecasting hand-foot-mouth disease incidence in Wuhan, China. Medicine (Baltimore). (2019) 98:e14195. doi: 10.1097/MD.0000000000014195

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ooi MH, Wong SC, Lewthwaite P, Cardosa MJ, Solomon T. Clinical features, diagnosis, and management of enterovirus 71. Lancet Neurol. (2010) 9:1097–105. doi: 10.1016/S1474-4422(10)70209-X

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Klein M, Chong P. Is a multivalent hand, foot, and mouth disease vaccine feasible? Hum Vaccin Immunother. (2015) 11:2688–704. doi: 10.1080/21645515.2015.1049780

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Xing W, Liao Q, Viboud C, Zhang J, Sun J, Wu JT, et al. Hand, foot, and mouth disease in China, 2008-12: an epidemiological study. Lancet Infect Dis. (2014) 14:308–18. doi: 10.1016/S1473-3099(13)70342-6

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Wu JT, Jit M, Zheng Y, Leung K, Xing W, Yang J, et al. Routine pediatric enterovirus 71 vaccination in China: a cost-effectiveness analysis. PLoS Med. (2016) 13:e1001975. doi: 10.1371/journal.pmed.1001975

PubMed Abstract | CrossRef Full Text | Google Scholar

16. He X, Dong S, Li L, et al. Using a Bayesian spatiotemporal model to identify the influencing factors and high-risk areas of hand, foot and mouth disease (HFMD) in Shenzhen city. PLoS Negl Trop Dis. (2020) 14:e0008085. doi: 10.1371/journal.pntd.0008085

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Yan XF, Gao S, Xia JF, Ye R, Yu H, Long JE. Epidemic characteristics of hand, foot, and mouth disease in Shanghai from 2009 to 2010: Enterovirus 71 subgenotype C4 as the primary causative agent and a high incidence of mixed infections with coxsackievirus A16. Scand J Infect Dis. (2012) 44:297–305. doi: 10.3109/00365548.2011.634433

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Chen Y, Badaruddin H, Lee VJ, Cutter J, Cook AR. The effect of school closure on hand, foot, and mouth disease transmission in singapore: a modeling approach. Am J Trop Med Hyg. (2018) 99:1625–32. doi: 10.4269/ajtmh.18-0099

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Song C, He Y, Bo Y, Wang J, Ren Z, Yang H. Risk assessment and mapping of hand, foot, and mouth disease at the county level in mainland china using spatiotemporal zero-inflated bayesian hierarchical models. Int J Environ Res Public Health. (2018) 15. doi: 10.3390/ijerph15071476

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yan S, Wei L, Duan Y, Li H, Liao Y, Lv Q, et al. Short-term effects of meteorological factors and air pollutants on hand, foot and mouth disease among children in Shenzhen city, China, 2009-2017. Int J Environ Res Public Health. (2019) 16. doi: 10.3390/ijerph16193639

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Qi H, Chen Y, Xu D, Su H, Zhan L, Xu Z, et al. Impact of meteorological factors on the incidence of childhood hand, foot, and mouth disease (HFMD) analyzed by DLNMs-based time series approach. Infect Dis Poverty. (2018) 7:7. doi: 10.1186/s40249-018-0388-5

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Takahashi S, Liao Q, Van Boeckel TP, et al. Hand, Foot, and Mouth Disease in China: Modeling Epidemic Dynamics of Enterovirus Serotypes and Implications for Vaccination. PLoS Med. (2016) 13:e1001958. doi: 10.1371/journal.pmed.1001958

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Luo K, Rui J, Hu S, et al. Interaction analysis on transmissibility of main pathogens of hand, foot, and mouth disease: a modeling study (a STROBE-compliant article). Medicine (Baltimore). (2020) 99:e19286. doi: 10.1097/MD.0000000000019286

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Lai CC, Jiang DS, Wu HM, Chen HH. A dynamic model for the outbreaks of hand, foot, and mouth disease in Taiwan. Epidemiol Infect. (2016) 144:1500–11. doi: 10.1017/S0950268815002630

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Wang J, Xiao Y, Cheke RA. Modelling the effects of contaminated environments on HFMD infections in mainland China. Biosystems. (2016) 140:1–7. doi: 10.1016/j.biosystems.2015.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Li Y, Huang M, Peng L. A multi-group model for estimating the transmission rate of hand, foot and mouth disease in mainland China. Math Biosci Eng. (2019) 16:2305–21. doi: 10.3934/mbe.2019115

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Koh WM, Badaruddin H, La H, Chen MI, Cook AR. Severity and burden of hand, foot and mouth disease in Asia: a modelling study. BMJ Glob Health. (2018) 3:e000442. doi: 10.1136/bmjgh-2017-000442

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bo YC, Song C, Wang JF Li XW. Using an autologistic regression model to identify spatial risk factors and spatial risk patterns of hand, foot and mouth disease (HFMD) in Mainland China. BMC Public Health. (2014) 14:358. doi: 10.1186/1471-2458-14-358

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang W, Du Z, Zhang D, Yu S, Huang Y, Hao Y. Assessing the impact of humidex on HFMD in Guangdong Province and its variability across social-economic status and age groups. Sci Rep. (2016) 6:18965. doi: 10.1038/srep18965

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhao J, Jiang F, Zhong L, Sun J, Ding J. Age patterns and transmission characteristics of hand, foot and mouth disease in China. BMC Infect Dis. (2016) 16:691. doi: 10.1186/s12879-016-2008-y

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Tian CW, Wang H, Luo XM. Time-series modelling and forecasting of hand, foot and mouth disease cases in China from 2008 to (2018). Epidemiol Infect. (2019) 147:e82. doi: 10.1017/S095026881800362X

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Municipality MBoS. Overview of Shenzhen City's Climate and Seasons. (2020).

Google Scholar

33. Bureau SS. Shenzhen City Statistical Yearbook. (2017).

Google Scholar

34. Yang Y, Sugimoto JD, Halloran ME, Basta NE, Chao DL, Matrajt L, et al. The transmissibility and control of pandemic influenza A (H1N1) virus. Science. (2009) 326:729–33. doi: 10.1126/science.1177373

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Zhuang ZC, Kou ZQ, Bai YJ, Cong X, Wang LH, Li C, et al. Epidemiological research on hand, foot, and mouth disease in mainland china. Viruses. (2015) 7:6400–11. doi: 10.3390/v7122947

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zhuang ZC, Kou ZQ, Bai YJ, Cong X, Wang LH, Li C, et al. The epidemiology of hand, foot and mouth disease in Asia: a systematic review and analysis. Pediatr Infect Dis J. (2016) 35:e285–300. doi: 10.1097/INF.0000000000001242

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Huang Z, Wang M, Qiu L, Wang N, Zhao Z, Rui J, et al. Seasonality of the transmissibility of hand, foot and mouth disease: a modelling study in Xiamen City, China. Epidemiol Infect. (2019) 147:e327. doi: 10.1017/S0950268819002139

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Tang L, Zhou Y, Wang L, Purkayastha S, Zhang L, He J, et al. A review of multi-compartment infectious disease models. Int Stat Rev. (2020) 88:462–513. doi: 10.1111/insr.12402

PubMed Abstract | CrossRef Full Text | Google Scholar