Front. Public Health Frontiers in Public Health Front. Public Health 2296-2565 Frontiers Media S.A. 10.3389/fpubh.2022.850369 Public Health Original Research Analysis of HFMD Transmissibility Among the Whole Population and Age Groups in a Large City of China Li Peihua 1 † Rui Jia 1 † Niu Yan 2 † Xie Fang 1 Wang Yifang 1 Li Zhuoyang 1 Liu Chan 1 Yu Shanshan 1 Huang Jiefeng 1 Luo Li 1 Deng Bin 1 Liu Weikang 1 Yang Tianlong 1 Li Qun 2 * Chen Tianmu 1 * 1State Key Laboratory of Molecular Vaccinology and Molecular Diagnostics, School of Public Health, Xiamen University, Xiamen, China 2Chinese Center for Disease Control and Prevention, Public Health Emergency Center, Beijing, China

Edited by: Ho Cheung William Li, The University of Hong Kong, Hong Kong SAR, China

Reviewed by: Weibing Wang, Fudan University, China; Daihai He, Hong Kong Polytechnic University, Hong Kong SAR, China

 *Correspondence: Tianmu Chen 13698665@qq.com; chentianmu@xmu.edu.cn Qun Li liqun@chinacdc.cn

This article was submitted to Children and Health, a section of the journal Frontiers in Public Health

†These authors have contributed equally to this work

 11 04 2022 2022 10 850369 07 01 2022 02 03 2022 Copyright © 2022 Li, Rui, Niu, Xie, Wang, Li, Liu, Yu, Huang, Luo, Deng, Liu, Yang, Li and Chen. 2022 Li, Rui, Niu, Xie, Wang, Li, Liu, Yu, Huang, Luo, Deng, Liu, Yang, Li and Chen

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

 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 Reff 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−7–2.31 ×10−4) while age group 3 had the lowest incidence (0–3.54 ×10−5). The differences in weekly incidence of EV71 between age groups were statistically significant (t12 = 7.563, P < 0.0001; t23 = 12.420, P < 0.0001; t13 = 16.996, P < 0.0001). The R2 of the SIR model Shenzhen city population-wide HFMD fit for each region was >0.5, and P < 0.001. Reff 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 Reff values between regions were judged by using analysis of variance (ANOVA) (F = 0.541, P = 0.744). SiIiRi-SjIjRj models between age groups had R2 over 0.7 for all age groups and P <0.001. The Reff 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.

 HFMD children mathematical model age group infectious diseases INV-005834 Bill and Melinda Gates Foundation10.13039/100000865 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 <italic>SIR</italic> 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).

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.

dSdt=brN-βSI-drSdIdt=βSI-γI-(dr+f)IdRdt=γI-drRN=S+I+R <italic>S</italic><sub><italic>i</italic></sub><italic>I</italic><sub><italic>i</italic></sub><italic>R</italic><sub><italic>i</italic></sub>-<italic>S</italic><sub><italic>j</italic></sub><italic>I</italic><sub><italic>j</italic></sub><italic>R</italic><sub><italic>j</italic></sub> Model

