Methods
The severity of infection was classified into tiers which require different levels of care: (i) clinically apparent, or symptomatic, cases; (ii) hospitalisations; (iii) complications, or severe cases when the child displays neurological or cardiopulmonary complications; and (iv) fatal cases. These manifestations were matched against similar definitions from the Global Burden of Disease 2013 (GBD2013) study16 that have existing estimates of disability weights, and the overall disease burden is the average of each level of severity, weighted by the proportion of cases entering into each level.
These severity tiers were arranged as a pyramid where lower echelons represent the majority of (milder) infections with decreasing frequency but increasing severity towards the apex11 17 (figure 1). The parameters which describe the severity of HFMD are the proportions of individuals in one state who reach the next higher severity level. This framework implicitly assumes that cases with higher severity are a proportion of cases at the preceding level. For example, those with significant CNS complications are modelled to be a subset of hospitalised cases, and similarly, all hospitalised cases are a subset of notified cases. We built two parallel pyramids for total HFMD (of all aetiologies) and specifically for EV-A71 due to the abundance of data on EV-A71 compared with other viruses, as well as this virus being thought to be responsible for more severe symptoms.
Figure 1The severity pyramid structure for hand, foot and mouth disease (HFMD) (total, all aetiologies) and EV-A71 in particular. Bordered boxes are data sources, while boxes within the pyramids are estimated quantities. Pale yellow boxes represent the models and outcomes of submodel 1, the infection rate model. This model uses serological data and differential equations to estimate the rate at which infection must occur to match with observed seropositivity level of EV-A71. Blue boxes represent the data, model and outcomes of submodel 2, the symptomatic incidence rate model. This model uses surveillance data from Singapore to infer the incidence rate of HFMD and EV-A71-specific HFMD. Together with submodel 1, the estimated symptomatic rate of EVA71 (cyan box) can be calculated. The orange boxes represent the data and models of submodel 3: the hierarchical model, which is used to estimate indices that quantifies the higher-level severities of HFMD. CNS, central nervous system; MOH, Ministry of Health.
Quantifying the disease burden of HFMD requires an estimate of the total number of cases, which was calculated from the infection rate (both symptomatic and asymptomatic) obtained from seroepidemiological studies. From a systematic review previously published by the authors,14 we extracted 10 studies of EV-A71 serology18–28 (online supplementary appendix A1, A2 and A3) from which the seropositive rate of EV-A71 by age was found using submodel 1 to model the required infection rate to achieve the observed seroprevalence.
To estimate clinically apparent cases, we used data on notified cases from Singapore, where HFMD is a legally notifiable disease and where active surveillance in preschools is conducted via daily temperature taking; this was used to calculate a lower bound for the overall asymptomatic rate of HFMD, in submodel 2. Indices for higher severity levels of HFMD and EV-A71, which include deaths, complications, hospitalisations and symptomatic cases, used information from 27 papers5 29–54 and 11 papers29 34 35 39 55–61 respectively from across Asia and were modelled in submodel 3. Finally, the severity indices estimated in the three submodels were fused to collectively estimate the overall disease burden of HFMD in Asia.
Submodel 1: infection rate model
The infection rate model estimates the age-specific proportion of children who must be infected by the EV-A71 virus each year to reach the observed level of seroprevalence. Serological data from China,18–20 Singapore,21 22 Taiwan,23 24 26 Thailand27 and Vietnam28 were extracted from the original papers (figure 2, left), from which it can be observed that seropositivity to EV-A71 starts low at infancy and reaches a plateau around 40%–60% by adolescence, a finding which is corroborated by cord blood seropositivity (indicative of maternal antibody levels) of 55% in Vietnam,28 47%–50% in Taiwan,25 26 38% in Singapore21 and 26% in China.62 To accommodate these patterns, we developed and fit models of the age-specific seroprevalence using seven dynamical models which make varying assumptions for the reason behind these phenomena. These assumptions include (i) an infection rate that drops to zero at age T, beyond which no further infection is possible; (ii) an exponentially decreasing infection rate which is asymptotically zero as age increases and (iii) antibodies for EV-A71 in humans will diminish over time like influenza63 64 and eventually reach an equilibrium from new infections and loss of immunity. The details of each of the models are explained further in online supplementary appendix A1.
Figure 2Data and results from the infection rate model (submodel 1). Left: seropositivity levels of EV-A71 with fitted curve. Red data points are from China, blue from Taiwan and orange from Singapore, green from Vietnam and brown from Thailand. The black line is the fitted model with 95% Bayesian credible intervals. Middle: required and observed rate of EV-A71 infection. The brown line shows the age-specific infection rate of EV-A71 implied from the seropositivity curve and is the required rate to get the observed seropositivity levels; the blue line shows the calculated average EV-A71 infection rate based on Singapore’s data from submodel 2. The difference between the two lines are the unobserved cases. Right: asymptomatic infection rate. Calculated from the percentage difference between the required and observed rates of infection in the middle.
Posterior distributions of parameters were estimated in R65 using Markov chain Monte Carlo (MCMC) and R package ‘deSolve’66 with non-informative prior distributions. The optimal model to represent infection rate was determined by the fit to the data and the model complexity and was evaluated using the deviance information criteria (DIC).
Submodel 2: symptomatic incidence rate model
Since October 2000, Singapore has adopted a comprehensive approach for HFMD surveillance including a mandatory case reporting system by medical practitioners and childcare centres,67 with children at childcares being screened once or twice daily for symptoms and then subsequently isolated for 10 days if symptomatic. The clearly defined population catchment in the city-state and the dual layers of doctor-driven and teacher-driven surveillance make the data on notified symptomatic presentations of HFMD collected by the Singapore’s Ministry of Health unusually complete compared with the rest of Asia. In calculating symptomatic incidence rate, we make the assumption that all symptomatic paediatric HFMD were clinically apparent and captured by the surveillance system in Singapore.
We constructed a matrix containing the incidence of cases in year (2003–2012) among age a (1–18), scaled by the number of individuals in the resident population for that age and year derived from the census.68 To separate the effects of age and time on incidence rate, singular value decomposition was performed on the incidence matrix to find the first left-singular vector which corresponds to the largest singular value. This vector is normalised to , which is equivalently the normalised first principal component of , and represents the systematic variation attributable to age. The second vector is a time-series index, , representing the stochastic variation attributable to HFMD incidence level each year. Values for years 2013–2015 were estimated using linear regression using data from Singapore’s Weekly Infectious Disease Bulletin.69 The modelled HFMD incidence rate is estimated from the product of the two vectors, .
The EV-A71-specific incidence is a proportion of all HFMD cases. We obtain , the year effect due to EV-A71, by taking it as a binomial sample from the total year effect, . The proportion of cases attributed to EV-A71 in 2008–2015, , was estimated from data obtained from the National Surveillance Program in Singapore which collects clinical specimens from patients in 107 clinics and 2 hospitals, KK Women’s and Children’s Hospital and National University Hospital (figure 3, right). Throat or vesicle swabs were collected from patients <12 years of age with clinical diagnosis of HFMD. Specimens were stored at room temperature or 4° (preferred) for no longer than 48 hours, before being sent to the National Public Health Laboratory for testing and typing. Data on prior to 2008 were estimated in a study by Ang et al.52 The average EV-A71 incidence rate by age is the average of the EV-A71 year effect, weighted by the age effect, or . The model is described in detail in online supplementary appendix B.
Figure 3Age and year trend of hand, foot and mouth disease (HFMD) symptomatic infection from the Incidence rate model (submodel 2). Left: the age effect. The black line reflects the overall proportion of incidence for each year of age, and red line represents the age effect for individual years from 2005 to 2012. Middle: the year effect. The black line is an index which reflects the number of notified cases for 2005–2012. The year effect for 2013–2015 is estimated from data (notified cases) obtained from the Weekly Infectious Disease Bulletin released by the Ministry of Health, Singapore, using linear regression. The year effect of EV-A71 (blue) is mathematically derived from the proportion of cases due to EV-A71 estimated by laboratory data (right). Right: laboratory data from Singapore’s surveillance system. This shows the proportion of cases due to each virus in each year from 2005 to 2015.
Submodel 3: hierarchical model of disease severity
The upper levels of the severity pyramid from across Asia were modelled as nested binomial distributions, with the probability of reaching state i or higher given the case reaches the lower state j denoted as . The states are D (death), S (severe disease, ‘CNS’), H (hospitalised) and C (symptomatic case) so that, for example, denotes the case hospitalisation rate. Standard rules of conditional probability mean that the product of these probabilities gives the risk of progressing through multiple levels; for instance, the overall case fatality rate (CFR) is . As a result, the likelihood function for data with missing tiers of the pyramid is specified via combinations of the basic parameters. For example, 4625 hospitalisations and 11 deaths were documented in the Sarawak 1997 outbreak,49 with no mention of severe cases. This data set will allow us to estimate , and thus make inferences on and through the multiplicative relationship, where is the individual-level estimate for paper p. Detailed description of the model and all data sources can be found in online supplementary appendix C1 (data), C2 and C3 (results) and C4 (model). The posterior distribution of parameters were estimated using the Metropolis Hastings algorithm (MCMC) in R.
Disability-adjusted life-year burden
To estimate the DALY lost across Asia to HFMD, we conservatively estimated the incidence of HFMD infection to be twice the incidence of EV-A71 using the incidence rate model (submodel 2), to allow for infection by other aetiological agents of HFMD which tend to be more common than EV-A71. The number of symptomatic cases was then calculated using the asymptomatic rate from the infection rate model (submodel 1). The numbers of more severe manifestations (hospitalised, CNS and deaths) were then estimated using the parameters from the hierarchical model (submodel 3). Again, posterior distributions were used to characterise uncertainty in these estimates.
Disability weights from GBD201316 for severe (0.133, 0.088–0.19), moderate (0.051, 0.032–0.074) and mild (0.006, 0.002–0.012) acute episodes of infectious disease were used as estimates for HFMD with complication, hospitalisation and symptomatic HFMD, respectively. To estimate the disability duration of each tier, we extracted commonly reported summary statistics such as the sample size, mean, median, range, CIs and variances of the disability duration of their samples from 31 papers and modelled them into Weibull distributions which best describe the statistical properties (online supplementary appendix D). The disability duration for fatalities was determined using the period life expectancy of Japan,70 which represents the expected lifespan of humans as healthcare improves in the current century.
We used the age-weighted approach15 which assigns greater weights to illnesses that occur during periods when the social role is higher, most notably between ages 10 and 55. As disability from HFMD usually occurs among children, this method provides a conservative estimate of DALY.
Under this framework, the formula for DALY for an individual is
where is the disability weight (1 for premature mortality), =0.03 is the discount rate, =0.166 is the age-weighting correcting coefficient, =0.04 is the parameter from the age-weighting function, is the age of onset and is the duration of disability or time lost due to premature mortality.