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

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.

Published in the journal:
. PLoS Med 13(2): e32767. doi:10.1371/journal.pmed.1001958

Category:
Research Article

doi: 10.1371/journal.pmed.1001958

## 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 χ^{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 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*^{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* ≠ *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.

## 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*_{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 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*_{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 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*_{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 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 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*_{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*–*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*–*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

##### Zdroje

1. Podin Y, Gias EL, Ong F, Leong Y-W, Yee S-F, Yusof MA, et al. Sentinel surveillance for human enterovirus 71 in Sarawak, Malaysia: lessons from the first 7 years. BMC Public Health. 2006;6:180. doi: 10.1186/1471-2458-6-180 16827926

2. Tseng F-C, Huang H-C, Chi C-Y, Lin T-L, Liu C-C, Jian J-W, et al. Epidemiological survey of enterovirus infections occurring in Taiwan between 2000 and 2005: analysis of sentinel physician surveillance data. J Med Virol. 2007;79:1850–1860. doi: 10.1002/jmv.21006 17935170

3. World Health Organization Western Pacific Region. Emerging disease surveillance and response. Hand, foot and mouth disease information sheet. 2012 Jul 10 [cited 20 Aug 2014]. Manila: World Health Organization Regional Office for the Western Pacific. Available: http://www.wpro.who.int/emerging_diseases/hfmd.information.sheet/en/

4. Chan Y-F, Sam I-C, Wee K-L, Abubakar S. Enterovirus 71 in Malaysia: a decade later. Neurol Asia. 2011;16:1–15.

5. Chan LG, Parashar UD, Lye MS, Ong FGL, Zaki SR, Alexander JP, et al. Deaths of children during an outbreak of hand, foot, and mouth disease in Sarawak, Malaysia: clinical and pathological characteristics of the disease. Clin Infect Dis. 2000;31:678–683. doi: 10.1086/314032 11017815

6. Ho M, Chen E-R, Hsu K-H, Twu S-J, Chen K-T, Tsai S-F, et al. An epidemic of enterovirus 71 infection in Taiwan. Taiwan Enterovirus Epidemic Working Group. N Engl J Med. 1999;341:929–935. doi: 10.1056/NEJM199909233411301 10498487

7. Fujimoto T, Chikahira M, Yoshida S, Ebira H, Hasegawa A, Totsuka A, et al. Outbreak of central nervous system disease associated with hand, foot, and mouth disease in Japan during the summer of 2000: detection and molecular epidemiology of enterovirus 71. Microbiol Immunol. 2002;46:621–627. doi: 10.1111/j.1348-0421.2002.tb02743.x 12437029

8. Chan KP, Goh KT, Chong CY, Teo ES, Lau G, Ling AE. Epidemic hand, foot and mouth disease caused by human enterovirus 71, Singapore. Emerg Infect Dis. 2003;9:78–85. doi: 10.3201/eid1301.020112 12533285

9. Van Tu P, Thao NTT, Perera D, Truong KH, Tien NTK, Thuong TC, et al. Epidemiologic and virologic investigation of hand, foot, and mouth disease, southern Vietnam, 2005. Emerg Infect Dis. 2007;13:1733–1741. doi: 10.3201/eid1311.070632 18217559

10. Zhang Y, Tan X-J, Wang H-Y, Yan D-M, Zhu S-L, Wang D-Y, et al. An outbreak of hand, foot, and mouth disease associated with subgenotype C4 of human enterovirus 71 in Shandong, China. J Clin Virol. 2009;44:262–267. doi: 10.1016/j.jcv.2009.02.002 19269888

11. Zhang Y, Zhu Z, Yang W, Ren J, Tan X, Wang Y, et al. An emerging recombinant human enterovirus 71 responsible for the 2008 outbreak of hand foot and mouth disease in Fuyang city of China. Virol J. 2010;7:94. doi: 10.1186/1743-422X-7-94 20459851

12. Ma E, Chan KC, Cheng P, Wong C, Chuang SK. The enterovirus 71 epidemic in 2008—public health implications for Hong Kong. Int J Infect Dis. 2010;14:e775–e780. doi: 10.1016/j.ijid.2010.02.2265 20599410

