Hand, Foot, and Mouth Disease in China: Modeling Epidemic Dynamics of Enterovirus Serotypes and Implications for Vaccination

Summary

Using time series susceptible–infected–recovered (TSIR) epidemic models, Saki Takahashi and colleagues investigate hand, foot, and mouth disease dynamics and the projected effect of vaccination in China.

Introduction

Hand, foot, and mouth disease (HFMD) is a common childhood illness caused by serotypes of the Enterovirus A species in the genus Enterovirus of the Picornaviridae family [1,2]. HFMD predominantly affects children younger than 5 y of age [3], and most patients exhibit a self-limiting illness that typically includes fever, skin eruptions on the hands and feet, and vesicles in the mouth. However, a small proportion of infections lead to the development of neurological and systemic complications that can be fatal, particularly in cases associated with the Enterovirus A71 (EV-A71) serotype [4]. Since 1997, EV-A71-associated HFMD epidemics have been increasingly reported across the Asia-Pacific region, chronologically in Malaysia [5], Taiwan of China [6], Japan [7], Singapore [8], Viet Nam [9], Mainland China [10,11], Hong Kong Special Administrative Region of China [12], South Korea [13], and Cambodia [14]. China reported 9 million cases of HFMD between 2008 and 2013, with the two serotypes EV-A71 and Coxsackievirus A16 (CV-A16) being responsible for 73% of these cases [15]. As such, HFMD is a growing public health concern that poses a considerable disease burden and economic impact in affected areas [16,17]. Recent concerns with disease due to outbreaks of Enterovirus D68 in the US and other countries have only increased the urgency to understand the local viral community dynamics of this group of infections [18].

Transmission of enteroviruses occurs by direct contact with the mucus, saliva, or feces of an infected individual, or through indirect contact via contaminated surfaces [19]. HFMD infection characteristically causes acute illness with a duration of approximately 1 wk [20–22]. There are no established antiviral treatments for HFMD, but three recent phase 3 clinical trials of inactivated monovalent EV-A71 vaccines manufactured in China were found to have high efficacy (90.0%–97.4%) against EV-A71-associated HFMD in infants and young children [23], and two of these vaccines were licensed in China in December 2015 [24]. However, these vaccines did not offer protection against HFMD caused by the CV-A16 serotype [25].

There is empirical evidence that infection with one serotype of a multi-serotype viral pathogen such as influenza or dengue virus can lead to transient immunity against other serotypes of the same infection [26–28]. There is also some historical and more recent evidence of interference between polioviruses (Enterovirus C species) and other enteroviruses [29–31]. For HFMD specifically, co-infection of a single individual with both serotypes is rarer than expected by chance [32,33], suggesting the existence of at least short-term cross-protection, and neutralization assays have shown partial cross-reactivity between the EV-A71 and CV-A16 serotypes [34,35], indicating a potential competitive interaction. However, the effects of multi-serotype interactions remain poorly understood to date. Mathematical models can be used to bridge this gap, allowing for the estimation of key parameters governing cross-protection and the evaluation of the epidemiological impact of vaccination. Previous modeling studies have used variations of the continuous-time susceptible–infected–recovered (SIR) model to study the seasonal patterns of HFMD and its basic and effective reproductive numbers [36–39]. These analyses have focused on estimating epidemiological parameters of either EV-A71-associated HFMD alone or all-cause HFMD, and have not accounted for the potential cross-protective effects of infection with a particular serotype on the dynamics of other serotypes.

To our knowledge, two important aspects of HFMD dynamics are yet to be addressed: first, the evaluation of the duration and strength of cross-protection between different serotypes and, second, potential serotype replacement by CV-A16 or other serotypes (e.g., CV-A6 or CV-A10) following vaccination against EV-A71, due to decreased EV-A71 incidence and an accompanying reduction in cross-immunity against CV-A16 and other serotypes, as well as potential increased circulation of these viruses in the population. Such serotype replacement has been observed, for example, for pneumococcal disease, with increases in the prevalence of non-vaccine serotypes following use of the pneumococcal heptavalent conjugate vaccine [40–42]. Vaccination against EV-A71 may potentially have important effects on the burden of disease caused by CV-A16 and other serotypes, and robust predictions of these effects are necessary before introduction of an EV-A71-containing vaccine into a population.