Considering the interactions between age groups, we created the age group SIR model (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.

dS1dt=brN-β11S1I1-β21S1I2-β31S1I3dI1dt=β11S1I1+β21S1I2+β31S1I3-γ′I1dR1dt=γ′I1 dS2dt=-β22S2I2-β12S2I1-β32S2I3dI2dt=β22S2I2+β12S2I1+β32S2I3-γ′I2dR2dt=γ′I2 dS3dt=-β33S3I3-β13S3I1-β23S3I2-drS3dI3dt=β33S3I3+β13S3I1+β23S3I2-γ′I3-drI3dR3dt=γ′I3-drR3 N=S1+S2+S3+I1+I2+I3+R1+R2+R3 Parameter Estimation

As shown in Table 1, there are eight parameters (β1, γ, 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 f1, f2, and f3 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. β11 denotes the transmission coefficient among individuals in age group 1; β22 denotes the transmission coefficient among individuals in age group 2 minus; β12 denotes the transmission coefficient from individuals in age group 1 to individuals in age group 2; β13 denotes the transmission coefficient from individuals in age group 1 to individuals in age group 2; and β14 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; β21 denotes the propagation coefficient from individuals within age group 2 to individuals in age group 1; β23 denotes the propagation coefficient from individuals within age group 2 to individuals in age group 3; β31 denotes the propagation coefficient from individuals within age group 3 to individuals in age group 1; and β32 denotes the propagation coefficient from individuals within age group 3 to individuals in age group 2. The Reff 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.

Parameter estimation.

Parameter Description Value Range Method
β Inter-individual propagation coefficient Seen text ≥0 Curve fitting
γ Relative recovery rate 1/14 - References
br Population birth rate 3.52 × 10−4 3.30 × 10−4-3.83 × 10−4 Analysis on the reported data
dr Population mortality rate 1.29 × 10−5 1.27 × 10−5-1.87 × 10−5 Analysis on the reported data
f Case fatality rate 0.0003 0–1 References
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 (Reff) 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 Reff 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.

ddt(I1I2I3)=(S1=(I1β11+I2β21+I3β31)S2=(I1β12+I2β22+I3β32)S3=(I1β13+I2β23+I3β33))−(I1γI2dr+I2γI3dr+I3γ)                      =defF−V

Step 2: Taking the derivatives of the vectors F and V with respect to the variable [I1 I2 I3], respectively, yields the corresponding Jacobi matrices F and V :

F=(S1β11S1β21S1β31S2β12S2β22S2β32S3β13S3β23S3β33‏)V=(γ000dr+γ000dr+γ‏)

and further calculate the V −1:

V-1=(1γ0001dr+γ0001dr+γ‏)

Step 3: Calculate the next generation matrix M = F V −1. The elements in row j and column i of the matrix M are the Reff, ij. of the intergroup interactions. The result of the calculation of Rij can be obtained by replacing Si=Ni into the equation.

M=(S1β11γS1β21dr+γS1β31dr+γS2β12γS2β22dr+γS2β32dr+γS3β13γS3β23dr+γS3β33dr+γ‏)

The elements in row j and column i of the matrix M are the Reff, ij. of the intergroup interactions. The result of the calculation of Rij can be obtained by replacing Si=Ni into the equation.

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

R0=λmax(M)

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 Rij and then the maximum eigenvalue of the numerical matrix M.

2) Use the definition method to calculate the overall R0.

R0=∑i=1n[P(x0∈Ni)∑i=1nRij]

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

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−2-7.53 × 10−4) and Futian having the lowest incidence rate (0–2.75 × 10−4). 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.

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

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

Comparison groups Baoan district Futian district Nanshan district Yantian district Luohu district Longgang district
Baoan district
Futian district 0
Nanshan district 0 0
Yantian district 0 0.0002 0
Luohu district 0 0.7209 0 0
Longgang district 0 0.0006 0 0.7404 0.0020