13. Ryu W-S, Kang B, Hong J, Hwang S, Kim J, Cheon D-S. Clinical and etiological characteristics of enterovirus 71-related diseases during a recent 2-year period in Korea. J Clin Microbiol. 2010;48:2490–2494. doi: 10.1128/JCM.02369-09 20463159

14. Seiff A. Cambodia unravels cause of mystery illness. Lancet. 2012;380:206.

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

16. Feng H, Duan G, Zhang R, Zhang W. Time series analysis of hand-foot-mouth disease hospitalization in Zhengzhou: establishment of forecasting models using climate variables as predictors. PLoS ONE. 2014;9:e87916. doi: 10.1371/journal.pone.0087916 24498221

17. Lee BY, Wateska AR, Bailey RR, Tai JHY, Bacon KM, Smith KJ. Forecasting the economic value of an Enterovirus 71 (EV71) vaccine. Vaccine. 2010;28:7731–7736. doi: 10.1016/j.vaccine.2010.09.065 20923711

18. US Centers for Disease Control and Prevention. Non-polio enterovirus: enterovirus D68. [cited 19 Oct 2014]. Atlanta: US Centers for Disease Control and Prevention. Available: http://www.cdc.gov/non-polio-enterovirus/about/ev-d68.html

19. McMinn PC. An overview of the evolution of enterovirus 71 and its clinical and public health significance. FEMS Microbiol Rev. 2002;26:91–107. 12007645

20. De W, Changwen K, Wei L, Monagin C, Jin Y, Cong M, et al. A large outbreak of hand, foot, and mouth disease caused by EV71 and CAV16 in Guangdong, China, 2009. Arch Virol. 2011;156:945–953. doi: 10.1007/s00705-011-0929-8 21305327

21. Liu MY, Liu W, Luo J, Liu Y, Zhu Y, Berman H, et al. Characterization of an outbreak of hand, foot, and mouth disease in Nanchang, China in 2010. PLoS ONE. 2011;6:e25287. doi: 10.1371/journal.pone.0025287 21980416

22. Ma E, Fung C, Yip SHL, Wong C, Chuang SK, Tsang T. Estimation of the basic reproduction number of enterovirus 71 and coxsackievirus A16 in hand, foot, and mouth disease outbreaks. Pediatr Infect Dis J. 2011;30:675–679. doi: 10.1097/INF.0b013e3182116e95 21326133

23. Li R, Liu L, Mo Z, Wang X, Xia J, Liang Z, et al. An inactivated enterovirus 71 vaccine in healthy children. N Engl J Med. 2014;370:829–837. doi: 10.1056/NEJMoa1303224 24571755

24. China Food and Drug Administration. [Announcement on licensed drugs approved by China Food and Drug Administration (No. 4 in 2016).] 2016 Jan 12 [cited 17 Jan 2016]. Available: http://www.sfda.gov.cn/WS01/CL0087/142000.html.

25. Zhu F-C, Meng F-Y, Li J-X, Li X-L, Mao Q-Y, Tao H, et al. Efficacy, safety, and immunology of an inactivated alum-adjuvant enterovirus 71 vaccine in children in China: a multicentre, randomised, double-blind, placebo-controlled, phase 3 trial. Lancet. 2013;381:2024–2032. doi: 10.1016/S0140-6736(13)61049-1 23726161

26. Pitzer VE, Patel MM, Lopman BA, Viboud C, Parashar UD, Grenfell BT. Modeling rotavirus strain dynamics in developed countries to understand the potential impact of vaccination on genotype distributions. Proc Natl Acad Sci U S A. 2011;108:19353–19358. doi: 10.1073/pnas.1110507108 22084114

27. Kamo M, Sasaki A. The effect of cross-immunity and seasonal forcing in a multi-strain epidemic model. Phys Nonlinear Phenom. 2002;165:228–241. doi: 10.1016/S0167-2789(02)00389-5

28. Zhuang L, Cressie N, Pomeroy L, Janies D. Multi-species SIR models from a dynamical Bayesian perspective. Theor Ecol. 2013;6:457–473. doi: 10.1007/s12080-013-0180-x

