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

Background 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. The disease has had a substantial burden throughout East and Southeast Asia over the past 15 y. China reported 9 million cases of HFMD between 2008 and 2013, with the two serotypes Enterovirus A71 (EV-A71) and Coxsackievirus A16 (CV-A16) being responsible for the majority of these cases. Three recent phase 3 clinical trials showed that inactivated monovalent EV-A71 vaccines manufactured in China were highly efficacious against HFMD associated with EV-A71, but offered no protection against HFMD caused by CV-A16. To better inform vaccination policy, we used mathematical models to evaluate the effect of prospective vaccination against EV-A71-associated HFMD and the potential risk of serotype replacement by CV-A16. We also extended the model to address the co-circulation, and implications for vaccination, of additional non-EV-A71, non-CV-A16 serotypes of enterovirus. Methods and Findings Weekly reports of HFMD incidence from 31 provinces in Mainland China from 1 January 2009 to 31 December 2013 were used to fit multi-serotype time series susceptible–infected–recovered (TSIR) epidemic models. We obtained good model fit for the two-serotype TSIR with cross-protection, capturing the seasonality and geographic heterogeneity of province-level transmission, with strong correlation between the observed and simulated epidemic series. The national estimate of the basic reproduction number, R 0, weighted by provincial population size, was 26.63 for EV-A71 (interquartile range [IQR]: 23.14, 30.40) and 27.13 for CV-A16 (IQR: 23.15, 31.34), with considerable variation between provinces (however, predictions about the overall impact of vaccination were robust to this variation). EV-A71 incidence was projected to decrease monotonically with higher coverage rates of EV-A71 vaccination. Across provinces, CV-A16 incidence in the post-EV-A71-vaccination period remained either comparable to or only slightly increased from levels prior to vaccination. The duration and strength of cross-protection following infection with EV-A71 or CV-A16 was estimated to be 9.95 wk (95% confidence interval [CI]: 3.31, 23.40) in 68% of the population (95% CI: 37%, 96%). Our predictions are limited by the necessarily short and under-sampled time series and the possible circulation of unidentified serotypes, but, nonetheless, sensitivity analyses indicate that our results are robust in predicting that the vaccine should drastically reduce incidence of EV-A71 without a substantial competitive release of CV-A16. Conclusions The ability of our models to capture the observed epidemic cycles suggests that herd immunity is driving the epidemic dynamics caused by the multiple serotypes of enterovirus. Our results predict that the EV-A71 and CV-A16 serotypes provide a temporary immunizing effect against each other. Achieving high coverage rates of EV-A71 vaccination would be necessary to eliminate the ongoing transmission of EV-A71, but serotype replacement by CV-A16 following EV-A71 vaccination is likely to be transient and minor compared to the corresponding reduction in the burden of EV-A71-associated HFMD. Therefore, a mass EV-A71 vaccination program of infants and young children should provide significant benefits in terms of a reduction in overall HFMD burden.


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][21][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][27][28]. There is also some historical and more recent evidence of interference between polioviruses (Enterovirus C species) and other enteroviruses [29][30][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][37][38][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][41][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.

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), I ðrÞ t;i 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 6 ¼ 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 χ 2 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 loglinear 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 2 ) 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 6 ¼ 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 χ 2 distribution with 1 degree of freedom. We estimated two time-varying epidemiological parameters for each serotype: R 0 (basic reproductive number) and R E (effective reproductive number), obtained by rearranging Eq 7 and defined as: R 0 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.

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 0 , 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 logtransformed 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 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 0 (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 Table 1. Median and interquartile range of R 0 by serotype and by province. regression in Eq 8 and was found to be generally low, between-week variation by province was assessed by the coefficient of variation inb s , and between-province variation was found to be generally high for both EV-A71 and CV-A16 (S1 Table;

Impact of EV-A71 Vaccination
Using our estimate of R 0 of approximately 25 and the basic equation for calculating the critical vaccination threshold p c = 1 − 1/R 0 [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 δ is the proportion of those infected with one serotype who are cross-protected against infection from all other serotypes for a fixed duration of k time steps, and the product k Á δ represents the average duration of cross-protection. 95% CIs for the individual parameters are derived from the profile likelihood using the χ 2 distribution with 1 degree of freedom. three types of EV-A71-containing vaccine on future EV-A71 and CV-A16 incidence (assuming R 0 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

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][57][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 b (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.

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 0 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 0 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][41][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 0 .
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][64][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][69][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. Calculated with α = 0.95 and province-specific maximum likelihood estimates of cross-protection, with provinces ordered by degree of cross-protection (high to low). (B) Maximum transient increase in yearly incidence of CV-A16 (deterministic simulation) by province (y-axis) following vaccine initiation compared to pre-vaccination equilibria, ignoring seasonality in β s , as a function of narrow EV-A71 vaccine coverage (x-axis). Calculated with α = 0.95 and province-specific maximum likelihood estimates of cross-protection, with provinces ordered by degree of cross-protection (high to low). Note that there is no change in the yearly incidence of CV-A16 after vaccination in provinces where the maximum likelihood estimates of cross-protection are zero. . (E and F) Output from deterministic simulation of incidence (y-axis) of (E) EV-A71 and (F) CV-A16 by year (x-axis) for 20 y before to 30 y following vaccine initiation (dotted black line, at year 0), normalized by serotype-specific yearly incidence in year −20 and ignoring seasonality in β s . Vaccination assumed to be narrow monovalent EV-A71 vaccine (administered at birth) achieving 90% coverage. . (E and F) Output from deterministic simulation of incidence (y-axis) of (E) EV-A71 and (F) CV-A16 by year (x-axis) for 20 y before to 30 y following vaccine initiation (dotted black line, at year 0), normalized by serotype-specific yearly incidence in year −20 and ignoring seasonality in β s . Vaccination assumed to be narrow monovalent EV-A71 vaccine (administered at birth) achieving 90% coverage.

Supporting Information
S for EV-A71 = 0.087 and S for CV-A16 = 0.044. Calculated with province-specific maximum likelihood estimates of cross-protection (k = 21 wk and δ = 1). (TIFF)  Table. Two-serotype TSIR model estimates for EV-A71 and CV-A16. Mean proportion of individuals that are susceptible to EV-A71 and CV-A16 ( S), reporting rate of EV-A71 and CV-A16 (ρ), mean weekly estimated transmission rate of EV-A71 and CV-A16 ( b), and coefficient of variation (CV) in estimated transmission rate of EV-A71 and CV-A16 by province. Calculated with α = 0.95 and province-specific maximum likelihood estimates of cross-protection. The 95% CIs for S are derived from the profile likelihood using the χ 2 distribution with 1 degree of freedom; the 95% CIs for ρ are derived from the standard errors for the coefficient 1/ρ in the OLS regression of cumulative births and cumulative cases; the maximum standard errors (se) of the log(β s ) terms across weeklyb s values are from the log-linear regression, representing variability around estimated β s . (DOCX) S2 Table. Spatial and temporal variation in R 0 of EV-A71 across provinces.b s and 95% CIs on each weekly β s value (row) for EV-A71 for each province (column) from the two-serotype model with α = 0.95 and province-specific maximum likelihood estimates of cross-protection. (DOCX) S3 Table. Spatial and temporal variation in R 0 of CV-A16 across provinces.b s and 95% CIs on each weekly β s value (row) for CV-A16 for each province (column) from the two-serotype model with α = 0.95 and province-specific maximum likelihood estimates of cross-protection. (DOCX) S4 Table. Omission of co-infection from the two-serotype model. We omit co-infection from this analysis, as co-infection of a single individual with both serotypes (EV-A71 and CV-A16) is rarer than expected by chance, given a reasonable sample size of EV-A71 and CV-A16 infections in a given year. We determined this using data from two published studies and χ 2 tests. (DOCX) S5 Table. One-serotype TSIR model estimates for EV-A71 and CV-A16. Mean proportion of individuals that are susceptible to EV-A71 and CV-A16 ( S), reporting rate of EV-A71 and CV-A16 (ρ), mean weekly estimated transmission rate of EV-A71 and CV-A16 ( b), and coefficient of variation (CV) in estimated transmission rate of EV-A71 and CV-A16 by province. Calculated with α = 0.95. The 95% CIs for S are derived from the profile likelihood using the χ 2 distribution with 1 degree of freedom; the 95% CIs for ρ are derived from the standard errors for the coefficient 1/ρ in the OLS regression of cumulative births and cumulative cases; the maximum standard errors (se) of the log(β s ) terms across weeklyb s values are from the log-linear regression, representing variability around estimated β s . (DOCX) S6 Table. Three-serotype TSIR model estimates for EV-A71 and CV-A16. Mean proportion of individuals that are susceptible to EV-A71 and CV-A16 ( S), reporting rate of EV-A71 and CV-A16 (ρ), mean weekly estimated transmission rate of EV-A71 and CV-A16 ( b), and coefficient of variation (CV) in estimated transmission rate of EV-A71 and CV-A16 by province. Calculated with α = 0.95 and province-specific maximum likelihood estimates of cross-protection from the two-serotype model. The 95% CIs for S are derived from the profile likelihood using the χ 2 distribution with 1 degree of freedom; the 95% CIs for ρ are derived from the standard errors for the coefficient 1/ρ in the OLS regression of cumulative births and cumulative cases; the maximum standard errors (se) of the log(β s ) terms across weeklyb s values are from the log-linear regression, representing variability around estimated β s . (DOCX) S7 Table. Within-province estimates of mean absolute yearly prediction error between observed and simulated weekly incidence. Shown for both EV-A71 and CV-A16 from 2010 to 2013, for the one-, two-, and three-serotype models with α = 0.95 and province-specific maximum likelihood estimates of cross-protection from the two-serotype model. Values with asterisks indicate lower PE for the serotype in the two-serotype model compared to the oneserotype model, and values with carets indicate lower PE for the serotype in the three-serotype model compared to the two-serotype model. (DOCX) S8 Table. Within-province estimates of the coefficient of determination (R 2 ) between observed and simulated weekly incidence. Shown for both EV-A71 and CV-A16 from 2010 to 2013, for the one-, two-, and three-serotype models with α = 0.95 and province-specific maximum likelihood estimates of cross-protection from the two-serotype model. (DOCX)

Editors' Summary
Background Hand, foot, and mouth disease (HFMD)-a common ailment that mainly affects young children-is caused by a group of enteroviruses (Enterovirus A species), which are close relatives of polioviruses (Enterovirus C species). Enteroviruses are divided into various viral serotypes (variants defined by molecules on their surface that are recognized by the immune system), and the most common serotypes that cause HFMD are Enterovirus A71 (EV-A71) and Coxsackievirus A16 (CV-A16). Enteroviruses spread from person to person through contact with the mucus or saliva produced when an infected individual coughs or sneezes, with the feces or the fluid from vesicles of an infected individual, and through contact with contaminated surfaces. Frequent handwashing and good hygiene practices can reduce the spread of HFMD. Symptoms of HFMD include fever, sore throat, a nonitchy red rash with small blisters on the hands and feet, and painful mouth ulcers. HFMD is usually a self-limiting illness, and most children recover within 7-10 days. A small proportion of patients infected with EV-A71 develop severe complications such as meningitis (infection of the membranes around the brain and spinal cord) or encephalitis (infection of the brain). Currently, there is no specific treatment for HFMD, and vaccines are not yet available for routine use.

Why Was This Study Done?
HFMD is increasingly common in East and Southeast Asia. China, for example, reported 9 million cases of HFMD between 2008 and 2013. Vaccination is a specific and effective way to reduce the burden of HFMD in China. In three clinical trials, inactivated monovalent EV-A71 vaccines made in China were highly efficacious against EV-A71-associated HFMD but provided no protection against CV-A16-associated HFMD (an inactivated monovalent vaccine contains a single virus strain that is unable to replicate; exposure to the vaccine "primes" the immune system to respond quickly when challenged with live virus, thereby preventing infection with that virus). So, before these vaccines can be used for routine vaccination of infants, it is important to know whether vaccination with EV-A71 will alter the burden of HFMD caused by other enterovirus serotypes. In particular, it is important to know whether infection with EV-A71 provides cross-protection against CV-A16 and whether infections with CV-A16 or other serotypes might increase following vaccination against EV-A71 because of increased circulation of these viruses in the population (serotype replacement). Here, the researchers use mathematical models to assess the effect of vaccination against EV-A71-associated HFMD in China and the potential risk of serotype replacement by CV-A16.

What Did the Researchers Do and Find?
The researchers used weekly data on HFMD incidence collected in 31 Chinese provinces between 2009 and 2013 to develop a two-serotype time series susceptible-infected-recovered epidemic model (a model in which individuals are born, become susceptible to a disease, become infected and infectious with the disease, and recover). The model accurately simulated the epidemic cycles of EV-A71-and CV-A16-associated HFMD for the 31 provinces and the seasonal transmission patterns in both northern and southern Chinese provinces. It provided an estimate of cross-protection following infection with EV-A71 or CV-A16 of ten weeks in 68% of the population (an average duration of cross-protection of 6.77 weeks). The estimated basic reproduction number (the average number of additional cases one case of an infectious disease generates in an otherwise uninfected population) across China for both serotypes was 25, which means that vaccination coverage levels of above 96% are required to achieve population-level immunity. Finally, the model predicted a decrease in EV-A71-associated HFMD incidence with higher rates of EV-A71 vaccination and suggested that CV-A16 incidence following EV-A71 vaccination would be comparable to or only slightly higher than its incidence before vaccination.

What Do These Findings Mean?
These findings suggest that herd immunity (indirect protection from infectious disease that occurs when most of a population has become immune to an infection, thereby providing some protection for the rare individuals who are not immune) is driving the dynamics of the HFMD epidemic caused by multiple enterovirus serotypes in China. Moreover, they suggest that the infection with EV-A71 or CV-A16 serotype can provide temporary immunity against each the other serotype and that serotype replacement by CV-A16 following EV-A71 vaccination is likely to be transient and minor compared to the reduction in the burden of EV-A71-associated HFMD produced by vaccination. The accuracy of these findings depends on the assumptions included in the model and the quality and quantity of data used to run the models. However, the researchers suggest that a mass EV-A71 vaccination campaign targeted at infants and young children should greatly reduce the burden of HFMD in China, provided a high vaccination uptake level is achieved.

Additional Information
This list of resources contains links that can be accessed when viewing the PDF on a device or via the online version of the article at http://dx.doi.org/10.1371/journal.pmed.1001958.
• The US Centers for Disease Control and Prevention provide information on hand, foot, and mouth disease (in English and Spanish), including a podcast on the condition • The UK National Health Service Choices website provides detailed information on hand, foot, and mouth disease • Further information about hand, foot, and mouth disease is provided by the World Health Organization (including up-to-date HFMD surveillance reports from China), the Nemours Foundation (in English and Spanish), and MedlinePlus (in English and Spanish) • The Government of the Hong Kong Special Administrative Region Department of Health Centre for Health Protection provides information on hand, foot, and mouth disease