In this analysis, we adopted a mathematical modeling approach to study the multi-serotype transmission dynamics of HFMD in the Asia-Pacific region. We linked models closely to surveillance time series, using a multi-serotype time series epidemic model to examine the dynamics of HFMD caused by EV-A71 and CV-A16 in the 31 provinces of Mainland China between 2009 and 2013. We evaluated the potential existence of cross-protection between the two serotypes and the risk of serotype replacement by CV-A16 following prospective EV-A71 vaccination. We also extended the model to assess the co-circulation of additional non-EV-A71, non-CV-A16 enteroviral serotypes.

Methods

Data

We used three sources of data in this analysis. First, weekly reports of HFMD incidence between 1 January 2009 and 31 December 2013 (a total of 262 wk) were obtained from a national surveillance system maintained by the Chinese Center for Disease Control and Prevention (China CDC) in Beijing, China. These reports were available at two spatial scales: the seven administrative regions and 31 provinces of China.

Second, we obtained time series of weekly laboratory-confirmed HFMD for a sub-sample of cases from each province, with virological test results classified as EV-A71, CV-A16, other non-EV-A71 and non-CV-A16 serotypes of enterovirus, or negative [15]. This information was aggregated to the regional scale on a monthly basis to reduce potential biases from the sampling scheme (S1 Fig). The proportions of each serotype were estimated, with all positive and negative tests included in the denominator. If there were no virological test data available for a given month, the proportions from the virological tests of the previous month were substituted. We applied these proportions to the reported case counts from the surveillance registry to estimate serotype-specific weekly incidence by province (Fig 1; S2 Dataset). Since infection with EV-A71 or CV-A16 accounted for the majority of total laboratory-confirmed HFMD cases between 2008 and 2013 (73%), we limited the scope of the main analysis to infection with either of these two serotypes. We used a time step of 1 wk: the incubation period of HFMD is 7 to 14 d and viral excretion persists for about 2 wk after symptom onset, so the generation time would be contained in this time frame. It has previously been shown that the estimation of seasonality is robust to the length of the chosen time step [43].

Third, we obtained yearly birth rates and population sizes between 2009 and 2013 by region and by province from the National Bureau of Statistics of China (S1 Dataset) [44]. Although reports of HFMD cases were available from 2008, the time frame of this analysis was limited to the period between 1 January 2009 and 31 December 2013 because of the sparsity of laboratory-confirmed cases during the first year of surveillance.

Host and Transmission Dynamics

The time series susceptible–infected–recovered (TSIR) model is a discrete-time version of the continuous-time SIR model in which individuals are born and enter the susceptible class of individuals, become infected and infectious with a disease, and recover and are removed thereafter [45]. Here, we used a multi-serotype extension of the TSIR model that allows for a cross-protective effect between serotypes [46]. The susceptible compartment of the TSIR model is defined by:

where at each time step t, S_t,i is the number of individuals that are susceptible to serotype i, B_t is the number of births (known from demographic data), It,i(r) is the reported number of individuals infected with serotype i, ρ_i is the reporting rate of serotype i (assumed to be constant over the entire time period), and CP_t,i represents the effect of cross-protection following infection with serotype i against all other serotypes. We assumed that there is no maternal immunity period for HFMD. Assuming that all individuals are born susceptible to HFMD and eventually become infected, we rearranged Eq 1 as follows:

where S¯i is the mean proportion of individuals that are susceptible to serotype i, N¯ is the mean number of individuals in the population over the time period of interest, and D_0,i is the unknown deviation around the mean number of individuals in the population that are susceptible to serotype i at the beginning of the time series. We reconstructed the time series of susceptible individuals by first extracting D_t,i, the residuals from the following linear regression:

The ρ_i values were estimated as the slope of this regression (using the ρ_{j ≠ i} parameters as iteratively estimated offset terms as per [46], since the I_t,i and CP_t,i terms both depend on ρ_i) and are quite low because of the mild nature of the infection [47], though the reporting rates of EV-A71 are generally higher than those of CV-A16 due to the severity of symptoms. The D_t,i values were then added to the mean number of susceptible individuals in the population (S¯i⋅N¯, with S¯i estimated using marginal profile likelihoods and its 95% confidence interval [CI] derived from the χ² distribution with 1 degree of freedom) to yield the complete time series of susceptible individuals:

In the multi-serotype TSIR model framework, HFMD transmission is characterized by the following stochastic frequency-dependent dynamics from the Poisson distribution:

where I_t,i is the number of infected individuals adjusted for reporting rate, λ_t+1,i is the expected number of individuals infected with serotype i at time t + 1, β_s,i is a seasonally varying transmission rate, α is a correction parameter accounting for nonseasonal heterogeneities in mixing [48,49] as well as for time discretization [50], and N_t is the total population size at time t. We linearized Eq 7 and estimated β_s,i with the following regression model:

Extensive simulation studies indicate that the log-linear model is a better regression model than the log-transformed linear regression model (S31 and S32 Figs). Additionally, the log-linear model is a tractable approach to approximating the negative binomial distribution, which is formally a more satisfying representation of the underlying Reed–Frost epidemic model [51]. Because weekly case counts are relatively high at the province level (and define the over-dispersion of the negative binomial distribution), we can closely reflect these dynamics using the log-linear model and avoid the need to explicitly account for over-dispersion.

Predictions for S_t+1,i and I_t+1,i were generated using Eqs 1 and 6, respectively, performing separate estimation on each of the 31 provinces of China and setting I_t+1,i = λ_t+1,i in the deterministic simulations. We used α = 0.95 throughout this analysis for consistency in comparing seasonal patterns across provinces; a sensitivity analysis varying α values between 0.91 and 0.99 is provided in S22–S30 Figs. Co-infection was omitted from the two-serotype model due to its relatively rare occurrence in practice (S4 Table). Model fit was assessed by comparing observed incidence against simulated incidence, and calculating the coefficient of determination (R²) and absolute yearly prediction error (PE):

Cross-Protection, Epidemiological Parameters, and Vaccine Simulation

The multi-serotype TSIR model allows for a cross-protective effect between serotypes. Supposing that infection with serotype i is fully immunizing against that serotype (e.g., exhibits SIR epidemic dynamics rather than susceptible–infected–recovered–susceptible [SIRS] epidemic dynamics) and that a proportion of those infected individuals also gain transient immunity against all serotypes j ≠ i, a fixed duration of cross-protection can be modeled as in [46]:

where δ is the proportion of those infected with one serotype (adjusted for reporting rate) who are cross-protected against infection from all other serotypes for a fixed duration of k time steps, and CP_t,i represents the effect of cross-protection following infection with serotype i. The product k ⋅ δ represents the average duration of cross-protection, incorporating individuals who do and do not experience cross-protection [46]. After the specified length of cross-protection, we assumed that individuals who had experienced cross-protection subsequently become re-susceptible to all other serotypes. 95% CIs for k and δ were derived from the profile likelihood using the χ² distribution with 1 degree of freedom.

We estimated two time-varying epidemiological parameters for each serotype: R₀ (basic reproductive number) and R_E (effective reproductive number), obtained by rearranging Eq 7 and defined as:

R₀ is the average number of secondary infectious persons resulting from one infectious person following his/her introduction into a completely susceptible population. R_E is the average number of secondary infectious persons resulting from one infectious person into a population that is partially immune; it is above 1 when incidence is increasing and below 1 when incidence is decreasing.

There is an inherent trade-off between having sufficient data to estimate the cross-protection parameters and having sufficient data to estimate the epidemiological parameters for the multi-serotype model. We chose to set aside the first year of the full laboratory-confirmed time series (2009) to back-fit the cross-protection parameters k and δ, and the remaining 4 y of data (2010–2013) were used for fitting the epidemiological parameters ρ_i and β_s,i. Thus, our estimate of k was necessarily restricted to be less than 1 y long, and we parsimoniously assumed the values of k and δ to be the same for infection with any serotype. Log-likelihood surfaces of the cross-protection parameters were generated by the following modified log-linear regression model, including the reporting rate as an additional offset term:

While our primary model is the two-serotype TSIR model for EV-A71 and CV-A16, we also explored a three-serotype model including the potential epidemic series encompassing non-EV-A71, non-CV-A16 serotypes of enterovirus (S6–S8 Tables).

To simulate the effects of introducing monovalent EV-A71 vaccination in the population, we assumed (based on findings from recent clinical trials) that vaccine-induced immunity against EV-A71 is narrow and has no immunological effect against CV-A16 infection [52]. In line with the high efficacy levels found in clinical trials, we also assumed a perfect vaccine efficacy of 100%. To capture the temporal dynamics in the simplest way, seasonality in transmission was omitted from vaccine simulations [53]. To model the effects of vaccination (assumed to be administered near birth), the size of each weekly birth cohort was reduced by the vaccine coverage. For example, if the achieved vaccine coverage was 90%, then 10% of the children born in a week would be susceptible to HFMD. We also assessed the impact of a broad monovalent EV-A71 vaccine that confers the same duration and strength of cross-protection against CV-A16 as natural infection with EV-A71, as well as the impact of a narrow bivalent EV-A71, CV-A16 vaccine that is equally protective against the two serotypes. All analysis was conducted using R version 3.0.2 (http://cran.r-project.org).

Ethical Approval

In May 2008, HFMD was added to the list of notifiable diseases in China. According to China’s law on the prevention and treatment of infectious diseases, personal identifiers should be collected for individual cases with diagnosis of a notifiable disease, for the purposes of public health surveillance and response. The National Health and Family Planning Commission of China decided that the collection of individual data for all notifiable diseases, including HFMD, according to the national surveillance protocol was part of an ongoing public health response and was thus exempt from institutional review board assessment.

The China CDC has strict regulations on how to protect patients’ privacy. The National Center for Public Health Surveillance and Information Services at the China CDC is responsible for the management of all disease surveillance data, and it anonymized the individual HFMD data by deleting the personal identifiers (such as patient name, parent name, home address, and telephone number) before the China CDC co-authors of this article, in the Division of Infectious Disease, were given access to the surveillance data for the purposes of research. The co-authors of this article did not participate in de-identifying the data and do not have the personal identifiers of the HFMD cases. Data were also aggregated by week, enterovirus serotype, and province by the China CDC co-authors before analysis.

Results

Model Fit

The two-serotype TSIR model provides a good fit to patterns of HFMD incidence in the 31 Chinese provinces. The median proportion of individuals susceptible (S¯) to EV-A71 was 0.067 (interquartile range [IQR]: 0.051, 0.153), and the median S¯ for CV-A16 was 0.061 (IQR: 0.051, 0.070); the mean reporting rate of both serotypes (ρ) was approximately 3% (S1 Table). Our national estimate of R₀, weighted by provincial population size, was 26.63 for EV-A71 (IQR: 23.14, 30.40) and 27.13 for CV-A16 (IQR: 23.15, 31.34), with considerable variation between provinces (Table 1). In comparison with separate one-serotype models (e.g., with no cross-protection) (S5 Table), the two-serotype model fit as well or better in terms of the mean absolute yearly PE in 20 of the 31 provinces for EV-A71, and 23 of the 31 provinces for CV-A16 (S7 Table).

In 27 of the 31 provinces, the maximum likelihood estimates from the two-serotype model indicate the existence of a transient cross-protection between the EV-A71 and CV-A16 serotypes (Table 2). The population-weighted national average suggests that infection with one serotype yields k = 9.95 wk (95% CI: 3.31, 23.40) of protection against infection with the other serotype in δ = 0.68 of the population (95% CI: 0.37, 0.96), resulting in an average duration of cross-protection of 6.77 wk (95% CI: 2.50, 10.03). We estimated likelihood surfaces of the cross-protection parameters by province to obtain 95% CIs (S15–S21 Figs).

The two-serotype model performed well in predicting weekly incidence of both EV-A71 and CV-A16 in the representative province of Beijing, capturing both the shape and scale of incidence from 2010 to 2013 (Fig 2). A sensitivity analysis varying α values in the two-serotype model for Beijing suggests that 0.95 is a feasible estimate of α (S22–S30 Figs) and that it is able to capture the yearly epidemic cycles while maintaining spatial stochasticity [54]; the actual value of α has quantitative implications for the impact of vaccination, though the qualitative conclusions given below are robust. As α approaches the value of 1, the model behavior becomes erratic because the underlying Reed–Frost epidemic model is neutrally stable and the TSIR approximation breaks down [50]. Two-serotype model fits of the observed and simulated time series for the remaining provinces are provided in S8–S14 Figs.

The time series were too short to do out-of-sample cross-validation, but the robustness of our inference procedure was tested by stochastically simulating incidence data, fitting the model, and recovering parameter values. The log-linear regression model recovered more reliable estimates of S¯ and β_s than the ordinary least squares (OLS) regression model with log-transformed incidence (S31 and S32 Figs). The estimated S¯ in the OLS regression is biased high, and the likelihood surface is flattened; the generalized linear models (GLMs) can recover reliable estimates of S¯, and the 95% profile likelihood intervals are more reasonable than those estimated from the OLS regression. The log-linear model was also able to robustly recover parameter values from simulated data incorporating various levels of cross-protection (S33 Fig).

Geographic Variations in Transmission

There are two distinct geographic trends of HFMD in China: in the north, there is an annual peak in incidence in June, and in the south, there are bi-annual peaks in May and in September–October [15]. The models captured these variations and accurately predicted the annual peaks in incidence in the northern provinces (e.g., Beijing [S11A Fig] and Tianjin [S11B Fig]) and the bi-annual peaks in the southern provinces (e.g., Chongqing [S13A Fig] and Sichuan [S13B Fig]) for EV-A71 and CV-A16. We also explored the seasonal patterns of R_E along a latitudinal gradient, which suggest that in the north both serotypes tend to have one period each year when incidence is increasing and R_E is above 1, while in the south there tend to be two distinct periods within a year during which incidence is increasing in the population (S3 Fig).

Based on the two-serotype model, we found considerable seasonal and geographic heterogeneities in the estimated values of R₀ (S2 and S3 Tables). We also disentangled the different sources of spatial and temporal variation in the transmission rate: within-week variation by province was estimated by the standard errors of the log(β_s) terms from the log-linear regression in Eq 8 and was found to be generally low, between-week variation by province was assessed by the coefficient of variation in β^s, and between-province variation was found to be generally high for both EV-A71 and CV-A16 (S1 Table; S2 Fig).

Impact of EV-A71 Vaccination

Using our estimate of R₀ of approximately 25 and the basic equation for calculating the critical vaccination threshold p_c = 1 − 1/R₀ [55], we estimate that a majority of provinces would require vaccine coverage levels upwards of 96% of infants and young children to end ongoing transmission of EV-A71-associated HFMD. We simulated the effects of nationally introducing three types of EV-A71-containing vaccine on future EV-A71 and CV-A16 incidence (assuming R₀ equals 25 for both serotypes and 90% vaccine coverage at birth as a base case): a broad monovalent EV-A71 vaccine (k_infection = k_vacc = 22 wk and δ_infection = δ_vacc = 1, the highest province-specific estimates of cross-protection from natural infection), a narrow monovalent EV-A71 vaccine conferring no cross-protection, and a narrow bivalent EV-A71, CV-A16 vaccine.

For all three vaccines, vaccination provides both a short- and long-term reduction in national EV-A71 incidence. While EV-A71 vaccination would not affect overall CV-A16 incidence in the long term, we found that a broad monovalent EV-A71 vaccine could cause a transient benefit in terms of reduced CV-A16 incidence due to cross-protective immunity (Fig 3A). The transient increase in CV-A16 in the second year following vaccine deployment reflects the end of a phenomenon known as a “honeymoon period,” during which susceptible individuals have slowly accumulated.

However, we found from our simulations that a narrow monovalent EV-A71 vaccine could cause a transient increase in national incidence and lead to a potential competitive release of CV-A16 (Fig 3B). This increase is due to the reduction in EV-A71 incidence that reduces the population levels of cross-immunity against CV-A16. The magnitude and duration of this increase following vaccination will depend on various factors including the achieved EV-A71 vaccine coverage (Fig 3D). Lastly, the deployment of a bivalent vaccine would lead to reductions in both EV-A71- and CV-A16-associated HFMD (Fig 3C).

Across individual provinces, we found that EV-A71 incidence decreases monotonically with higher vaccination coverage, and ongoing EV-A71 transmission approaches zero in the range of the critical vaccination threshold (S5A Fig). We also found that post-vaccination CV-A16 incidence increases slightly or remains comparable to its pre-vaccination levels with use of a narrow monovalent vaccine, and that the magnitude of this transient increase in CV-A16 incidence scales with EV-A71 vaccine coverage (S5B Fig). Further investigation of other combinations of cross-protection parameters led to the same qualitative conclusions (S6 and S7 Figs.).

Additional Serotype Circulation

We extended the two-serotype model to address the co-circulation of non-EV-A71, non-CV-A16 serotypes of enterovirus as an additional third epidemic series, though it is likely to be a combination of two or more enteroviruses including CV-A6 and CV-A10 [56–58]. In comparison with the two-serotype model, the three-serotype model (S6 Table) fit as well or better in terms of the mean absolute yearly PE in 20 of the 31 provinces for EV-A71, and 22 of the 31 provinces for CV-A16 (S7 Table). We did not see qualitative differences in our results with the inclusion of this additional serotype epidemic series. Multi-serotype models generate more multi-collinearity among parameters and introduce statistical identifiability issues, since a number of parameter value combinations are similarly supported by the likelihood. There is an inverse relationship between S¯ and β¯ (and therefore incidence) that is mediated by cross-protection (S4 and S34 Figs). Therefore, we parsimoniously used province-level cross-protection estimates from the two-serotype model to parametrize the three-serotype model.

Discussion

Findings

The two-serotype TSIR model captured the observed dual epidemic cycles of EV-A71- and CV-A16-associated HFMD for the 31 provinces in China and the varying seasonal patterns of transmission in the northern and southern provinces of the country. Comparing values of the mean absolute yearly PE, we found that the two-serotype model with cross-protection fit the incidence data better on average than the separate one-serotype models; the public health message of HFMD as a highly transmissible infection remained consistent between models. Additionally, estimates of the cross-protection parameters across provinces suggest that HFMD serotypes provide a temporary immunizing effect against each other. However, our results suggest considerable geographic variability in this effect, pointing to the substantial value of further study. Our estimated national R₀ of approximately 25 is consistent with the relatively low median age at HFMD infection [15], though any loss of immunity (see below) could also keep the age at primary infection low.

The simulations of EV-A71 vaccination showed that incidence would decrease monotonically with higher vaccine coverage, and suggest that a mass EV-A71 vaccination program targeted at infants and young children in China could greatly reduce the burden of HFMD caused by EV-A71. The high R₀ values of HFMD serotypes imply that disease elimination is unlikely and translate to necessary coverage levels of above 96% to achieve population-level immunity. Additionally, by assuming 100% vaccine efficacy, this estimate of required vaccination coverage is yet an underestimation. A considerable increase in the prevalence of non-vaccine serotypes has been observed in other pathogens such as pneumococcal disease [40–42], but our two-serotype model results suggest that the incidence of CV-A16 following narrow monovalent EV-A71 vaccination would be expected to remain comparable to its pre-vaccination levels. However, a bivalent or polyvalent vaccine would be desirable in light of the co-circulation of several enteroviruses associated with HFMD [59,60] and because the extent of indirect protection afforded by vaccination would not be major at these relatively high levels of R₀.

This analysis is novel in that it evaluates the existence and magnitude of cross-protection between different HFMD serotypes in a mathematical modeling context. Our framework allowed us to look for potential circulating serotype replacement by CV-A16 following vaccination against EV-A71. However, further studies are required to better understand the strength and breadth of immunity conferred by the EV-A71 vaccine relative to the levels of immunity following natural infection and to test for potential loss of immunity over time. While we acknowledge that 5 y may be too short a time series to definitively understand serotype interactions, based on our findings we are cautiously optimistic that monovalent EV-A71 vaccination is sufficient and will not significantly increase the risk of serotype replacement by CV-A16 in the population.

Caveats and Directions for Future Work

Using this relatively simple mathematical model, we were able to detect a robust signature of herd immunity driving the outbreak dynamics of HFMD. Although unlikely to dominate the observed incidence dynamics, significant evidence for antigenic drift in EV-A71 indicates some potential loss of immunity, which could affect subsequent immune escape, particularly in the context of vaccination [61,62]. In general, protection (via a rise in neutralizing antibody titers) against a wide spectrum of EV-A71 sub-genotypes is observed shortly after vaccination regardless of whether children already had antibodies [63–65], but longer follow-up is required to assess whether immunity against the sub-genotypes other than the vaccine serotype wanes (more rapidly) over time, as seen in influenza and dengue virus. The current time series were too short to evaluate the extent of loss of immunity following infection (captured, for example, by SIRS epidemic models [66,67]); however, the quality of the fit underlines that our assumption of strong immunity is reasonable as a first step.

The scope of the data also did not allow us to explicitly disentangle infection with alternative serotypes such as CV-A6 and CV-A10, which are known to be currently co-circulating in China and throughout the Asia-Pacific region [68–70] and might show further interactions. The increasingly multi-collinear nature of cross-protection parameters when additional serotypes are included in the model poses an analytical challenge. However, if cross-protection were substantially stronger than our estimates, we might expect to see the epidemics be out of phase. Understanding the interactions between co-circulating serotypes would allow for more nuanced estimates of the risk of serotype replacement by non-vaccine serotypes. In this analysis, we assumed stationarity in the mean transmission rate between years because of the short time series, but there may be confounding effects of year-to-year variation in transmissibility on estimates of cross-immunity. However, our results are conservative in terms of assessing potential serotype replacement because such exogenous forcing would increase correlations in the incidence of different serotypes, which would in turn reduce apparent cross-protection.

Additionally, the absence of age structure in this formulation of the TSIR model did not allow us to assess the degree of loss of homologous and heterologous immunity over time or the age-focusing of vaccination. Subsequent models could be further refined to allow for heterogeneity in mixing and transmission by age, and to model different vaccination implementation strategies. Furthermore, serological surveys would be useful to better understand population-level susceptibility to HFMD.

The richness of the spatial scale of these data allow for future work to better understand spatial correlations in incidence and synchrony in the timing of epidemics [71]. These data will also be useful for evaluating the spatial scale (national, regional, provincial, etc.) at which vaccination efforts should be coordinated and deployed. A better understanding of the dynamics of serotype invasion and interaction is necessary for accurately predicting the effects of introducing a new HFMD vaccine formulation with another serotype of enterovirus in China.

Since the recent outbreaks of HFMD have been exclusive to the Asia-Pacific region, it will also be crucial to understand the potential competitive exclusion and invasion dynamics of HFMD and other enteroviral diseases in other parts of the world due to cross-protective immunity from other enteroviruses. Recent concerns with disease due to outbreaks of Enterovirus D68 in the US and other countries only increase the urgency of understanding the spatial and viral dynamics of this group of infections. A more accurate reflection of the sero-epidemiological and antigenic landscape [72] across HFMD and other enteroviral diseases will be a key component of efforts to reduce the clinical burden of HFMD and its associated economic costs in affected areas.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 38

Attachment 39

Attachment 40

Attachment 41

Attachment 42

Attachment 43

Attachment 44

Attachment 45