29. Hovi T, Roivainen M. Peptide antisera targeted to a conserved sequence in poliovirus capsid VP1 cross-react widely with members of the genus Enterovirus. J Clin Microbiol. 1993;31:1083–1087. 8388885

30. Samuelson A, Forsgren M, Johansson B, Wahren B, Sällberg M. Molecular basis for serological cross-reactivity between enteroviruses. Clin Diagn Lab Immunol. 1994;1:336–341. 7496972

31. Deng C, Yang C, Wan J, Zhu L, Leng Q. Irregular poliovirus vaccination correlates to pulmonary edema of hand, foot, and mouth disease. Clin Vaccine Immunol. 2011;18:1589–1590. doi: 10.1128/CVI.05132-11 21752953

32. Yan X-F, Gao S, Xia J-F, Ye R, Yu H, Long J-E. Epidemic characteristics of hand, foot, and mouth disease in Shanghai from 2009 to 2010: enterovirus 71 subgenotype C4 as the primary causative agent and a high incidence of mixed infections with coxsackievirus A16. Scand J Infect Dis. 2012;44:297–305. doi: 10.3109/00365548.2011.634433 22176514

33. Li Y, Zhu R, Qian Y, Deng J, Sun Y, Liu L, et al. Comparing Enterovirus 71 with Coxsackievirus A16 by analyzing nucleotide sequences and antigenicity of recombinant proteins of VP1s and VP4s. BMC Microbiol. 2011;11:246. doi: 10.1186/1471-2180-11-246 22050722

34. Lin Y, Wen K, Pan Y, Wang Y, Che X, Wang B. Cross-reactivity of anti-EV71 IgM and neutralizing antibody in series sera of patients infected with Enterovirus 71 and Coxsackievirus A 16. J Immunoassay Immunochem. 2011;32:233–243. doi: 10.1080/15321819.2011.559297 21574094

35. Chou A-H, Liu C-C, Chang J-Y, Jiang R, Hsieh Y-C, Tsao A, et al. Formalin-inactivated EV71 vaccine candidate induced cross-neutralizing antibody against subgenotypes B1, B4, B5 and C4A in adult volunteers. PLoS ONE. 2013;8:e79783. doi: 10.1371/journal.pone.0079783 24278177

36. Li Y, Zhang J, Zhang X. Modeling and preventive measures of hand, foot and mouth disease (HFMD) in China. Int J Environ Res Public Health. 2014;11:3108–3117. doi: 10.3390/ijerph110303108 24633146

37. Roy N. Mathematical modeling of hand-foot-mouth disease: quarantine as a control measure. Int J Adv Sci Eng Technol Res. 2012;1:34–44.

38. Zhu Y, Xu B, Lian X, Lin W, Zhou Z, Wang W. A hand-foot-and-mouth disease model with periodic transmission rate in Wenzhou, China. Abstr Appl Anal. 2014;2014:e234509. doi: 10.1155/2014/234509

39. Wang Y, Feng Z, Yang Y, Self S, Gao Y, Longini IM, et al. Hand, foot, and mouth disease in China: patterns of spread and transmissibility. Epidemiol Camb Mass. 2011;22:781–792. doi: 10.1097/EDE.0b013e318231d67a

40. Weinberger DM, Malley R, Lipsitch M. Serotype replacement in disease after pneumococcal vaccination. Lancet. 2011;378:1962–1973. doi: 10.1016/S0140-6736(10)62225-8 21492929

41. Hanage WP, Finkelstein JA, Huang SS, Pelton SI, Stevenson AE, Kleinman K, et al. Evidence that pneumococcal serotype replacement in Massachusetts following conjugate vaccination is now complete. Epidemics. 2010;2:80–84. doi: 10.1016/j.epidem.2010.03.005 21031138

42. Lloyd-Smith JO. Vacated niches, competitive release and the community ecology of pathogen eradication. Philos Trans R Soc Lond B Biol Sci. 2013;368. doi: 10.1098/rstb.2012.0150

43. Metcalf CJE, Bjørnstad ON, Grenfell BT, Andreasen V. Seasonality and comparative dynamics of six childhood infections in pre-vaccination Copenhagen. Proc Biol Sci. 2009;276:4111–4118. doi: 10.1098/rspb.2009.1058 19740885