α = 0.05.

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 (t12 = 7.563, P < 0.0001; t23 =12.420, P < 0.0001; t13 =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.

Fitting results for EV71 incidence rates by age group.

 Model Fitting Effect <italic>SIR</italic> 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 R2 >0.5 and p < 0.001 for all districts (Figure 4). The quartiles of reported incidence in Shenzhen city were P25 = 1.78 × 10−5, P50 = 5.26 × 10−5, P75 = 1.19 × 10−4. The total incidence reported in Shenzhen city had an increasing trend year by year, with the highest peak in 2017 (4.51 × 10−4), while the remaining districts showed a trend of increasing and then decreasing trend, except for Futian and Nanshan districts.

 <italic>S</italic><sub><italic>i</italic></sub><italic>I</italic><sub><italic>i</italic></sub><italic>R</italic><sub><italic>i</italic></sub>-<italic>S</italic><sub><italic>j</italic></sub><italic>I</italic><sub><italic>j</italic></sub><italic>R</italic><sub><italic>j</italic></sub> Model

We also developed the SiIiRi-SjIjRj 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 R2 values for each group are listed in the listed fit effect table, and it can be seen that R2 is over 0.7 for all age groups, and p < 0.001 for all (Table 3). The quartiles for age group 1 were P25 = 1.08 × 10−5, P50 = 3.17 × 10−5, P75 = 6.58 × 10−5, age group 2 were P25 = 4.98 × 10−6, P50 = 1.23 × 10−5, P75 = 2.62 × 10−5, and age group 3 were P25 = 8.26 × 10−7, P50 = 2.26 × 10−6, P75 = 4.97 × 10−6. The transmission coefficients β for each age group at each time period were obtained by fitting, and the fitting results are shown in Table 4.

The P values and fitted R2 values for each group.

Age groups R2 P
0–3 years old group 0.818 <0.0001
3–5 years old group 0.715 <0.0001
>5 years old group 0.726 <0.0001

The weekly transmission coefficients β for each age group.

Year Weeks β11 β22 β33 β12 β13 β21 β23 β31 β32
2010 1–8 0.228 0.245 0.020 0.000 0.000 0.015 0.097 0.030 0.003
8–29 0.127 0.120 0.004 0.014 0.007 0.069 0.014 0.149 0.002
29–52 0.128 0.149 0.023 0.001 0.000 0.038 0.023 0.016 0.001
2011 1–5
5–33 0.018 0.096 0.009 0.015 0.009 0.190 0.000 0.045 0.005
34–52 0.260 0.221 0.129 0.009 0.023 0.157 0.011 0.052 0.436
2012 1–8
9–36 0.118 0.161 0.018 0.000 0.000 0.044 0.025 0.197 0.000
37–52 0.079 0.069 0.046 0.001 0.002 0.021 0.001 0.004 0.001
2013 1–5
6–12 0.226 0.194 0.024 0.021 0.022 0.182 0.001 0.091 0.075
13–40 0.076 0.035 0.214 0.108 0.000 0.138 0.000 0.046 0.085
41–52 0.427 0.294 0.225 0.142 0.031 0.251 0.027 0.083 0.456
2014 1–8
9–45 0.254 0.050 0.039 0.114 0.018 0.000 0.003 0.046 0.003
45–52 0.222 0.145 0.004 0.020 0.010 0.114 0.004 0.123 0.681
2015 1–10
11–48 0.001 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.001
48–53 0.127 0.008 0.007 0.099 0.006 0.248 0.015 0.006 0.456
2016 1–23
24–52 0.242 0.017 0.008 0.099 0.018 0.000 0.012 0.000 0.110
2017 1–27 0.043 0.055 0.001 0.050 0.020 0.151 0.005 0.006 0.076
27–53 0.000 0.018 0.001 0.026 0.004 1.036 0.016 0.002 0.007
Assessment of Transmissibility

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

The transimissibility of the first stage.

The median Reff of Shenzhen city was 4.001 (Interquartile range [IQR]: 0.021–8.929), with Reff values <1 in the second cycle of 2011; the median Reff for Baoan was 4.069 (IQR: 2.096–10.139), with Reff values all over 1; the median Reff of Longgang District was 4.126 (IQR: 2.935–8.922), with Reff values all over 1 and a gradual upward trend from 2010 to 2016; the median Reff of Luohu District had a median Reff of 4.298 (IQR: 1.730- 9.314), with Reff values over 1; the median Reff of Nanshan District was 3.939 (IQR: 1.937- 8.002), with Reff values all over 1; the median Reff of Futian District was 4.777 (IQR: 2.924- 11.105), with a decreasing trend of Reff over 1; the median Reff in Yantian District was 5.131 (IQR: 1.730- 9.314), with Reff over 1. Reff 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 Reff 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 Reff 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 Reff for the 3–5 years age group fluctuated within a certain range and had the largest value in the second half of 2017. Reff values for the over 5 years age showed an increasing and then decreasing trend, peaking at the end of 2013 (Figure 7A).

The transimissibility of the second stage (A) Group1 indicates the Reff value of 0–2years old age group; Group 2 indicates the Reff value of 3–5 years old age group; Group 3 indicates the Reff 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 (Reff: M = 1.787, IQR = 0–9.146), followed by that of Group 1 to Group 2 (Reff12: M = 0.287, IQR = 0–1.988) and finally Group 1 to Group 3 (Reff: M = 0.287, IQR = 0–1.988 median: 0.287, range: 0–1.988). The Reff 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 Reff 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 Reff year by year.

In the over 5 years age group, the Reff 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 Reff 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, Reff 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 Reff 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 Reff values in other districts experienced ups and downs between 1 and 11 with no significant trend.

When comparing the Reff of the three age groups, it is clear that EV71 is most transmissible in the 0–2 years old group, while the Reff 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 Reff 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 Reff 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 Reff 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 Reff 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 β21 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 Hii YL Rocklov J Ng N. Short term effects of weather on hand, foot and mouth disease. PLoS One. (2011) 6:e16796. 10.1371/journal.pone.001679622993785 Chua KB Kasri AR. Hand foot and mouth disease due to enterovirus 71 in Malaysia. Virol Sin. (2011) 26:221–8. 10.1007/s12250-011-3195-821847753 Fujimoto T Iizuka S Enomoto M Abe K Yamashita K Hanaoka N . Hand, foot, and mouth disease caused by coxsackievirus A6, Japan, (2011). Emerg Infect Dis. (2012) 18:337–9. 10.3201/eid1802.11114722304983 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.19271036 Hand Foot Mouth Disease Situation Update 2018. (2018). Available online at: https://www.who.int/westernpacific/health-topics/hand-foot-and-mouth-disease-(hfmd) 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 Health GPDo. Prevalence of National Legal Infectious Diseases in Guangdong Province. (2017). Xiao X Gasparrini A Huang J Liao Q Liu F Yin F . 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. 10.1016/j.envint.2016.11.02128069250 Li XW Ni X Qian SY Wang Q Jiang RM Xu WB . Chinese guidelines for the diagnosis and treatment of hand, foot and mouth disease (2018 edition). World J Pediatr. (2018) 14:437–47. 10.1007/s12519-018-0189-830280313 Abdennebi B Bendisari K Mansouri B Abada M. Intracranial tumor developing from the vitelline sac. Apropos of a case] Neurochirurgie. (1988) 34:133–6.2457175 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. 10.1097/MD.000000000001419530732135 Ooi MH Wong SC Lewthwaite P Cardosa MJ Solomon T. Clinical features, diagnosis, and management of enterovirus 71. Lancet Neurol. (2010) 9:1097–105. 10.1016/S1474-4422(10)70209-X20965438 Klein M Chong P. Is a multivalent hand, foot, and mouth disease vaccine feasible? Hum Vaccin Immunother. (2015) 11:2688–704. 10.1080/21645515.2015.104978026009802 Xing W Liao Q Viboud C Zhang J Sun J Wu JT . Hand, foot, and mouth disease in China, 2008-12: an epidemiological study. Lancet Infect Dis. (2014) 14:308–18. 10.1016/S1473-3099(13)70342-624485991 Wu JT Jit M Zheng Y Leung K Xing W Yang J . Routine pediatric enterovirus 71 vaccination in China: a cost-effectiveness analysis. PLoS Med. (2016) 13:e1001975. 10.1371/journal.pmed.100197527050103 He X Dong S Li L . 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. 10.1371/journal.pntd.000808532196496 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. 10.3109/00365548.2011.63443322176514 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. 10.4269/ajtmh.18-009930350767 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. 10.3390/ijerph1507147630002344 Yan S Wei L Duan Y Li H Liao Y Lv Q . 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. 10.3390/ijerph1619363931569796 Qi H Chen Y Xu D Su H Zhan L Xu Z . 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. 10.1186/s40249-018-0388-529391070 Takahashi S Liao Q Van Boeckel TP . Hand, Foot, and Mouth Disease in China: Modeling Epidemic Dynamics of Enterovirus Serotypes and Implications for Vaccination. PLoS Med. (2016) 13:e1001958. 10.1371/journal.pmed.100195826882540 Luo K Rui J Hu S . 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. 10.1097/MD.000000000001928632176053 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. 10.1017/S095026881500263026567705 Wang J Xiao Y Cheke RA. Modelling the effects of contaminated environments on HFMD infections in mainland China. Biosystems. (2016) 140:1–7. 10.1016/j.biosystems.2015.12.00126738934 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. 10.3934/mbe.201911531137214 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. 10.1136/bmjgh-2017-00044229564154 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. 10.1186/1471-2458-14-35824731248 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. 10.1038/srep1896526743684 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. 10.1186/s12879-016-2008-y27871252 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. 10.1017/S095026881800362X30868999 Municipality MBoS. Overview of Shenzhen City's Climate and Seasons. (2020). Bureau SS. Shenzhen City Statistical Yearbook. (2017). Yang Y Sugimoto JD Halloran ME Basta NE Chao DL Matrajt L . The transmissibility and control of pandemic influenza A (H1N1) virus. Science. (2009) 326:729–33. 10.1126/science.117737319745114 Zhuang ZC Kou ZQ Bai YJ Cong X Wang LH Li C . Epidemiological research on hand, foot, and mouth disease in mainland china. Viruses. (2015) 7:6400–11. 10.3390/v712294726690202 Zhuang ZC Kou ZQ Bai YJ Cong X Wang LH Li C . The epidemiology of hand, foot and mouth disease in Asia: a systematic review and analysis. Pediatr Infect Dis J. (2016) 35:e285–300. 10.1097/INF.000000000000124227273688 Huang Z Wang M Qiu L Wang N Zhao Z Rui J . Seasonality of the transmissibility of hand, foot and mouth disease: a modelling study in Xiamen City, China. Epidemiol Infect. (2019) 147:e327. 10.1017/S095026881900213932115010 Tang L Zhou Y Wang L Purkayastha S Zhang L He J . A review of multi-compartment infectious disease models. Int Stat Rev. (2020) 88:462–513. 10.1111/insr.1240232834402