44. National Bureau of Statistics of China. Annual data. [cited 6 Jul 2015]. Available: http://www.stats.gov.cn/english/Statisticaldata/AnnualData/

45. Finkenstädt BF, Bjørnstad ON, Grenfell BT. A stochastic model for extinction and recurrence of epidemics: estimation and inference for measles outbreaks. Biostat Oxf Engl. 2002;3:493–510. doi: 10.1093/biostatistics/3.4.493

46. Reich NG, Shrestha S, King AA, Rohani P, Lessler J, Kalayanarooj S, et al. Interactions between serotypes of dengue highlight epidemiological impact of cross-immunity. J R Soc Interface. 2013;10:20130414. doi: 10.1098/rsif.2013.0414 23825116

47. Zhu F-C, Liang Z-L, Meng F-Y, Zeng Y, Mao Q-Y, Chu K, et al. Retrospective study of the incidence of HFMD and seroepidemiology of antibodies against EV71 and CoxA16 in prenatal women and their infants. PLoS ONE. 2012;7:e37206. doi: 10.1371/journal.pone.0037206 22662137

48. Finkenstädt BF, Grenfell BT. Time series modelling of childhood diseases: a dynamical systems approach. J R Stat Soc Ser C Appl Stat. 2000;49:187–205.

49. Bjørnstad ON, Finkenstädt BF, Grenfell BT. Dynamics of measles epidemics: estimating scaling of transmission rates using a time series SIR model. Ecol Monogr. 2002;72:169–184.

50. Glass K, Xia Y, Grenfell BT. Interpreting time-series analyses for continuous-time biological models—measles as a case study. J Theor Biol. 2003;223:19–25. 12782113

51. Bjørnstad ON, Finkenstädt BF, Grenfell BT. Dynamics of measles epidemics: scaling noise, determinism, and predictability with the TSIR model. Ecol Monogr. 2002;72:185–202.

52. Zhu F, Xu W, Xia J, Liang Z, Liu Y, Zhang X, et al. Efficacy, safety, and immunogenicity of an enterovirus 71 vaccine in China. N Engl J Med. 2014;370:818–828. doi: 10.1056/NEJMoa1304923 24571754

53. Grassly NC, Fraser C. Seasonal infectious disease epidemiology. Proc R Soc Lond B Biol Sci. 2006;273:2541–2550. doi: 10.1098/rspb.2006.3604

54. Bartlett MS. The critical community size for measles in the United States. J R Stat Soc Ser A. 1960;123:37–44. doi: 10.2307/2343186

55. Keeling MJ, Rohani P. Modeling infectious diseases in humans and animals. Princeton (New Jersey): Princeton University Press.

56. He Y-Q, Chen L, Xu W-B, Yang H, Wang H-Z, Zong W-P, et al. Emergence, circulation, and spatiotemporal phylogenetic analysis of coxsackievirus A6- and coxsackievirus A10-associated hand, foot, and mouth disease infections from 2008 to 2012 in Shenzhen, China. J Clin Microbiol. 2013;51:3560–3566. doi: 10.1128/JCM.01231-13 23966496

57. Hongyan G, Chengjie M, Qiaozhi Y, Wenhao H, Juan L, Lin P, et al. Hand, foot and mouth disease caused by coxsackievirus A6, Beijing, 2013. Pediatr Infect Dis J. 2014;33:1302–1303. doi: 10.1097/INF.0000000000000467 25037037

58. Tian H, Zhang Y, Sun Q, Zhu S, Li X, Pan Z, et al. Prevalence of multiple enteroviruses associated with hand, foot, and mouth disease in Shijiazhuang City, Hebei province, China: outbreaks of coxsackieviruses A10 and B3. PLoS ONE. 2014;9:e84233. doi: 10.1371/journal.pone.0084233 24392117

59. Di B, Zhang Y, Xie H, Li X, Chen C, Ding P, et al. Circulation of Coxsackievirus A6 in hand-foot-mouth disease in Guangzhou, 2010–2012. Virol J. 2014;11:157. doi: 10.1186/1743-422X-11-157 25178398

60. Stewart CL, Chu EY, Introcaso CE, Schaffer A, James WD. Coxsackievirus A6-induced hand-foot-mouth disease. JAMA Dermatol. 2013;149:1419–1421. doi: 10.1001/jamadermatol.2013.6777 24172861

61. Bible JM, Pantelidis P, Chan PKS, Tong CYW. Genetic evolution of enterovirus 71: epidemiological and pathological implications. Rev Med Virol. 2007;17:371–379. doi: 10.1002/rmv.538 17487831

62. Huang S-W, Hsu Y-W, Smith DJ, Kiang D, Tsai H-P, Lin K-H, et al. Reemergence of enterovirus 71 in 2008 in Taiwan: dynamics of genetic and antigenic evolution from 1998 to 2008. J Clin Microbiol. 2009;47:3653–3662. doi: 10.1128/JCM.00630-09 19776232

63. Zhang H, An D, Liu W, Mao Q, Jin J, Xu L, et al. Analysis of cross-reactive neutralizing antibodies in human HFMD serum with an EV71 pseudovirus-based assay. PLoS ONE. 2014;9:e100545. doi: 10.1371/journal.pone.0100545 24964084

64. Mao Q, Cheng T, Zhu F, Li J, Wang Y, Li Y, et al. The cross-neutralizing activity of enterovirus 71 subgenotype c4 vaccines in healthy Chinese infants and children. PLoS ONE. 2013;8:e79599. doi: 10.1371/journal.pone.0079599 24260259

65. Huang M-L, Chiang P-S, Chia M-Y, Luo S-T, Chang L-Y, Lin T-Y, et al. Cross-reactive neutralizing antibody responses to enterovirus 71 infections in young children: implications for vaccine development. PLoS Negl Trop Dis. 2013;7:e2067. doi: 10.1371/journal.pntd.0002067 23459633

66. Dushoff J, Plotkin JB, Levin SA, Earn DJD. Dynamical resonance can account for seasonality of influenza epidemics. Proc Natl Acad Sci U S A. 2004;101:16915–16916. doi: 10.1073/pnas.0407293101 15557003

67. Grassly NC, Fraser C, Garnett GP. Host immunity and synchronized epidemics of syphilis across the United States. Nature. 2005;433:417–421. doi: 10.1038/nature03072 15674292

68. Xu M, Su L, Cao L, Zhong H, Dong N, Xu J. Enterovirus genotypes causing hand foot and mouth disease in Shanghai, China: a molecular epidemiological analysis. BMC Infect Dis. 2013;13:489. doi: 10.1186/1471-2334-13-489 24148902

69. Wei S-H, Huang Y-P, Liu M-C, Tsou T-P, Lin H-C, Lin T-L, et al. An outbreak of coxsackievirus A6 hand, foot, and mouth disease associated with onychomadesis in Taiwan, 2010. BMC Infect Dis. 2011;11:346. doi: 10.1186/1471-2334-11-346 22168544

70. Lu Q-B, Zhang X-A, Wo Y, Xu H-M, Li X-J, Wang X-J, et al. Circulation of Coxsackievirus A10 and A6 in hand-foot-mouth disease in China, 2009–2011. PLoS ONE. 2012;7:e52073. doi: 10.1371/journal.pone.0052073 23272213

71. Xia Y, Bjørnstad ON, Grenfell BT. Measles metapopulation dynamics: a gravity model for epidemiological coupling and dynamics. Am Nat. 2004;164:267–281. doi: 10.1086/422341 15278849

72. Smith DJ, Lapedes AS, de Jong JC, Bestebroer TM, Rimmelzwaan GF, Osterhaus ADME, et al. Mapping the antigenic and genetic evolution of influenza virus. Science. 2004;305:371–376. doi: 10.1126/science.1097211 15218094

##### Štítky

Interní lékařstvíČlánek vyšel v časopise

### PLOS Medicine

2016 Číslo 2

Nejčtenější v tomto čísle

Tomuto tématu se dále věnují…

- Hand, Foot, and Mouth Disease in China: Modeling Epidemic Dynamics of Enterovirus Serotypes and Implications for Vaccination
- Estimated Effects of Different Alcohol Taxation and Price Policies on Health Inequalities: A Mathematical Modelling Study
- Transforming Living Kidney Donation with a Comprehensive Strategy
- 2015 Reviewer Thank You